# Project raster when some of extent out of range

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 <- 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)
``````

``````# 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 rgdal::rawTransform(projfrom, projto, nrow(xy), xy[, 1], xy[,  :
36 projected point(s) not finite

rasterTemplateLongLat <- projectRaster(rasterTemplate, to = e)

``````

Any suggestions?

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 <- 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)
``````

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.

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