[R-sig-Geo] How to merge SpatialPolygonsDataFrame?

Joan Maspons j.maspons at creaf.uab.cat
Tue Jun 21 19:34:50 CEST 2011


Hello,
I'm new in GIS from R and I have some problems. I have some points with many 
variables and I want to interpolate and apply a mask from country borders:

valuesES<- data.frame(LON=runif(n, min=-9, max=3), LAT=runif(n, min=36, 
max=43), val1=runif(n), val2=runif(n), valn=runif(n))

createMap1<- function(valuesES, fileName=NULL, country, ...){
  require(raster)
  r<-raster(nrows=300, ncols=300, xmn=min(values$LON)-1, 
xmx=max(values$LON)+1, ymn=min(values$LAT)-1, ymx=max(values$LAT)+1)
  r<-rasterize(as.matrix(values[,1:2]), r)
  r[cellFromXY(r, values[,1:2])]<- values[[3]]
  tps<- Tps(values[,1:2], values[,3])
  p <- raster(r)
  p <- interpolate(p, tps)
  c<- getData("GADM", country=country, level=0)
  rc<- rasterize(c, r)
  pc<- mask(p, rc)
#   png(fileName)
  plot(pc, ...)
#   dev.off()
}

createMap1(valuesES[,1:3], country="ESP")

This seems to work but I'm not sure about the interpolation method. The 
problem arise when I try to use more than one country as a mask:

createMap<- function(values, fileName, countries, ...){
  require(raster)
  require(fields)
  require(maptools) # unionSpatialPolygons
  gpclibPermit() # polygon geometry computations in maptools

  r<-raster(nrows=300, ncols=300, xmn=min(values$LON)-1, 
xmx=max(values$LON)+1, ymn=min(values$LAT)-1, ymx=max(values$LAT)+1)
  rv<-rasterize(as.matrix(values[,1:2]), r)
  rv[cellFromXY(r, values[,1:2])]<- values[,3]
  tps<- Tps(values[,1:2], values[,3])
  p <- raster(r)
  p <- interpolate(p, tps)

  country<- list()
  pol<- list()
  pols<- list()
  data<- data.frame()
  for (j in 1:length(countries)){
    country[[j]]<- getData("GADM", country=countries[j], level=0)
    data<- rbind(data, country[[j]]@data)
    pol[[j]]<- Polygon(country[[j]]) # not to plotable
    pols[[j]]<- Polygons(list(pol[[j]]), countries[j]) # not to plotable
  }
  sPol<- SpatialPolygons(pols, 1:length(countries)) # not to plotable
  row.names(data)<- data$ISO
  sPolDF<- SpatialPolygonsDataFrame(sPol, data) # not to plotable
  spplot(sPolDF, "ISON") # not to plotable

#   uSPolDF<- unionSpatialPolygons(sPolDF, countries) 
#   spplot(uSPolDF, "ISON")

  reuro<- rasterize(sPolDF, r)

  pam<- mask(p, reuro)
  plot(pam, ...)

}

n<- 100
countries<- c("AUT", "BEL", "BGR", "CZE", "DEU", "DNK", "EST", "ESP", "FIN", 
"FRA", "GRC", "HUN", "IRL", "ITA", "LTU", "LUX", "LVA", "NLD", "POL", "PRT", 
"ROU", "SWE", "SVN", "SVK")

values<- data.frame(LON=runif(n, min=-9, max=31), LAT=runif(n, min=36, 
max=70), val1=runif(n), val2=runif(n), valn=runif(n))
valuesES<- data.frame(LON=runif(n, min=-9, max=3), LAT=runif(n, min=36, 
max=43), val1=runif(n), val2=runif(n), valn=runif(n))

for (i in 3:(ncol(values)-2)){
  createMap(values[, c(1:2, i)], fileName=names(values)[i], 
countries=countries, main=paste(substr(names(values)[i], 1, 
nchar(names(values)[i])-4)))
}

I've tried to create a SpatialPolygonsDataFrame containing all the countries 
but it don't display anything. Where is the mistake? Could you give me a hand?

> sessionInfo()
R version 2.11.1 (2010-05-31) 
i686-pc-linux-gnu 

locale:
 [1] LC_CTYPE=ca_ES.UTF-8          LC_NUMERIC=C                  
LC_TIME=ca_ES.UTF-8          
 [4] LC_COLLATE=ca_ES.UTF-8        LC_MONETARY=ca_ES.UTF-8       
LC_MESSAGES=ca_ES.UTF-8      
 [7] LC_PAPER=ca_ES.UTF-8          LC_NAME=ca_ES.UTF-8           
LC_ADDRESS=ca_ES.UTF-8       
[10] LC_TELEPHONE=ca_ES.UTF-8      LC_MEASUREMENT=ca_ES.UTF-8    
LC_IDENTIFICATION=ca_ES.UTF-8

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] gpclib_1.5-1    maptools_0.7-34 lattice_0.18-8  foreign_0.8-40  fields_6.3      
spam_0.23-0     raster_1.8-22  
[8] sp_0.9-81       rkward_0.5.3   

loaded via a namespace (and not attached):
[1] grid_2.11.1  tools_2.11.1


Yours,

Joan Maspons 
CREAF (Centre de Recerca Ecològica i Aplicacions Forestals)
Universitat Autònoma de Barcelona, 08193 Bellaterra (Barcelona), Catalonia
Tel +34 93 581 2915            j.maspons at creaf.uab.cat
http://www.creaf.uab.cat



More information about the R-sig-Geo mailing list