1

我用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<-( , 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

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

1 回答 1

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
于 2015-09-03T21:19:19.700 回答