====== Maps in R ======
====== Maps ======
The basic library in R for dealing with spatial data is
**''sp''**. Some useful function can be found also in libraries
**''maps''**, **''mapproj''**,
**''maptools''** and **''rgdal''**.
For producing map displays we will use the library **''maptools''**.
Some maps can be found in **''mapdata''**, many more on WWW -- for example
[[https://www.diva-gis.org/gdata|DIVA-GIS]]
or [[http://gadm.org/|GADM]],
[[https://www.cdc.gov/epiinfo/support/downloads/shapefiles.html|CDC/shape]].
For colors we will consider options from **''RColorBrewer''**.
The descriptions of maps are stored in
**''shape''**/[[http://www.esri.com/news/arcuser/0401/topo.html|ESRI]] files.
A map description consists of at least three files:
**''dbf''** (data), **''shp''** (shapes), **''shx''** (index).
For Slovenian communes download {{notes:zip:sicommunes.zip|Slovenija}};
additional data can be obtained from
[[http://www.stat.si/StatWeb/en|SURS]] and [[https://www.e-prostor.gov.si/access-to-geodetic-data/free-geodetic-data/|Free geodetic data]].
===== Example: Slovenia: proportion of foreigners =====
> setwd("C:/Users/.....")
> library(maptools)
> pop <- read.table("population.csv", sep = ";", as.is = TRUE,
+ col.names = c("commune", "slovenians", "others"))
> pop$total <- pop$slovenians+pop$others
> pop$proportion <- pop$others/pop$total
> pop[order(pop$proportion), ]
commune slovenians others total proportion
27 Dobje 954 0 954 0.000000000
191 Velika Polana 1478 0 1478 0.000000000
161 Sveti Andraž v Slov. goricah 1208 2 1210 0.001652893
164 Sveti Tomaž 2106 4 2110 0.001895735
15 Cankova 1903 4 1907 0.002097535
........
68 Koper/Capodistria 48171 3744 51915 0.072117885
26 Divača 3578 279 3857 0.072336012
117 Osilnica 369 33 402 0.082089552
150 Sežana 11720 1108 12828 0.086373558
87 Loška dolina 3590 340 3930 0.086513995
56 Izola/Isola 14482 1553 16035 0.096850639
> SIcom <- readShapeSpatial("OB/OB.shp")
> comOrder <- rank(SIcom$OB_UIME)
> SIcom$proportion <- pop$proportion[comOrder]
> spplot(SIcom,"proportion")
===== Example: Romania: unemployment May 2010 =====
Data: [[http://www.insse.ro/cms/files/arhiva_buletine2010/bsl_5.pdf|buletine2010/bsl_5.pdf]] / page 128. {{notes:zip:romania.zip|Romania files}}.
Shape file: [[https://gadm.org/download_country_v3.html|GADM]].
County Ilfov is not in the shape file.
> setwd("C:/Users/.....")
> library(maptools)
> gpclibPermit()
> ro <- read.csv2("unemployedRO.txt",skip=3,colClasses = "character")
> ROcou <- readShapeSpatial("epi/ro.shp")
> names(ROcou)
> cbind(ro[[1]],levels(ROcou$ADMIN_NAME))
> un <- sapply(ro[,-1],as.numeric)
> rownames(un) <- ro[[1]]
> rate <- un[1:41,"Rate.T"]
> brks <- quantile(rate, seq(0,1,1/7))
> library(RColorBrewer)
> pal <- brewer.pal(8,'Greens')
> plot(ROcou,col=pal[findInterval(rate,brks,all.inside=TRUE)])
> title("Unemployment May 2010")
> text(coordinates(ROcou),labels=as.character(ROcou$ADMIN_NAME),cex=0.4)
> legend("topright",
+ legend=c("2.5-6.2","6.2-7.8","7.8-8.2","8.2-8.7","8.7-9.8","9.8-11.2","11.2-14.1"),
+ fill=pal,cex=0.75)
===== Example: US counties =====
Shape data from [[http://gadm.org/country|GADM]].
> USout <- readShapeSpatial("USA/USA_adm0.shp")
> USsta <- readShapeSpatial("USA/USA_adm1.shp")
> UScou <- readShapeSpatial("USA/USA_adm2.shp")
> pdf("UScounties.pdf",width=11.7,height=8.3,paper="a4r")
> plot(UScou,xlim=c(-124,-67),ylim=c(23,48),col="wheat",bg="skyblue",border="red",lwd=0.05)
> plot(USsta,xlim=c(-124,-67),ylim=c(23,48),lwd=0.2,border="blue",add=TRUE)
> plot(USout,xlim=c(-124,-67),ylim=c(23,48),lwd=0.2,add=TRUE)
> title("Central US"); dev.off()
The obtained picture **''UScounties.pdf''** is very large (around 70 MB). It is available on a separate [[http://www.educa.fmf.uni-lj.si/didikta/gradivo/UScounties.pdf|file]].