[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