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(po...@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@stat.math.ethz.ch
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
