Re projecting rasters (projectRaster)
Dear Bart, You have data NE of the Caribbean, and project these to the whole world. The below illustrates what goes wrong: library(raster) library(maptools) data(wrld_simpl) # 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") a[] <- 1:ncell(a) # project to longlat ageo <- projectRaster(a, raster()) plot(ageo) plot(wrld_simpl, add=T) # Oops a copy of the values in Australia. This happens because the coordinates of those locations transform to the same coordinates in the Azimuthal Equidistant projection (going "through the world"; projectRaster uses the reverse of the projection you ask it to do, i.e. it transforms coordinates from the output raster to those of the input raster). I had not considered that someone would do this type of transfer (from a small part to the entire world) but there is nothing wrong with this and I will fix projectRaster for these cases (which should also make the function much faster too for these cases). Robert
On Fri, Nov 26, 2010 at 7:49 AM, bart <bart at njn.nl> wrote:
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)
_______________________________________________ R-sig-Geo mailing list R-sig-Geo at stat.math.ethz.ch https://stat.ethz.ch/mailman/listinfo/r-sig-geo