[R-sig-Geo] Re projecting rasters (projectRaster)

bart bart at njn.nl
Fri Nov 26 16:49:11 CET 2010


Hi All

I have several rasters in an equal distance projection spread around the
world. My goal is to get them all in one mollweide projection. But there
happens something weird, after using projectRaster i get my raster
replicated in two locations on the world. See the example below. I think
it is not the reprojection because if i do that manually without the
raster i do only get one location/ continues raster.

Bart


require(rgdal)
require(raster)
require(maps)

# create example raster
h<-1000000
a<-raster(matrix(1:100, nrow=10),xmn=-h, xmx=h, ymn=-h,ymx=h,
crs="+proj=aeqd +lon_0=-52.577189655172 +lat_0=21.5126379310345")

# create world wide molllweide raster
grd <- GridTopology(c(-135,-67.5), c(89.99999,44.99999), c(4,4))
polys <- as.SpatialPolygons.GridTopology(grd)
grid <- SpatialPolygonsDataFrame(polys,
data=data.frame(dummy=rep("dummy", length(polys at plotOrder)),
row.names=sapply(slot(polys, "polygons"), function(i) slot(i, "ID"))))
proj4string(grid) <-  CRS("+proj=longlat +datum=WGS84")
#gridProj <- spTransform(grid, CRS("+proj=moll +datum=WGS84"))
gridProj <- spTransform(grid, CRS("+proj=moll +datum=WGS84"))
moll.grid <- spsample(gridProj, type="regular",
offset=bbox(gridProj)[,1]+12500, cellsize=50000)
gridded(moll.grid) <- TRUE
world.raster <- raster(moll.grid)

#show raster before reprojecting
plot(a)
#raster after reprojecting  now in two locations
plot(projectRaster(a, world.raster, method="ngb"))

# add world to show locations (Australia and the atlantic)

mapw<-data.frame(map("world", plot =  FALSE)[c("x","y")])
names(mapw)<-c("location.long", "location.lat")
id<- cumsum(is.na(mapw$location.long))[!is.na(mapw$location.long)]
map<-cbind(mapw[!is.na(mapw$location.lat),], id)
coordinates(map) <- ~location.long+location.lat
proj4string(map)<-CRS("+proj=longlat +datum=WGS84")
map<-as.data.frame(spTransform(map, CRS("+proj=moll +datum=WGS84")))
p<-function(x, col="black")
{
         if(nrow(x)>1)
                 for(i in 1:(nrow(x)-1))
                         lines(x[i:(i+1),]$location.long,
x[i:(i+1),]$location.lat, col=col, lwd=.5)
}
map<-split(map, map$id)
         lapply(map, p)



# only do reprojection
tmp<-as.data.frame(as(a,"SpatialGridDataFrame"))
coordinates(tmp)<- ~s1+s2
proj4string(tmp)<-CRS("+proj=aeqd +lon_0=-52.577189655172
+lat_0=21.5126379310345")
tmp<-spTransform(tmp, CRS("+proj=moll +datum=WGS84"))
plot(tmp)



More information about the R-sig-Geo mailing list