R Raster – Resolving Warnings for Projected Points Not Finite

coordinate systemrrasterrgdalsp

1. Question

I have encountered a warning using the projectRaster() function in the raster package in R. A full reproducible example is pasted below.

   Warning message:
   In rgdal::rawTransform(projto_int, projfrom, nrow(xy), xy[, 1],  :
   33940 projected point(s) not finite

My question is: Is this warning a problem I need to fix if I am working on terrestrial data? In other words is data "lost". This would be a major problem for me if it is. If that is the case, do you know if there a way I can fix it?

I have searched for a solution to this problem online and found some mention of it here, here & here but I think none offer a suitable answer to this problem.

2. Load the raster library

  library(raster)

3. First create a (long lat) map of the world, with a constant in each grid cell

I am putting a constant in each grid cell to see if I can diagnose the problem of where the warning is likely to affect the data, if at all.

 rastertest.longlat<-raster(ncol=360, nrow=180)
 values(rastertest.longlat)<-c(rep(1,n=180*360))

4. If you reproject the (long lat) map onto an equal area grid, it produces the warning message

   rastertest.eck4<-projectRaster(rastertest.longlat, res=c(100000,100000), crs="+proj=eck4",method="ngb", over=T)


Warning message:
 In rgdal::rawTransform(projto_int, projfrom, nrow(xy), xy[, 1],  :
33940 projected point(s) not finite

I think this message is basically is saying that the the re-projection failed for some of the grid cells.

5. But if you plot out the two maps, it doesn't look like this warning poses any problems for the data

That is you don't see any white gaps in the data plotted. My guess is that the lost cells are non terrestrial cells, at the top and the edge of the world extent. Any ideas?

par(mfrow=c(2,1))

plot(rastertest.longlat, col="blue")
data(wrld_simpl)
wrld <- spTransform(wrld_simpl, CRS('+proj=longlat'))
plot(wrld, add=TRUE)

plot(rastertest.eck4, col="blue") 
wrld <- spTransform(wrld_simpl, CRS('+proj=eck4'))
plot(wrld, add=TRUE)

enter image description here

Best Answer

Short answer: it's fine, and you can thank raster for proceeding to completion instead of failure, and for letting you know that some data were lost.

Long answer:

it will depend on the projection, and in this case it's probably just on the "the edges". what the edge is and how it manifests for a given instance of a given projection family is the "depends" part.

You can see that it's not the centre points of the cells that are lost: 

tpoints <- rgdal::project(coordinates(rastertest.longlat), "+proj=eck4")
sum(is.na(tpoints))
#[1] 0

But it probably is the corners, and possibly the edges of very cell. This perhaps shows that raster projects based on the extent of cells, not just their centre points. 

 rgdal::project(as.matrix(expand.grid(x = c(-180, 0, 180), y = c(-90,0, 90))), "+proj=eck4")

I admit I expected that to be where the missing values come from, so maybe projectRaster is extending out a little further north and south? Set values there for latitude outside the -90/90 range and you start getting the warning. I'll follow up if I get a chance to explore more.

Finally, you should probably use an explicit ellipsoid or datum parameter, i.e. "+proj=eck4 +ellps=WGS84".

Related Question