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

finalPolygonCRS + over() function sp package in R

我这样创建一个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))

但是,我收到以下错误消息:

Warning: Unhandled error in observer: identicalCRS(x, y) is not 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")

我收到以下错误:

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

  1. 为什么会出现 CRS are not identical 错误?它们应该是相同的,因为 proj4stringspTransform 是相同的代码?!
  2. 为什么当我尝试将 proj4string 分配给 SpatialPoints 时会出现此错误,但在多边形上执行此操作时却不会?

我认为这个问题可以通过更加小心地在整个代码中指定 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