CODEWITHSUNDEEP
rgeospatialspatialr-sp
four-eyes
four-eyes

Reputation: 12394

finalPolygonCRS + over() function sp package in R

I create a SpatialPolygon this way

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

And SpatialPoints this way:

#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

I wanto run the over() function:

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

But, I recieve the following error message:

Warning: Unhandled error in observer: identicalCRS(x, y) is not TRUE

When I add this line to my 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")

I get the following error:

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. Why is the CRS are not identical error raised? They should be the same since the proj4string and the spTransform is the same code?!
  2. Why is this error raised when I try to assign the proj4string to the SpatialPoints, but not when I do it on the Polygon?

Upvotes: 1

Views: 1076

Answers (1)

MichaelChirico
MichaelChirico

Reputation: 34733

I think this problem can be solved by being more careful about how you specify your CRS throughout your code.

The following produces no errors for me:

#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

Upvotes: 1

Related Questions