finalPolygonCRS + пакет функций over () в 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))

Но я получаю следующее сообщение об ошибке:

Предупреждение: необработанная ошибка в наблюдателе: идентичный CRS (x, y) НЕ ИСТИНА

Когда я добавляю эту строку в свой 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")

Я получаю следующую ошибку:

Предупреждение в _9 _ (_ 10_, 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

  1. Почему возникает ошибка CRS are not identical? Они должны быть одинаковыми, поскольку proj4string и spTransform - это один и тот же код ?!
  2. Почему эта ошибка возникает, когда я пытаюсь назначить proj4string SpatialPoints, но не когда я делаю это на полигоне?

person Stophface    schedule 03.09.2015    source источник
comment
Вы можете добавить dput(x) и dput(y) к своему вопросу?   -  person MichaelChirico    schedule 03.09.2015
comment
Есть ли причина, по которой вы просто не используете CRS("+init=epsg:4326")? Я не могу проверить, работает ли это без образцов данных, но, похоже, это должно быть, потому что он указан здесь: пространственная ссылка. org / ref / epsg / 4326   -  person MichaelChirico    schedule 03.09.2015
comment
@MichaelChirico использует CRS("+init=epsg:4326") вместо какой строки ?!   -  person Stophface    schedule 03.09.2015
comment
везде, где вы пытаетесь использовать Mercator, например объявление cornersEPSG4326   -  person MichaelChirico    schedule 03.09.2015
comment
@MichaelChirico Я думал, мне нужно использовать proj4string для этого ?! Проверяя здесь Spacereference.org/ref/epsg/4326, это будет то, что Spacereference.org/ref/epsg/4326/proj4?   -  person Stophface    schedule 03.09.2015
comment
Хм, не уверен, что слежу за вами, но, например, не работает cornersEPSG4326 <- SpatialPoints(coords=cbind(x,y), proj4string = CRS("+init=epsg:4326"))?   -  person MichaelChirico    schedule 03.09.2015
comment
и аналогично для использования CRS("+init=epsg:3587")   -  person MichaelChirico    schedule 03.09.2015
comment
@MichaelChirico ага. запуталась с цифрами. Хотя в коде все правильно. Хорошо, если вы посмотрите здесь cran.r-project.org/web/packages /sp/sp.pdf в нем говорится, что при создании SpatialPoints вам необходимо определить proj4string=CRS(as.character(NA)). Если вы посмотрите proj4string здесь spacereference.org/ref/epsg/4326, я написал то же самое! Или я совершенно не прав?   -  person Stophface    schedule 03.09.2015
comment
Позвольте нам продолжить это обсуждение в чате.   -  person MichaelChirico    schedule 03.09.2015
comment
возможный дубликат Добавление CRS в sp кажется несовместимым   -  person Edzer Pebesma    schedule 04.09.2015


Ответы (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
person MichaelChirico    schedule 03.09.2015