在 sp 中添加 CRS 似乎不一致

Question

我想使用 R 中 sp 包中的 over() 函数。

我分配了一个 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")

一切似乎都很好。

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"

让我们检查 allPoints 的 CRS。

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"

所以，当我运行 over() 函数现在

pointsInPolygon <- over(allPoints, finalPolygon))

我收到错误：

identicalCRS(x, y) is not TRUE

我想我知道这里的问题是什么，但我不知道如何解决。

仔细一看，allPoints多了几个字，就是+init=epsg:3857。我读到 here sp package 只是比较 CRS 插槽中的 strings 是否相同。好吧，正如您所见，它们在相同的 CRS 中（空间参考的创建完全相同），但由于我创建它们的过程，字符串略有不同。

当我使用

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

在 allPoints 上它把我扔回去了：

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

over() 函数可以正常工作，但我得到的结果没有任何意义。

如何解决这个问题？！

Answer 1

您应该通过

创建finalPolygon

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

因为它的文档指出 CRS 默认设置为 NA。相反，您将下一条语句中的 CRS 设置为 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")，这与

不同

> 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

这是您用于 cornersEPSG3857 的内容。尽管这两个字符串可能表示相同的 CRS（即具有相同的语义），但它们在句法上有所不同，sp 和底层 PROJ.4 库（GDAL 的一部分，通过包 rgdal 连接） ) 具有比较两个句法不同的 proj4strings 语义的功能。

解决这个问题的方法是定义一次新的CRS，例如通过

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

并在整个脚本中使用它。

（顺便说一句，最后的警告是为了确保用户不会认为他们只是重写了 CRS；你确实重写了它，它确实也解决了你的问题，但不是很干净)

在 sp 中添加 CRS 似乎不一致

Adding CRS in sp seems inconsistent

r

spatial

geospatial

sp