问题描述
我想使用R中sp软件包中的over()函数.
I want to use the over() function from the sp package in 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"
让我们检查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()函数
So, when I run the over() function now
pointsInPolygon <- over(allPoints, finalPolygon))
我得到了错误:
我想我知道问题出在哪里，但我不知道如何解决.
I think I know what the problem is here, but I do not know how to solve it.
如果仔细观察，allPoints还有几个词-即+init=epsg:3857.我在此处阅读，sp package只是比较CRS插槽中的 strings 是否相同.好吧，它们位于您所看到的相同的CRS中(Spatial引用的creatin完全相同)，但是由于我创建它们的过程，字符串略有不同.
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.
当我使用
#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上，这让我反感:
over()函数可以正常工作，但是我得到的结果没有任何意义.
The over() function works then, but what I get back does not make sense.
如何解决这个问题?
推荐答案
您应该已经创建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")，与
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
这是您用于cornersEPSG3857的内容.尽管这两个字符串可能表示相同的CRS(即，具有相同的语义)，但是它们在语法上有所不同，并且sp或基础的PROJ.4库(GDAL的一部分，通过包rgdal进行接口)都没有用于比较的函数.两个语法上不同的proj4strings的语义.
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.
此问题的解决方案是一次定义新的CRS，例如由
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.
(最后的警告是BTW，以确保用户仅覆盖CRS时不会认为他们在重新投影；您 did 覆盖了它，然后 did 也可以解决您的问题，但是用一种不太干净的方式
(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)
这篇关于在sp中添加CRS似乎不一致的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持！