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 DIVA-GIS or GADM, CDC/shape. For colors we will consider options from RColorBrewer.

The descriptions of maps are stored in shape/ESRI files. A map description consists of at least three files: dbf (data), shp (shapes), shx (index). For Slovenian communes download Slovenija; additional data can be obtained from SURS and 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: buletine2010/bsl_5.pdf / page 128. Romania files.

Shape file: 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 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 file.