Demography, R, and Mac
Come meet me at EPC to discuss Spatial Demography!
Thursday, 7 June 2018, 9:00-10:30 Session 1 Motherhood, Labor Market and Wages: Contextual Determinants of Childbearing in Spain: A Spatial Panel Study
Friday, 8 June 2018, 11:00-12:30 Session 65 Fertility Trends and Prospects: The Geography of Fertility Rates in Low and Middle-Income Countries: Analysis of Cross-Sectional Surveys from 74 Countries
You can download session 9 files for constructing the population pyramids of Georgia here: https://github.com/rladies/meetup-presentations_tbilisi and specify your working directory with setwd(“/Users/mydomain/myfolder/”)
#set working directory
mypath<-"/Users/DrSpengler/The rectification of the Vuldrini/"
#upload shape files
georgia <- readOGR("./GEO_adm/","GEO_adm0")
## OGR data source with driver: ESRI Shapefile ## Source: "./GEO_adm/", layer: "GEO_adm0" ## with 1 features ## It has 70 fields
# plot(georgia, lwd=1.5)
georgia1 <- readOGR("./GEO_adm/","GEO_adm1")
## OGR data source with driver: ESRI Shapefile ## Source: "./GEO_adm/", layer: "GEO_adm1" ## with 12 features ## It has 16 fields
# plot(georgia1)
georgia2 <- readOGR("./GEO_adm/","GEO_adm2")
## OGR data source with driver: ESRI Shapefile ## Source: "./GEO_adm/", layer: "GEO_adm2" ## with 69 features ## It has 18 fields
# plot(georgia2)
gwat <- readOGR("./GEO_wat/" , "GEO_water_lines_dcw")
## OGR data source with driver: ESRI Shapefile ## Source: "./GEO_wat/", layer: "GEO_water_lines_dcw" ## with 559 features ## It has 5 fields
# plot(gwat)
gpop <- raster("./GEO_pop/geo_pop.grd")
# plot(gpop)
galt <- raster("./GEO_msk_alt/GEO_msk_alt.grd")
# plot(galt)
plot(georgia, lwd=1.5) #n1
plot(georgia1, lwd=1.5) #n2
plot(georgia2, lwd=1.5) #n3
plot(georgia, lwd=1.5) #n4 plot(gwat, lwd=1.5, col="blue", add=T) #n4
plot(gpop) #n5 plot(georgia, lwd=1.5, add=T) #n5
plot(galt, lwd=1.5) #n6
tur <- readOGR("./TUR_adm" , "TUR_adm0")
## OGR data source with driver: ESRI Shapefile ## Source: "./TUR_adm", layer: "TUR_adm0" ## with 1 features ## It has 70 fields ## Integer64 fields read as strings: ID_0 OBJECTID_1
arm <- readOGR("./ARM_adm" , "ARM_adm0")
## OGR data source with driver: ESRI Shapefile ## Source: "./ARM_adm", layer: "ARM_adm0" ## with 1 features ## It has 70 fields ## Integer64 fields read as strings: ID_0 OBJECTID_1
rus <- readOGR("./RUS_adm" , "RUS_adm0")
## OGR data source with driver: ESRI Shapefile ## Source: "./RUS_adm", layer: "RUS_adm0" ## with 1 features ## It has 70 fields ## Integer64 fields read as strings: ID_0 OBJECTID_1
aze <- readOGR("./AZE_adm" , "AZE_adm0")
## OGR data source with driver: ESRI Shapefile ## Source: "./AZE_adm", layer: "AZE_adm0" ## with 1 features ## It has 70 fields ## Integer64 fields read as strings: ID_0 OBJECTID_1
plot(georgia, lwd=1.5, col="white", bg="lightblue") plot(georgia1, add=T, lty=2) plot(tur, add=T, col="white") plot(arm, add=T, col="white") plot(rus, add=T, col="white") plot(aze, add=T, col="white")
x.loc <- c(44.32002, 46.35746, 44.40421, 42.18156, 40.71662)
y.loc <- c(43.42472, 40.87209, 40.82228, 40.90945, 41.99276)
nb.lab <- c("Russia", "Azerbaijan", "Armenia", "Turkey", "Black Sea")
plot(georgia, lwd=1.5, col="white", bg="lightblue")
plot(georgia1, add=T, lty=2)
plot(tur, add=T, col="white")
plot(arm, add=T, col="white")
plot(rus, add=T, col="white")
plot(aze, add=T, col="white")
text(x.loc, y.loc, nb.lab)
plot(gwat, col="blue") # plot(georgia1[1,], lwd=1, col="lightblue", border="black", add=T) plot(georgia2, lwd=0.5, border="black", lty=3, add=T) plot(georgia1, border="black", lty=2, add=T) plot(georgia, lwd=1.5, add=T)
head(georgia1)
## ID_0 ISO NAME_0 ID_1 NAME_1 VARNAME_1 NL_NAME_1 HASC_1 CC_1 ## 0 81 GEO Georgia 1034 Abkhazia Sokhumi GE.AB ## 1 81 GEO Georgia 1035 Ajaria Batumi GE.AJ ## 2 81 GEO Georgia 1036 Guria Ozurgeti GE.GU ## 3 81 GEO Georgia 1037 Imereti Kutaisi GE.IM ## 4 81 GEO Georgia 1038 Kakheti Telavi GE.KA ## 5 81 GEO Georgia 1039 Kvemo Kartli Rustavi GE.KK ## TYPE_1 ENGTYPE_1 VALIDFR_1 VALIDTO_1 REMARKS_1 ## 0 Avtonomiuri Respublika Autonomous Republic 1994 Present ## 1 Avtonomiuri Respublika Autonomous Republic 1994 Present ## 2 Region Region 1994 Present ## 3 Region Region 1994 Present ## 4 Region Region 1994 Present ## 5 Region Region 1994 Present ## Shape_Leng Shape_Area ## 0 6.643211 0.9744622 ## 1 3.055014 0.3074264 ## 2 2.880653 0.2092665 ## 3 4.214567 0.6783179 ## 4 6.820519 1.2485036 ## 5 5.219352 0.6807876
coords2<- coordinates(georgia2[2:6,]) admin2 <- c(as.character(georgia2$NAME_2[1:5])) admin2
## [1] "Gagra" "Gali" "Gudauta" "Gulripshi" "Ochamchire"
dt <- read.csv("/Users/ac1y15/Google Drive/blog/RLadies_Georgia_files/Session_3/Data_Extract_From_Subnational_Malnutrition/3f075abc-c51c-40c5-afb1-f8fbcfa30f23_Data.csv", header=T)
dt.1 <- subset(dt, dt$type==1&dt$select==1)
head(dt.1)
## Admin.Region.Name select order ## 6 1 1 ## 7 Georgia, Adjara Aut. Rep. 1 2 ## 16 Georgia, Guria 1 3 ## 26 Georgia, Imereti 1 4 ## 31 Georgia, Kakheti 1 5 ## 36 Georgia, Kvemo Kartli 1 6 ## Admin.Region.Code type ## 6 1 ## 7 GEO_Adjara_Aut._Rep._GE.AR_1297_GEO002 1 ## 16 GEO_Guria_GE.GU_1298_GEO003 1 ## 26 GEO_Imereti_GE.IM_1299_GEO004 1 ## 31 GEO_Kakheti_GE.KA_1300_GEO005 1 ## 36 GEO_Kvemo_Kartli_GE.KK_1301_GEO006 1 ## Series.Name ## 6 ## 7 Prevalence of overweight, weight for height (% of children under 5) ## 16 Prevalence of overweight, weight for height (% of children under 5) ## 26 Prevalence of overweight, weight for height (% of children under 5) ## 31 Prevalence of overweight, weight for height (% of children under 5) ## 36 Prevalence of overweight, weight for height (% of children under 5) ## Series.Code YR2000 YR2005 YR2009 ## 6 NA NA NA ## 7 SN.SH.STA.OWGH.ZS NA 28.1 NA ## 16 SN.SH.STA.OWGH.ZS NA 7.9 NA ## 26 SN.SH.STA.OWGH.ZS 9.9 21.5 NA ## 31 SN.SH.STA.OWGH.ZS 7.0 19.6 13.2 ## 36 SN.SH.STA.OWGH.ZS 9.5 28.2 19.1
library(classInt) nclassint <- 3 #number of colors to be used in the palette cat <- classIntervals(dt.1$YR2005, nclassint,style = "quantile") #style refers to how the breaks are created
colpal <- brewer.pal(nclassint,"Greens") #sequential
color.palette <- findColours(cat,colpal)
is.na(color.palette)
## [1] TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE ## [12] FALSE
bins <- cat$brks
lb <- length(bins)
color.palette[c(1, 10)] <- "gray"
value.vec <- c(round(bins[-length(bins)],2))
value.vec.tail <- c(round(bins[-1],2))
plot(georgia1, col=color.palette, border=T, main="Prevalence of overweight, \nweight for height (% of children under 5)")
legend("topright",fill=c("gray", "#E5F5E0", "#A1D99B", "#31A354"),legend=c("NA",paste(value.vec,":",value.vec.tail)),cex=1.1, bg="white", bty = "n")
# map.scale(41, 41, 2, "km", 2, 100)
map.scale(x=40.1, y=41.2, relwidth=0.1 , metric=T, ratio=F, cex=0.8)
SpatialPolygonsRescale(layout.north.arrow(2), offset= c(40.1, 41.6), scale = 0.5, plot.grid=F)
Understanding what kind of data you have (polygons or points?) and what you want to map is pivotal to start your mapping.
setwd(paste(mypath))
fr.prj <- readOGR(".", "FRA_adm2")
## OGR data source with driver: ESRI Shapefile ## Source: ".", layer: "FRA_adm2" ## with 96 features ## It has 18 fields
## NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
map(fr.prj)
## Warning in SpatialPolygons2map(database, namefield = namefield): database ## does not (uniquely) contain the field 'name'.
head(fr.prj@data)
## ID_0 ISO NAME_0 ID_1 NAME_1 ID_2 NAME_2 VARNAME_2 ## 0 76 FRA France 989 Alsace 13755 Bas-Rhin Unterelsaá ## 1 76 FRA France 989 Alsace 13756 Haut-Rhin Oberelsaá ## 2 76 FRA France 990 Aquitaine 13757 Dordogne <NA> ## 3 76 FRA France 990 Aquitaine 13758 Gironde Bec-D'Ambes ## 4 76 FRA France 990 Aquitaine 13759 Landes Landas ## 5 76 FRA France 990 Aquitaine 13760 Lot-Et-Garonne <NA> ## NL_NAME_2 HASC_2 CC_2 TYPE_2 ENGTYPE_2 VALIDFR_2 VALIDTO_2 ## 0 <NA> FR.BR <NA> Département Department 17900226 Unknown ## 1 <NA> FR.HR <NA> Département Department 17900226 Unknown ## 2 <NA> FR.DD <NA> Département Department 17900226 Unknown ## 3 <NA> FR.GI <NA> Département Department 17900226 Unknown ## 4 <NA> FR.LD <NA> Département Department 17900226 Unknown ## 5 <NA> FR.LG <NA> Département Department 17900226 Unknown ## REMARKS_2 Shape_Leng Shape_Area ## 0 <NA> 4.538735 0.5840273 ## 1 <NA> 3.214178 0.4198797 ## 2 <NA> 5.012795 1.0389622 ## 3 <NA> 9.200047 1.1489822 ## 4 <NA> 5.531231 1.0372815 ## 5 <NA> 4.489830 0.6062017
# load or create data set.seed(100) myvar <- rnorm(1:96)
# manipulate data for the plot france.geodata <- data.frame(id=rownames(fr.prj@data), mapvariable=myvar) head(france.geodata)
## id mapvariable ## 1 0 1.12200636 ## 2 1 0.05912043 ## 3 2 -1.05873510 ## 4 3 -1.31513865 ## 5 4 0.32392954 ## 6 5 0.09152878
# fortify prepares the shape data for ggplot france.dataframe <- fortify(fr.prj) # convert to data frame for ggplot
## Regions defined for each Polygons
head(france.dataframe)
## long lat order hole piece id group ## 1 7.847912 49.04728 1 FALSE 1 0 0.1 ## 2 7.844539 49.04495 2 FALSE 1 0 0.1 ## 3 7.852439 49.04510 3 FALSE 1 0 0.1 ## 4 7.854333 49.04419 4 FALSE 1 0 0.1 ## 5 7.855955 49.04431 5 FALSE 1 0 0.1 ## 6 7.856299 49.03776 6 FALSE 1 0 0.1
#now combine the values by id values in both dataframes france.dat <- join(france.geodata, france.dataframe, by="id") head(france.dat)
## id mapvariable long lat order hole piece group ## 1 0 1.122006 7.847912 49.04728 1 FALSE 1 0.1 ## 2 0 1.122006 7.844539 49.04495 2 FALSE 1 0.1 ## 3 0 1.122006 7.852439 49.04510 3 FALSE 1 0.1 ## 4 0 1.122006 7.854333 49.04419 4 FALSE 1 0.1 ## 5 0 1.122006 7.855955 49.04431 5 FALSE 1 0.1 ## 6 0 1.122006 7.856299 49.03776 6 FALSE 1 0.1
# Plot 3
p <- ggplot(data=france.dat, aes(x=long, y=lat, group=group))
p <- p + geom_polygon(aes(fill=mapvariable)) +
geom_path(color="white",size=0.1) +
coord_equal() +
scale_fill_gradient(low = "#ffffcc", high = "#ff4444") +
labs(title="Our map",fill="My variable")
# plot the map
p
nclassint <- 5 #number of colors to be used in the palette cat <- classIntervals(myvar, nclassint,style = "jenks") #style refers to how the breaks are created colpal <- brewer.pal(nclassint,"RdBu") color <- findColours(cat,rev(colpal)) #sequential bins <- cat$brks lb <- length(bins)
plot(fr.prj, col=color,border=T)
legend("bottomleft",fill=rev(colpal),legend=paste(round(bins[-length(bins)],1),":",round(bins[-1],1)),cex=1, bg="white")
rm(list = ls(all=T)) #clear workspace library(dplyr) library(readr) library(stringr) library(tidyr) library(readxl) library(classInt) library(RColorBrewer) library(maptools) #to read shapefiles
2. Download the data files (note they are not ready for use but need some cleaning as there are more areas in the excel files than polygons in the shape file). I copy here the code as I have used it in my script but it’s available at RPubs thanks to David Robinson.
download.file("http://www2.census.gov/prod2/statcomp/usac/excel/LND01.xls", "LND01.xls")
download.file("http://www2.census.gov/prod2/statcomp/usac/excel/POP01.xls", "POP01.xls")
according to metadata, this is Land Area in 2010 and resident population in 2010:
us_county_area <- read_excel("LND01.xls")
transmute(CountyCode = as.character(as.integer(STCOU)), Area = LND110210D)
us_county_population <- read_excel("POP01.xls")
transmute(CountyCode = as.character(as.integer(STCOU)),Population = POP010210D)
3. Adjust data
election_url <- "https://raw.githubusercontent.com/Prooffreader/election_2016_data/master/data/presidential_general_election_2016_by_county.csv"
county_data <- read_csv(election_url)
group_by(CountyCode = as.character(fips))
mutate(TotalVotes = sum(votes))
ungroup()
mutate(name = str_replace(name, ".\\. ", ""))
filter(name %in% c("Trump", "Clinton", "Johnson", "Stein"))
transmute(County = str_replace(geo_name, " County", ""),
State = state,
CountyCode = as.character(fips),
Candidate = name,
Percent = vote_pct / 100,
TotalVotes)
spread(Candidate, Percent, fill = 0)
inner_join(us_county_population, by = "CountyCode")
inner_join(us_county_area, by = "CountyCode")
you can save the data into a csv file:
# write_csv(county_data, "county_election_2016.csv")
You can download the cleaned datafile here: data_election_2016_by_county
4. Upload data and shape files
setwd("/Users/...")
dt <- read.csv("new_county_election_2016.csv", header=T)
us <- readShapePoly("./USA_adm/USA_adm2.shp")
us0 <- readShapePoly("./USA_adm/USA_adm0.shp")
us.m <- us[-c(which(us$NAME_1=="Alaska")),] #get rid of Alaska
us.d <- us.m[-c(67:71),]
5. Prepare the color palette(s)
nclassint <- 5 #number of colors to be used in the palette cat.T <- classIntervals(dt$Trump[-c(67:71)], nclassint,style = "jenks") #style refers to how the breaks are created colpal.T <- brewer.pal(nclassint,"Reds") color.T <- findColours(cat.T,colpal.T) #sequential bins.T <- cat.T$brks lb.T <- length(bins.T)
5. Plot the maps with map basic…
# pdf("Where are the trump voters.pdf")
# plot(us.d, col=color.T, border=F)
# plot(us0,add=T, lwd=0.1)
# legend("bottomleft",fill=colpal.T,legend=paste(round(bins[-length(bins.T)],1),":",round(bins.T[-1],1)),cex=1, bg="white")
# dev.off()
… or ggplot2
library(ggplot2)
library(scales)
theme_set(theme_bw())
ggplot(county_data, aes(Population / Area, Trump)) +
geom_point() +
geom_point(data=county_data[which(county_data$State=="Texas"),], aes(x=Population/Area, y=Trump), colour="red")+
scale_x_log10() +
scale_y_continuous(labels = percent_format()) +
xlab("Population density (ppl / square mile)") +
ylab("% of votes going to Trump") +
geom_text(aes(label = County), vjust = 1, hjust = 1, check_overlap = TRUE) +
geom_smooth(method = "lm") +
ggtitle("Population density vs Trump voters by county (Texas Counties in red)")
This is the code to plot in red points according to State (in red) and to add red labels to those points. The check_overlap=T avoids overlapping labels.
# ggplot(county_data, aes(Population / Area, Trump)) +
# geom_point() +
# geom_point(data=county_data[which(county_data$State=="California"),], aes(x=Population/Area, y=Trump), colour="red")+
# scale_x_log10() +
# scale_y_continuous(labels = percent_format()) +
# xlab("Population density (ppl / square mile)") +
# ylab("% of votes going to Trump") +
# geom_text(data=county_data[which(county_data$State=="California"),], aes(label = ifelse(Trump&gt;.5, as.character(dt$County), "" )), color= "red",size=5,vjust = 1, hjust = 1, check_overlap = TRUE) +
# geom_smooth(method = "lm") +
# ggtitle("Population density vs Trump voters by county (California in red)")
Using data on single men (21-35 years old) and women (16-30 years old) you can map where the unbalances in the marriage market are caused by “excess” of male or female population
Using data on single men (21-35 years old) and women (16-30 years old) you can map where the unbalances in the marriage market are caused by “excess” of male or female population.
Not every shape file is as nice as those provided in libraries. Sometimes we have to deal with historical maps, which have been hand-drawn, re-touched and what not. To work with geo-referenced data it is essential to have a variable in both shape file and dataframe with unique coding that has exactly the same number of areas and the same ordering in both files.
A quick way to check if shapefile and dataframe have the same number of areas:
nrow(df) == length(shape.file$Code)
In the shapefile, one can also select a couple of areas big enough so that they can easily be located, and plot them as “control” areas.
For instance, I want to select the area with code “15078” in the shapefile:
>which(shape.file$Code=="15078",arr.Ind=T)
[1] 271
which is the area in the 271-th position (same way shape.file$Code[271] gives the code of area 271).
plot(shape.file)
plot(shape.file[c(271,898),],col="red",border="red",add=T)
this is an easy way to locate your “control” area(s).
Ideally, you should have some variable that is identical to the one in the shapefile, a codification of some sort, providing a unique Code, the name of the area or some factors that allow you to locate the area in space.
An easy way to check if both shape file and data frame have the same ordering of geographical areas is to test it:
>code.sh <- cbind(c(1:length(shape.file$Code)),as.vector(shape.file$Code))
>code.df <- cbind(c(1:nrow(df)),df$Code)
>code.df==code.sh
[,1] [,2]
[1,] TRUE TRUE
[2,] TRUE TRUE
[3,] TRUE TRUE
What if it’s not?
First option: the inelegant solution
Manually change the order of the areas in a csv file according to the exact order they have in the shape file. It’s easy as you can create an ordinal index for the shapefile codes, paste it in excel, and assign it with a vlookup function.
Second option: the smart R match
In R there is a function called match that returns a vector of the positions of first matches of the first argument in its second:
>my.match <- match(df$Code, shape.file$Code)
NB: to use match the two variables providing the code for the areas have to have the very same unique and identical codes, or else funny stuff happens. To check that everything is in its right place, you can plot the two “control” spatial polygons we chose in the beginning, using their position in the dataframe rather than in the shapefile:
>plot(shape.file)
>plot(shape.file[c(which(df$Code=="305"),which(df$Code=="15078")),],col="orange",add=T)
The map of GOT world with rivers, roads, lakes, the Wall, and main cities:
Neighborhood relations according to Sphere of Influence pretty much coincide with roads and rivers (package spdep):
Paste some images to locate the (surviving) Stark family members, using rasterImage from the png library:
