R中的finalPolygonCRS + over()函数sp包答案

【问题标题】：finalPolygonCRS + over() function sp package in RR中的finalPolygonCRS + over()函数sp包
【发布时间】：2015-11-29 09:06:15
【问题描述】：

我以这种方式创建SpatialPolygon

#make spatial points
#assign Original CRS WGS84 EPSG:4326 (DATUM --> on sphere)
cornersEPSG4326 <- SpatialPoints(coords=cbind(x,y), proj4string = CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"))

#transform to EPSG3857 (Web Mercator PROJECTION)
cornersEPSG3857 <- spTransform(cornersEPSG4326, CRS("+init=epsg:3857"))

#create Polygon
bbox <- Polygon(cornersEPSG3857)

#create PolygonsObject
myPolygon <- Polygons(list(bbox),1)

#create SpatialPolygonsObject
finalPolygon <- SpatialPolygons(list(myPolygon))

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

然后SpatialPoints 这样：

#make spatial points
#assign Original CRS WGS84 EPSG:4326 (DATUM --> on sphere)               
spatialEPSG4326 <- SpatialPoints(coords=cbind(x,y), proj4string = CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"))

#transform to EPSG3857 (Web Mercator PROJECTION)
spatialEPSG3857 <- spTransform(spatialEPSG4326, CRS("+init=epsg:3857"))

allPointsSpatial <- spatialEPSG3857

我想运行over()函数：

pointsInPolygon <- over(allPointsSpatial, finalPolygon)
print(length(pointsInPolygon))

但是，我收到以下错误消息：

警告：观察者中未处理的错误：相同CRS(x, y) 不是 TRUE

当我将此行添加到我的SpatialPoints时

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

我收到以下错误：

proj4string<-(*tmp*, value = ) 中的警告：新的 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 没有重新投影。对于重投影，使用函数 spTransform 在包rgdal

为什么会出现CRS are not identical 错误？它们应该是相同的，因为proj4string 和spTransform 是相同的代码？！
为什么当我尝试将 proj4string 分配给 SpatialPoints 时会出现此错误，但在 Polygon 上执行时却没有？

【问题讨论】：

您可以在问题中添加dput(x) 和dput(y) 吗？
你有什么理由不只是使用CRS("+init=epsg:4326")？如果没有样本数据，我无法测试这是否有效，但它似乎应该是因为它列在这里：spatialreference.org/ref/epsg/4326
@MichaelChirico 使用这个CRS("+init=epsg:4326") 而不是哪一行？！
您尝试使用墨卡托的任何地方，例如声明cornersEPSG4326
@MichaelChirico 我想我需要为此使用proj4string？！在这里查看spatialreference.org/ref/epsg/4326 那就是spatialreference.org/ref/epsg/4326/proj4 ？

标签： r geospatial spatial sp

【解决方案1】：

我认为这个问题可以通过更加小心地在整个代码中指定 CRS 来解决。

以下内容对我没有错误：

#test on the unit square
test.box<-as.matrix(rbind(c(0,0),
                          1:0,
                          c(1,1),
                          0:1,
                          c(0,0)))
#store CRS for clarity
EPSG4326<-CRS("+init=epsg:4326")
EPSG3857<-CRS("+init=epsg:3857")

#re-run code:
cornersEPSG4326 <- 
  SpatialPoints(coords=test.box, proj4string = EPSG4326)
cornersEPSG3857 <- spTransform(cornersEPSG4326, EPSG3857)
bbox <- Polygon(cornersEPSG3857)
myPolygon <- Polygons(list(bbox),1)
finalPolygon <- SpatialPolygons(list(myPolygon))
proj4string(finalPolygon) <- EPSG3857

spatialEPSG4326 <- 
  SpatialPoints(coords=test.box, proj4string = EPSG4326)
spatialEPSG3857 <- spTransform(spatialEPSG4326, EPSG3857)
allPointsSpatial <- spatialEPSG3857

#output: no errors
pointsInPolygon <- over(allPointsSpatial, finalPolygon)
> cat(length(pointsInPolygon))
5

【讨论】：