0

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 – mdsumner Aug 19 '16 at 0:01
  • 1
    Please don't add "thank you" as an answer. Instead, accept the answer that you found most helpful. - From Review – Andre Silva Aug 19 '16 at 20:14
1

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.

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.