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

[Robert J. Hijmans]

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
>