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