开发者

GADM-Maps cross-country comparison graphics

开发者 https://www.devze.com 2023-02-13 19:37 出处:网络
Maybe due to the fact I\'m relatively new to R, I have problems using the gadm-Mapfiles on http://www.gadm.org/.

Maybe due to the fact I'm relatively new to R, I have problems using the gadm-Mapfiles on http://www.gadm.org/.

I try to draw a map with several countries and compare them to each other (using different colors).

This is what I do

library('sp')
##
load(url('http://biogeo.ucdavis.edu/data/gadm2/R/ARG_adm0.RData')) 
# loads an Object "gadm" with shape of Argentinia
arg <- gadm # is there a more convenient way to do this in one line?
load(url('http://biogeo.ucdavis开发者_开发技巧.edu/data/gadm2/R/CHL_adm0.RData'))
# loads an Object "gadm" with shape of Chile
chl <-gadm
load(url('http://biogeo.ucdavis.edu/data/gadm2/R/BOL_adm0.RData'))
# loads an Object "gadm" with shape of Bolivia
bol <- gadm
##
spplot(c(arg, chl, bol))
# output: unable to find an inherited method for function "spplot", for signature "list"

Here are my problems:

  1. (This question is probably caused by my newbieness) Is there a more convenient way to load the shapefiles? I find it quite stupid to rename the gadm-Object all the time. Maybe there is even a way where R only downloads the data once and then has them stored in the workspace/somewhere locally?
  2. How can I convince R to plot all those countries on ONE map?

Thank you in advance!

[edit]

some nice functions With the help of Gavin Simpson, I was able to create some nice functions that reduce the whole map-merging to one line:

## you will need the sp-package
library('sp')

## load a file from GADM (you just have to specify the countries "special part" of the file name, like "ARG" for Argentina. Optionally you can specify which level you want to have
loadGADM <- function (fileName, level = 0, ...) {
    load(url(paste("http://biogeo.ucdavis.edu/data/gadm2/R/", fileName, "_adm", level, ".RData", sep     = "")))
    gadm
}

## the maps objects get a prefix (like "ARG_" for Argentina)
changeGADMPrefix <- function (GADM, prefix) {
    GADM <- spChFIDs(GADM, paste(prefix, row.names(GADM), sep = "_"))
    GADM
}

## load file and change prefix
loadChangePrefix <- function (fileName, level = 0, ...) {
    theFile <- loadGADM(fileName, level)
    theFile <- changeGADMPrefix(theFile, fileName)
    theFile
}

## this function creates a SpatialPolygonsDataFrame that contains all maps you specify in "fileNames".
## E.g.: 
## spdf <- getCountries(c("ARG","BOL","CHL"))
## plot(spdf) # should draw a map with Brasil, Argentina and Chile on it.
getCountries <- function (fileNames, level = 0, ...) {
    polygon <- sapply(fileNames, loadChangePrefix, level)
    polyMap <- do.call("rbind", polygon)
    polyMap
}

When you find this page, make sure you read this answer: https://stackoverflow.com/a/33264548/263589


The simplest solution for this is

library(raster)
bol <- getData('GADM', country='BOL', level=1)

plot(bol)

The R data were added to the GADM website to support this function. Also note that the file format has changed such that the other functions described on this page no longer work.

To combine different countries

per <- getData('GADM', country='PER', level=1)
bp <- bind(bol, per)
plot(bp)


For problem 1, this is R so you can roll your own load() function that does what you want, for example:

loadGADM <- function(file, ...) {
    load(file, ...)
    gadm
}

And use it as:

> ls()
character(0)
> loadGADM <- function(file, ...) {
+     load(file, ...)
+     gadm
+ }
> arg <- loadGADM(url('http://gadm.org/data/rda/ARG_adm0.RData'))
> ls()
[1] "arg"      "loadGADM"

This is a local solution when you know that the object loaded will be called gadm - you could improve the function to not need this, e.g.:

loadGADM <- function(file, ...) {
    f <- load(file, ...)
    get(f)
}

which works because load() returns the character strings of the names of the loaded objects.

For problem 2, you need to rbind() the three sp objects together, not concatenate them. However, this doesn't work for these objects and the Polygon IDs are non-unique:

> sa <- rbind(arg, chl, bol)
Error in validObject(res) : 
  invalid class "SpatialPolygons" object: non-unique Polygons ID slot values

I'm working on this and will update if I figure out the work around. The solution is to change the Polygons ID slot values using spChFIDs(). Here we append "arg_" etc to the rownames of the the objects such that these are no all unique:

arg <- spChFIDs(arg, paste("arg", row.names(arg), sep = "_"))
chl <- spChFIDs(chl, paste("chl", row.names(chl), sep = "_"))
bol <- spChFIDs(bol, paste("bol", row.names(bol), sep = "_"))
sa <- rbind(arg, chl, bol)

Then we can plot the combined sp object:

plot(sa) ## beware might take a long time to plot...


There some minor improvements that I can offer.

First, if you want to draw map graphics which shade countries according to some criteria, you have to load that value in the respective GADM file, e.g. from another data.frame. Here is how I did it:

## load a file from GADM (you just have to specify the countries "special part" of the file name, like "ARG" for Argentina. Optionally you can specify which level you want to have

loadGADM <- function (fileName, level = 0, ...) {
    load(paste("K:/Rdata/gadm/", fileName, "_adm", level, ".RData", sep = ""))
    gadm$coef <- land$coef[land[["abr"]]==fileName]
    gadm
}

If you draw maps of larger areas, your maps will suffer from distortion. Use the rgdal package (and its dependencies) to apply a map projection:

newproj <- CRS("+proj=wintri ellps=WGS84 +lon_0=15E")
spdf.wintri <- spTransform(spdf, newproj)

In order to plot the property coef, the package sp offers a nice method:

karte <- spplot(spdf.wintri,"coef",
        xlim=c(-13,46),ylim=c(33,72),
        col.regions = rainbow(100, start = 2/6, end = 4/6),
        main="Country Dummies")
0

精彩评论

暂无评论...
验证码 换一张
取 消