0

I have a raster that is in a North America Albers equal area projection. I would like to project this raster to unprojected long/lat, but I get an error. I understand that this is likely because some of the northern-most cells would transform to > 90 degrees, but I haven't been able to think of a way to work around this problem.

I've put together a reproducible example:

library(raster)
library(maptools)
data(wrld_simpl)

# generate a raster of north america
bb <- rbind(c(-166, 78), c(-175, 48), c(-108, 0), c(-55, 12), c(5, 78), c(-166, 78))
bb <- SpatialPoints(bb, proj4string = CRS('+proj=longlat +datum=WGS84'))
bb <- rgeos::gConvexHull(bb)

rasterTemplate <- raster(ext = extent(wrld_simpl), res = c(1,1))
rasterTemplate <- rasterize(wrld_simpl, rasterTemplate)
rasterTemplate[!is.na(rasterTemplate)] <- 1

rasterTemplate <- crop(rasterTemplate, bb)
rasterTemplate <- mask(rasterTemplate, bb)
rasterTemplate <- trim(rasterTemplate)

# here is our equal area projection
EAproj <- '+proj=aea +lat_1=20 +lat_2=60 +lat_0=40 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs'

# project to the EA projection
e <- projectExtent(rasterTemplate, crs = EAproj)
rasterTemplate <- projectRaster(rasterTemplate, to = e)
rasterTemplate <- trim(rasterTemplate)

plot(rasterTemplate, legend = FALSE)

enter image description here

# project raster to long/lat
e <- projectExtent(rasterTemplate, crs = '+proj=longlat +datum=WGS84')

Error in if (maxy == miny) { : missing value where TRUE/FALSE needed
In addition: Warning message:
In rgdal::rawTransform(projfrom, projto, nrow(xy), xy[, 1], xy[,  :
  36 projected point(s) not finite

rasterTemplateLongLat <- projectRaster(rasterTemplate, to = e)

Any suggestions?

1

Yes that is strange, perhaps a bug somewhere. A workaround is to do it manually:

Set the target extent and resolution you want with rr and then get the coordinates, project them to the source and extract.

Use method = "bilinear" for non-integer data.

  library(raster)
library(maptools)
data(wrld_simpl)

# generate a raster of north america
bb <- rbind(c(-166, 78), c(-175, 48), c(-108, 0), c(-55, 12), c(5, 78), c(-166, 78))
bb <- SpatialPoints(bb, proj4string = CRS('+proj=longlat +datum=WGS84'))
bb <- rgeos::gConvexHull(bb)

rasterTemplate <- raster(ext = extent(wrld_simpl), res = c(1,1))
rasterTemplate <- rasterize(wrld_simpl, rasterTemplate)
rasterTemplate[!is.na(rasterTemplate)] <- 1

rasterTemplate <- crop(rasterTemplate, bb)
rasterTemplate <- mask(rasterTemplate, bb)
rasterTemplate <- trim(rasterTemplate)

# here is our equal area projection
EAproj <- '+proj=aea +lat_1=20 +lat_2=60 +lat_0=40 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs'

# project to the EA projection
e <- projectExtent(rasterTemplate, crs = EAproj)
rasterTemplate <- projectRaster(rasterTemplate, to = e)
rasterTemplate <- trim(rasterTemplate)


rr <- raster::raster(raster::extent(-176, 4, -1, 79 ), res = 1, crs = '+proj=longlat +datum=WGS84')

## do it manually
rr[] <- raster::extract(rasterTemplate, rgdal::project(coordinates(rr), projection(rasterTemplate)), 
                        method = "simple")

plot(rr)

enter image description here

btw, fasterize is much faster than rasterize, but you need sf and end up with a mish-mash of format-specific code across packages. Or use the stars package.

2
  • This is great, thanks! I was able to additionally infer the extent by extracting coordinates of cells with non-NA values, and transforming those.
    – Pascal
    Apr 1 '20 at 17:52
  • And yes, fasterize is excellent, I was just focusing on simplicity here for the reproducible example.
    – Pascal
    Apr 1 '20 at 17:53

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.