I have a lat/lon combination and want to check whether the point is inside a polygon (sp::Polygon
class)
Consider this example:
UKJ32 <- sp::Polygon(cbind(c(-1.477037449999955, -1.366895449999959, -1.365159449999965, -1.477037449999955),
c(50.923958250000027, 50.94686525000003, 50.880069750000018, 50.923958250000027))) %>%
list() %>%
sp::Polygons(ID="UKJ32 - Southampton")
I would now like to test whether the points in df
are in this polygon (and if so, return the Polygon ID).
tibble(lon = c(-1.4, 10), lat = c(50.9, 10))
Can someone tell me how I get to the result
tibble(lon = c(-1.4, 10), lat = c(50.9, 10), polyg_ID = 'UKJ32')
CodePudding user response:
If you wish to stick to sp
, there is a point.in.polygon()
function in sp
package:
UKJ32 <- sp::Polygon(cbind(c(-1.477037449999955, -1.366895449999959, -1.365159449999965, -1.477037449999955),
c(50.923958250000027, 50.94686525000003, 50.880069750000018, 50.923958250000027))) |>
list() |>
sp::Polygons(ID="UKJ32 - Southampton")
a <- tibble::tibble(lon = c(-1.4, 10), lat = c(50.9, 10))
sp::point.in.polygon(a$lon, a$lat, UKJ32@Polygons[[1]]@coords[,1], UKJ32@Polygons[[1]]@coords[,2])
#> [1] 1 0
Created on 2022-10-16 with reprex v2.0.2
CodePudding user response:
The {sp}
package is by now somewhat dated - after having lived a long & fruitful life - and most of current action happens in context of its successor, the {sf}
package.
Assigning some kind of a polygon feature - either an id or a metric - to a points dataset is a frequent use case. It at present often done via a sf::st_join()
call. For an example in action consider this earlier answer https://stackoverflow.com/a/64704624/7756889
I suggest that you try to move your workflow to the more current {sf}
package; you will find it easier to keep up with recent development.
And even if this were not possible for whatever reason - use sp::Polygons()
with utmost caution. I carries no information about coordinate reference system - which is a fancy way of saying it has no way of interpreting the coordinate numbers. Are they decimal degrees, or meters? Could be feet or fathoms for all that I know.
Strictly speaking you should not be allowed to proceed with a point-in-polygon operation calculation without this information.