2

我想over()使用sp.R

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

让我们检查一下CRS。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"

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

pointsInPolygon <- over(allPoints, finalPolygon))

我得到错误：

相同CRS(x, y) 不为真

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

如果你仔细看，allPoints还有几个字——即+init=epsg:3857。我在这里读到，sp package简单地比较插槽中的字符串是否CRS相同。嗯，它们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我扔回去了：

proj4string<-( , value = )中的警告*tmp*：新的 CRS 已分配给具有现有 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 不重新投影。对于重投影，使用包 rgdal 中的函数 spTransform

然后该over()功能起作用，但是我得到的没有意义。

如何解决这个问题？！

4

1 回答 1

3

你应该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 的一部分，通过 package 接口rgdal）都没有比较两个语法上不同的 proj4string 的语义的功能。

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

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

并在整个脚本中使用它。

（最后的警告是顺便说一句，以确保用户在仅覆盖 CRS 时不会认为他们重新投影；您确实覆盖了它，它也确实解决了您的问题，但是以一种不太干净的方式）

于 2015-09-04T12:00:15.253 回答