Adding CRS in sp seems inconsistent

Question

I want to use the over() function from the sp package in R.

I assigne a CRS.

#say that polygon is EPSG3857 (Web Mercator PROJECTION)
proj4string(finalPolygon) <- CRS("+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext  +no_defs")

and all seems good.

str(finalPolygon)
> ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
> .. .. ..@ projargs: chr "+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +no_defs"

Lets check the CRS of allPoints.

str(allPoints)
>..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
> .. .. ..@ projargs: chr "+init=epsg:3857 +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +no_def"

So, when I run the over() function now

pointsInPolygon <- over(allPoints, finalPolygon))

I get the error:

identicalCRS(x, y) is not TRUE

I think I know what the problem is here, but I do not know how to solve it.

If you look closely, allPoints has a few words more - namely +init=epsg:3857. I read here that the sp package simply compares if the strings in the CRS slot are identical. Well, they are in the same CRS as you can see (the creatin of the Spatial reference is exactly the same), but the strings differ slightly due to the process how I created them.

When I use

#say that points is EPSG3857 (Web Mercator PROJECTION)
proj4string(spatialEPSG3857) <- CRS("+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext  +no_defs")

on the allPoints it throws me this back:

Warning in proj4string<-(*tmp*, value = ) : A new CRS was assigned to an object with an existing CRS: +init=epsg:3857 +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +no_defs without reprojecting. For reprojection, use function spTransform in package rgdal

The over() function works then, but what I get back does not make sense.

How to tackle this problem?!

Edzer Pebesma Edzer Pebesma · Accepted Answer · 2015-09-04T12:00:15

You should have created finalPolygon by

finalPolygon <- SpatialPolygons(list(myPolygon), proj4string = CRS(proj4string(cornersEPSG3857)))

as its documentation states that the CRS is set to NA by default. Instead, you set the CRS in the next statement to CRS("+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs") which is different from

> CRS("+init=epsg:3857")
CRS arguments:
 +init=epsg:3857 +ellps=WGS84 +proj=merc +a=6378137 +lat_ts=0.0
+lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +no_defs

which is what you used for cornersEPSG3857. Although these two strings may represent the same CRS (i.e. have the same semantics), they differ syntactically, and neither sp nor the underlying PROJ.4 library (part of GDAL, interfaced through package rgdal) have a function to compare semantics of two syntactically different proj4strings.

The solution to this problem is to define the new CRS once, e.g. by

crs3857 =  CRS("+init=epsg:3857")

and use that throughout your entire script.

(The warning at the end is BTW there to make sure users don't think they reproject when they merely overwrite a CRS; you did overwrite it, and it did solve your problem too, but in a not so clean way)

Adding CRS in sp seems inconsistent

1 Answers