# Finding lat/lon point in coord polygon

I wrote a function to determine which ith polygon with coordinates contains a given lon/lat point, using the polygons from the tz_world.shp file. I pass in lon = -77, lat = 42, which is in New York State. The ID from the 155th entry is 'America/New_York', which should be the correct answer. But, the use of point.inpoly returns a zero-row length object. Here is my code.

``````i = 155

timeZoneList = split( timeZonesShpFile, timeZonesShpFile\$TZID )

ps = lapply( timeZoneList, Polygon )

p1 = lapply(seq_along(ps), function(i) Polygons(list(ps[[i]]), ID = names( timeZoneList )[i] ) )

my_spatial_polys = SpatialPolygons( p1, proj4string = CRS("+proj=longlat +datum=WGS84") )

polyNames = sapply(slot(my_spatial_spdf, 'polygons'), function(i) slot(i, 'ID'))

pt = SpatialPointsDataFrame(cbind(lon,lat), data.frame(row=1), proj4string=CRS("+proj=longlat +datum=WGS84" ))

thisPoly = my_spatial_polys[ i ]

thisList = getSpPPolygonsIDSlots( thisPoly )

thisDf = data.frame( row = 1, row.names = thisList )

thisSpdf = SpatialPolygonsDataFrame( thisPoly, thisDf )

pIp = point.in.poly( pt, thisSpdf )

> pIp

[1] coordinates row         row.1
<0 rows> (or 0-length row.names)
``````

The screen-view of the coords shows typical bounding lon/lat points for NYS, which should surround my supplied lon/lat point.

• See ?over this has all been done, ?raster:: extract provides a bit easier use Commented Aug 19, 2016 at 0:01

I recommend to use the built-in tools, which admittedly are a bit confusing due to historical incremental development.

The most important functions are `sp::over` and `sp::extract`, and `extract` for these objects (polys and raw points) drives `over` anyway. Under the hood this uses `point.in.polygon` for single polygon to single point tests, but it's optimized to be as fast as possible by cunningly avoiding unnecessary tests. It's also designed to work with "holes" and to not return a match for a point that falls inside a polygon hole, it's complicated and hard to make all this happen. So I'm not discouraging you just saying it's already available.

If I understand what you want, the easiest way is with `raster::extract`.

``````## first we need a SpatialPolygonsDataFrame and some points
library(raster)
xy <- data.frame(lon = c(-77, -80), lat = c(35, 35.1))
dwg <- shapefile(system.file("shapes", "sids.shp", package = "maptools"))
``````

This returns a data.frame with the 'point.ID' and 'poly.ID' of dwg and xy, where they match. These values are the row numbers in the sense of `dwg[i, ]` and `xy[i, ]`. Note that I'm blithely ignoring the projection metadata, which you shouldn't.

``````extract(dwg, xy)
``````

It also returns all the attributes from the data in `dwg`, so you don't have to index it back yourself, but you can. To really strip this down, cast the dwg to 'SpatialPolygons' to drop the data attributes, just depends on how it's best in your workflow.

e.g.

`````` extract(geometry(dwg), xy)
#  point.ID poly.ID
# 1        1      91
# 2        2      89
``````

Please provide a reproducible example for questions like this, it's much easier for an answer to be provided if you do that. I had to do all that first. It's helpful in many ways, 1) the answerer doesn't have to think about a relevant example 2) the answerer can try the solution that comes to mind and be sure that it's actually going to work. 3) the example is "ready-made" and easy for others to try.