Data source: EUROSTAT, personal elaboration of data. Population by age and sex at 1st of July 2020.
Regional age-sex composition in Germany by NUTS2 60 to 80+ years old
Category Archives: maptools
The Spatial side of Demography at EPC…
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
Georgia Mapping in R
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
Plot neighbouring countries
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 maps
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")
add labels for the countries
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)
let’s add everything (or almost everything) together
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)
check georgia@data
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
print labels on the map
labels for admin 2
coords2<- coordinates(georgia2[2:6,]) admin2 <- c(as.character(georgia2$NAME_2[1:5])) admin2
## [1] "Gagra" "Gali" "Gudauta" "Gulripshi" "Ochamchire"
Upload data from World Bank
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
Map the prevalence overweight w/h
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 and SAVE map:
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)
Plot maps with base mapping tools and ggmap in R
Plot maps with ‘base’ mapping tools in R
Understanding what kind of data you have (polygons or points?) and what you want to map is pivotal to start your mapping.
- First you need a shapefile of the area you want to plot, such as metropolitan France. There are various resources where to get them from: DIVA-GIS and EUROSTAT are those that I use the most. It’s always important to have a .prj file included, as your final map ‘should’ be projecte. I say “should” as sometimes it is just not possible, especially if you work with historical maps.
- Upload libraries
Load and prepare data
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
Use ggmap
# 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
Use plot basic
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")
