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)")
