0

I have a shapefile obtained from this link.

library(rgdal)
library(rgeos)

dir <- tempdir()
unzip('~/Downloads/4106902.zip', exdir = paste0(dir, '/4106902'))

shape <- readOGR(paste0(dir, '/4106902'), '4106902')
shape <- gBuffer(shape, byid = T, width = 0) #correcting orphaned hole problem

When I check the CRS associated with the file I get:

> proj4string(shape)
[1] "+proj=utm +zone=22 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"

This is equivalent to the EPSG:32722, according to spatialreference.org.
I need to locate a point obtained from Google Geocode API in this shapefile, so I'll read the latlong values and reproject them to the EPSG:32722:

google <- CRS("+init=epsg:4326")
point <- SpatialPoints(cbind(-49.25848, -25.45543),   
                       proj4string = google)
point.rep <- spTransform(point, CRS('+init=epsg:32722')) 

> summary(point.rep)
Object of class SpatialPoints
Coordinates:
          min       max
lng  675099.7  675099.7
lat 7183477.2 7183477.2
Is projected: TRUE 
proj4string :
[+init=epsg:32722 +proj=utm +zone=22 +south +datum=WGS84 +units=m +no_defs
+ellps=WGS84 +towgs84=0,0,0]
Number of points: 1

However, when I try to match the locations, I receive an error message stating that the CRS are not identical.

> over(point.rep, shape)
Error in .local(x, y, returnList, fn, ...) : 
  identicalCRS(x, y) is not TRUE

If I try to reproject the coordinates from google using the CRS extracted directly from the shapefile, instead of using the equivalent EPSG code, I get wrong lat value from the reprojection:

point.rep2 <- spTransform(point, proj4string(shape))

> summary(point.rep2)
Object of class SpatialPoints
Coordinates:
           min        max
lng   675099.7   675099.7
lat -2816522.8 -2816522.8
Is projected: TRUE 
proj4string :
[+proj=utm +zone=22 +datum=WGS84 +units=m +no_defs +ellps=WGS84
+towgs84=0,0,0]
Number of points: 1

> over(point.rep2, shape)
   ID_ NAME1_ NAME2_ PARTS_ POINTS_ LENGTH_ AREA_
1 <NA>   <NA>   <NA>     NA      NA      NA    NA

My current workaround is to force the CRS of my points after the reprojection, allowing me to find the correct location of the point in my shapefile:

> proj4string(point.rep) <- proj4string(shape) 
Warning message:
In ReplProj4string(obj, CRS(value)) :
  A new CRS was assigned to an object with an existing CRS:
+init=epsg:32722 +proj=utm +zone=22 +south +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
without reprojecting.
For reprojection, use function spTransform

> over(point.rep, shape)
              ID_ NAME1_ NAME2_ PARTS_ POINTS_  LENGTH_    AREA_
1 410690205010318   <NA>   <NA>      1      27 1.661884 0.115704

Am I missing something? Why do I fail do reproject my point using the CRS from the shapefile?

0

Indeed, you have objects with different projection strings.

> identicalCRS(shape, point.rep)
[1] FALSE

Inspecting them makes it clearer:

> shape@proj4string
CRS arguments:
 +proj=utm +zone=22 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
> point.rep@proj4string
CRS arguments:
 +init=epsg:32722 +proj=utm +zone=22 +south +datum=WGS84 +units=m +no_defs +ellps=WGS84
+towgs84=0,0,0 

you will note that after the reprojection they do not have the same identical CRS string. Reprojecting the two of them to CRS("+proj=longlat +ellps=WGS84 +datum=WGS84") worked out for me:

> crs.geo <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84")
> shape <- spTransform(shape, crs.geo)
> point.rep <- spTransform(point, crs.geo)
> identicalCRS(shape, point.rep)
[1] TRUE

and so your over() command should work out. I noted though that the coordinates you create are out of your shape's bbox, so it finds no overlaying between the two objects.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.