Thursday, 3 December 2015

Re-projecting raster in R: gives warning that projected point(s) not finite?



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.



  library(raster)


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



   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.



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



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


No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...