How to Find Which Polygon a Point Belong to via Sf

How to find which polygon a point belong to via sf

Use st_point() on the lon/lat then it can work with other sf functions.

Example:

find_precinct <- function(precincts, point) {
  precincts %>%
    filter(st_contains(geometry, point) == 1) %>%
    `[[`("WARDS_PREC")
}

ggmap::geocode("nicollet mall, st paul") %>%
  rowwise() %>%
  mutate(point = c(lon, lat) %>%
           st_point() %>%
           list(),
         precinct = find_precinct(msvc_precincts, point)
         )

finding points with polygon in sf package

I think you are looking for st_filter. This way you'll get only the points that are found within the polygon. The good points object now contains only two points, instead of three (like in the OP).

good_points <- st_filter(point.sf, poly)

# Simple feature collection with 2 features and 0 fields
# Geometry type: POINT
# Dimension:     XY
# Bounding box:  xmin: 136.2 ymin: 36.5 xmax: 136.5 ymax: 36.5
# CRS:           NA
#            geometry
# 1 POINT (136.2 36.5)
# 2 POINT (136.5 36.5)

find the nearest polygon for a given point

If you're ok to convert to sf objects, you can find the nearest polygon to each point outside the polys using something on this lines:

set.seed(1)
library(raster)
library(rgdal)
library(rgeos)
library(sf)
library(mapview)
  
p <- shapefile(system.file("external/lux.shp", package="raster"))
p2 <- as(1.5*extent(p), "SpatialPolygons")
proj4string(p2) <- proj4string(p)
pts <- spsample(p2, n=10, type="random")
  
## Plot to visualize
plot(p, col=colorRampPalette(blues9)(12))
plot(pts, pch=16, cex=.5,col="red", add = TRUE)

# transform to sf objects
psf   <- sf::st_as_sf(pts) %>% 
    dplyr::mutate(ID_point = 1:dim(.)[1])
polsf <- sf::st_as_sf(p)

# remove points inside polygons
in_points  <- lengths(sf::st_within(psf,polsf))
out_points <- psf[in_points == 0, ]

# find nearest poly
nearest <- polsf[sf::st_nearest_feature(out_points, polsf) ,]  %>% 
    dplyr::mutate(id_point = out_points$ID)
nearest

> Simple feature collection with 6 features and 6 fields
> geometry type:  POLYGON
> dimension:      XY
> bbox:           xmin: 5.810482 ymin: 49.44781 xmax: 6.528252 ymax: 50.18162
> geographic CRS: WGS 84
>   ID_1       NAME_1 ID_2           NAME_2 AREA                       geometry id_point
> 1    2 Grevenmacher    6       Echternach  188 POLYGON ((6.385532 49.83703...        1
> 2    1     Diekirch    1         Clervaux  312 POLYGON ((6.026519 50.17767...        2
> 3    3   Luxembourg    9 Esch-sur-Alzette  251 POLYGON ((6.039474 49.44826...        5
> 4    2 Grevenmacher    7           Remich  129 POLYGON ((6.316665 49.62337...        6
> 5    3   Luxembourg    9 Esch-sur-Alzette  251 POLYGON ((6.039474 49.44826...        7
> 6    2 Grevenmacher    6       Echternach  188 POLYGON ((6.385532 49.83703...        9
> 

#visualize to check
mapview::mapview(polsf["NAME_2"]) + mapview::mapview(out_points)

HTH !

Is a POINT within a MULTIPOLYGON using sf in R

@Tony Machado, you are running into a resolution issue here. The borderline of the polygon (i.e. Canada) from naturalearth comes in a resolution that may cut the borderline.
If you plot your problem at hand, you will see that some of the cities are close to the border.

library(ggplot2)
ggplot() + 
  geom_sf(data = can, fill = "lightblue") + 
  geom_sf(data = addresses_sf, color = "red", size = 2)

You can either

• get a higher resolution boundary (e.g. add a scale = "large" in your ne_countries() call); or work
• with a buffer.

For the buffering it might be worth reprojecting your WGS84 to UTM and work with a finer distance measure. In WGS84 the dist is in arc degrees. That can be painful.
But from below, you can see that you increase your includes when we buffer the Canadian borderline. Unfortunately, this catches also "Detroit" :(

Hope this gets you going.

sf::st_intersection(addresses_sf, st_buffer(can, dist = 0.001))
Simple feature collection with 5 features and 3 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: -83.1107 ymin: 42.11125 xmax: -70.84368 ymax: 45.41126
Geodetic CRS:  WGS 84
# A tibble: 5 x 4
  rowid country name               geometry
* <int> <chr>   <chr>           <POINT [°]>
1     2 canada  Canada (-75.70591 45.41126)
2     3 canada  Canada   (-83.06181 42.235)
3     4 canada  Canada (-70.84368 45.38441)
4     5 canada  Canada  (-83.1107 42.11125)
5     6 usa     Canada (-74.86597 44.74441)

Checking if a point falls within polygon Shapefile

I have my own approach. Would that fulfill your requirements? I can't tell you what specifically is wrong with your code, but this one is also a bit cleaner:

library(sf)
tt <- read_sf('./Downloads/taxi_zones/taxi_zones.shp')

pnts <- data.frame(
  "x" = c(-73.97817, -74.00668, 0, 500),
  "y" = c(40.75798, 40.73178, 0, 400)
)

pnts_sf <- st_as_sf(pnts, coords = c('x', 'y'), crs = st_crs(4326))

pnts_trans <- st_transform(pnts_sf, 2163)
tt_trans <- st_transform(tt, 2163)

pnts_trans <- pnts_sf %>% mutate(
  intersection = as.integer(st_intersects( pnts_trans,tt_trans)))

The result would be

                    geometry intersection
1 POINT (-73.97817 40.75798)          161
2 POINT (-74.00668 40.73178)          158
3                POINT (0 0)           NA
4            POINT (500 400)           NA

point in polygon with `sf`

The question as posted is somewhat unclear - what exactly is the problem?

The way I read it is that you have the point-in-polygon problem already solved (via the sf::st_join() call). The id of each polygon is stored in the ID column of your points_sf_joinded object (the point id is in id column, which may be kind of confusing btw.)

See the result of your last plot, slightly amended by me (I have removed the plotting of original coordinates not assigned to a polygon, and color coded the polygon IDs).

I have also added a geom_sf() call for the polygons to be plotted.

ggplot() +
  geom_sf(data = buffers) + # plot the polygons first
  geom_sf(data=points_sf_joined, aes(color = as.factor(ID))) +
  scale_color_brewer("polygon ID", palette = "Dark2") + 
  theme_minimal()

Check whether point coordinate lies within polygon

The class of problems you describe is called point in polygon.

You can handle it via sf::st_join() - it will add polygon columns to your point dataset. It uses left join as default = preserves rows of left hand object. It is therefore the best to start with your points object, and add to that characteristics of your polygons when spatially aligned (and NAs when not).

Two things to keep in mind:

• the CRS system has to be aligned (which one you choose matters rarely; it would take an extreme edge case - but it has to be aligned; use st_transform for that)
• depending on the class of your points object (sfc vs. sf) you may need to call st_as_sf() first

Consider something along these lines:

pip <- points %>%
  st_join(heat_df)

R sf point Popups inaccessible after adding polygon via simplevis

Use addPolylines instead of addPolygons, then the polygon feature won't be capturing the click events (which is what I think is happening) because its only lines instead of an area-covering polygon.

---

Related Topics