8

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

8

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".

  • Thanks! I just re-set the extent of the original raster to less/greater than +/- 70 degrees: rastertest.longlat <- raster(ncol=4320, nrow=2160, xmn=-180, xmx=180, ymn=-70, ymx=70) and then reprojected using: rastertest.eck4<-projectRaster(rastertest.longlat, res=c(100000,100000), crs="+proj=eck4,method="ngb", over=T). The warning completely disappears, so I think all the missing values must occur +/- ~70 degrees. – mike Aug 4 '16 at 1:45
  • @mike I am sorry but this is not an acceptable solution, losing data and information is not a solution. Qgis and Arcgis can do this operation without a problem, there must be some problem with projectRaster. I hope the developer can chime in. – Herman Toothrot Oct 23 '18 at 7:27
  • (-‸ლ) I welcome you to demonstrate exactly what QGIS and Arc do in this case, speculating is not really helpful. We have access to exactly what is going on in R here, there's no actual problem. – mdsumner Nov 12 '18 at 1:00
5

This is not a full answer to my original question, in terms of all the gritty details, but it does provide the interested reader somewhere to go.

In general data is lost and distorted during reprojection of rasters from longlat to equal area projections. This can be problematic for global analysis. If you can get your data in vector format, it is better to reproject using that format instead.

Here is one reference on the general problem. And here another on trying to quantify the loss. Hope that helps.

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.