# R: Error in checking point inside polygon

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

• 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. Commented Apr 9, 2018 at 13:47
• Or use `over` with `sp` class objects and save the pointless (no pun intended) `fortify` step? Commented Apr 9, 2018 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. Commented Apr 9, 2018 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. Commented Apr 9, 2018 at 17:12
• I think we can't do anything without your shapefile. Can you share it? Commented Apr 9, 2018 at 17:24

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>
``````
• 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 Commented Apr 10, 2018 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. Commented Apr 10, 2018 at 13:40