Given a Shapefile, understand in which polygon a point is

Question

Given a set of points (latitude and longitude) in the area of Milan,Italy I want to decide to which "neighborhood" each point belongs Here is what I have tried

library(shapefiles)
library(maptools)
library(maps)
library(sp)
library(rgdal)

mn.zip.map <- readShapePoly("NILZone.shp") #shape file that contains the shapes of the neighborhoods

#latlon inizialization and declaration

latlon$lat<-data$latitude
latlon$lon<-data$longitude

#at this point I have a data frame with N rows and 2 columns

coordinates(latlon) = ~lat+lon

So my first guess was:

over(latlon, mn.zip.map)

Of course this was a failure, when I tried

summary(mn.zip.map)
Object of class SpatialPolygonsDataFrame
Coordinates:
      min     max
x 1503202 1521750
y 5025952 5042528
Is projected: NA 
proj4string : [NA]
Data attributes: ....

     summary(latlon)
Object of class SpatialPoints
Coordinates:
          min       max
lat 45.391091 45.534113
lon  9.046796  9.274482
Is projected: NA 
proj4string : [NA]
Number of points: 155226

So what I tried was:

 map<-readOGR(dsn="C:...",layer="NILZone") #wHERE I FOUND THE RIGHT PARAMETERS TO CRS

 coordinates(latlon) = ~lat+lon
 proj4string(mn.zip.map) <- CRS("+proj=tmerc +lat_0=0 +lon_0=9 +k=0.9996       +x_0=1500000 +y_0=0 +ellps=intl +units=m +no_defs ")
 mn.zip.map <- spTransform(mn.zip.map, CRS("+proj=longlat"))
 proj4string(latlon) <- CRS(proj4string(mn.zip.map))
 head(over(latlon, mn.zip.map))

     FID_1 FID_1_1 ID_NIL  NIL AreaHA AreaMQ
 1    NA      NA     NA <NA>     NA     NA
 2    NA      NA     NA <NA>     NA     NA
 3    NA      NA     NA <NA>     NA     NA
 4    NA      NA     NA <NA>     NA     NA
 5    NA      NA     NA <NA>     NA     NA
 6    NA      NA     NA <NA>     NA     NA

Any guess of why this happens??

Thank you very much, sorry for the very long question!

EDIT: here you can download an archive.zip containing several files of the shapes and other information of the neighborhoods of Milan and a file latlon.txt which has the coordinates of 30 points that I'm interested to classify into a neighborhood

https://www.dropbox.com/s/l935nm3edhn5cc7/Archive.zip?dl=0

@G.Cocca I changed the parameters of CRS(...) now the regions should be overlapping!

mn.zip.map@bbox
    min       max
 x  9.040939  9.278398
 y 45.386074 45.535303

 latlon@bbox
      min       max
 lat 45.391091 45.534113
 lon  9.046796  9.274482

But still no result!

Your example is not reproducible and you will get better answers if people have data to play with. Can you make it reproducible, either by making your data available or by making your example with data included in the R packages you are using? — SlowLoris, May 12 '16 at 16:26
There is nothing wrong with what you did. Simply your polygons and your points do not overlay. See `mn.zip.map@bbox` and `latlon@bbox` — G. Cocca, May 12 '16 at 18:52
Thanks, now I'm wrapping my mind on this, since I have several files associated whit my .shp, might it be that the command CRS("+proj=utm +zone=32 +datum=WGS84") is wrong and the right arguments to pass are in a .prj file or somwhere else? How do I make sure I'm using the right ones? — mariob6, May 12 '16 at 19:31
I changed proj4string(mn.zip.map) <- CRS(..) and now as I posted above there seems to be an overlap to me! — mariob6, May 12 '16 at 20:54
Your axis orders are flipped from each other! The normal axis order is Cartesian (x y) or (lon lat), i.e. `coordinates(latlon) = ~lon+lat`. — Mike T, May 13 '16 at 00:28

score 1 · Answer 1 · answered May 15 '16 at 01:14

To turn a comment or two into an answer, here are some tips. To do any spatial analysis between two datasets, they must share the same coordinate reference system (CRS).

So either all in +proj=longlat or in a projected CRS (i.e. one of the datasets is tmerc which is a Transverse Mercator projection). As shown in the question, spTransform is a great way to translate geometry data from one CRS to another.

Lastly, if using +proj=longlat, the axis order is Cartesian (x y) or (long lat). So, for example, use coordinates(latlon) = ~lon+lat. This is a very common "gotcha".

Given a Shapefile, understand in which polygon a point is

1 Answers1