0

I am trying to check if location point is inside the polygon and it works for certain location point. But for certain location point output says point is outside polygon even if it is inside. 

Rmap <- readShapePoly('MAP.shp')#rome shape file
mappoints <- fortify(Rmap)

data_raw<-data.frame(Longitude=12.477,Latitude=41.91020)
ind<-point.in.polygon( data_raw$Longitude, data_raw$Latitude,    mappoints$long, mappoints$lat)

ggplot(mappoints) + aes(long, lat) + geom_polygon(fill="grey80",aes   (group=group))  + scale_fill_manual(values = pretty_palette)+ geom_point(data=data_raw,aes(x=Longitude,y=Latitude))

Output: ind is always 0 which means point is outside polygon. Original plot is

plot of point and polygon from shapefile

Spasan
  • 13
  • 2
  • I am pretty sure that you're comparing two points with each other, instead of a point with a polygon. As seen in the help for "point.in.polygon", the last two arguments should be arrays: point.in.polygon(1:10,1:10,c(3,5,5,3),c(3,3,5,5)), as you want to check whether the point is within the extent of the polygon. – RJJoling Apr 09 '18 at 13:47
  • Or use over with sp class objects and save the pointless (no pun intended) fortify step? – Spacedman Apr 09 '18 at 14:57
  • @rjjoling you might be right. But when I plot the same data point and polygon coordinates, point overlapped with polygon. Doesn't that mean point is within the extent of polygon? Thanks for help. Would appreciate if you can further provide opinion on this one. – Spasan Apr 09 '18 at 17:05
  • @Spacedman I've tried with over too aith spatial polygons. And it returns NA for those point which were returned as 0 using point.polygon() even if they were inside polygon. Return value from over() and point.polygon() indicates point is outside polygon. I am not sure if it's my data or anything else in code tjat has issue. – Spasan Apr 09 '18 at 17:12
  • I think we can't do anything without your shapefile. Can you share it? – Spacedman Apr 09 '18 at 17:24
  • @Spacedman Here is the file link – Spasan Apr 10 '18 at 00:47

1 Answers1

1

You've not read the help for points.in.polygon correctly. The output from fortify describes all the polygons in the set, but the function only expects a single ring.

Using sp functions works perfectly. First DO NOT use readShapePoly because it doesn't read the coordinate system. Use rgdal or raster::shapefile which work properly:

> library(rgdal)
> Rmap = readOGR(".","roma")  # replace "." with the path to the folder

Now make points from your data frame. Everything is in lat-long coordinates so we can set the points coordinate system from the shapefile:

> data_raw<-data.frame(Longitude=12.477,Latitude=41.91020)
> coordinates(data_raw) = ~Longitude + Latitude
> proj4string(data_raw)=proj4string(Rmap)

Now we can over them in one of two ways, the first returns rows from the data or the polygons, the other returns row numbers for the polygons that you can use as indexes:

> over(data_raw, Rmap)
  ID_0 ISO NAME_0 ID_1 NAME_1 ID_2 NAME_2 ID_3 NAME_3 CCN_3 CCA_3  TYPE_3
1  112 ITA  Italy    8  Lazio   39   Roma 2507   Roma     0  <NA> Commune
  ENGTYPE_3 NL_NAME_3 VARNAME_3
1   Commune      <NA>      <NA>

> over(data_raw, as(Rmap,"SpatialPolygons"))
 1 
98 
> Rmap@data[98,]
   ID_0 ISO NAME_0 ID_1 NAME_1 ID_2 NAME_2 ID_3 NAME_3 CCN_3 CCA_3  TYPE_3
97  112 ITA  Italy    8  Lazio   39   Roma 2507   Roma     0  <NA> Commune
   ENGTYPE_3 NL_NAME_3 VARNAME_3
97   Commune      <NA>      <NA>
Spacedman
  • 63,755
  • 5
  • 81
  • 115
  • this helps and works better for most coordinates but for some coordinates like long:12.45561 and lat:41.90656 still returns NA although point is inside polygon – Spasan Apr 10 '18 at 13:21
  • There is a tiny hole in the middle of your shapefile. If you plot region 98 you'll see it a bit better: plot(Rmap[98,],col="red") - this is Rome so... its the Vatican. The Vatican is not included in your shapefile. You do not see this when you use ggplot because something is broken. I do not use ggplot for maps. – Spacedman Apr 10 '18 at 13:40