0

I'm trying to use st_within to find points within a polygon, but I can't get a shapefile into the right format. I'm using the standard shp file JPN_adm1.shp from https://purl.stanford.edu/np872wp5062 or https://gadm.org/download_country.html

JPN <- st_read("JPN_adm1.shp")
JPN <- JPN[1,] ## select one prefecture, Aichi
JPN <- st_transform(JPN, crs="WGS84")

tokyo.df <- data.frame("name"="Tokyo","lat"=35.507,"lon"=139.208)
some_point.sf <- st_as_sf(tokyo.df, coords = c("lon","lat"), crs="WGS84")
some_point.sf <- st_transform(some_point.sf, crs="WGS84")
subset <- st_join(some_point.sf, JPN, join = st_within)

The code runs, but it doesn't exclude Tokyo which isn't in Aichi. What's wrong?

1 Answer 1

2

st_join is doing a "left" join on the points to the polygon, so every point will have an output record. If its not within any polygon, you still get an output record but its all NAs:

> st_join(some_point.sf, JPN, join = st_within)
Simple feature collection with 1 feature and 13 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 139.208 ymin: 35.507 xmax: 139.208 ymax: 35.507
Geodetic CRS:  WGS 84
   name ID_0  ISO NAME_0 ID_1 NAME_1 HASC_1 CCN_1 CCA_1 TYPE_1 ENGTYPE_1
1 Tokyo   NA <NA>   <NA>   NA   <NA>   <NA>    NA  <NA>   <NA>      <NA>
  NL_NAME_1 VARNAME_1               geometry
1      <NA>      <NA> POINT (139.208 35.507)

If you want to drop unmatched records from the "left" side of the join, use left=FALSE:

> st_join(some_point.sf, JPN, join = st_within, left=FALSE)
Simple feature collection with 0 features and 13 fields
Bounding box:  xmin: NA ymin: NA xmax: NA ymax: NA
Geodetic CRS:  WGS 84
 [1] name      ID_0      ISO       NAME_0    ID_1      NAME_1    HASC_1   
 [8] CCN_1     CCA_1     TYPE_1    ENGTYPE_1 NL_NAME_1 VARNAME_1 geometry 
<0 rows> (or 0-length row.names)

and then there's no row for Tokyo.

Also worth checking with a point you know to be inside the polygon to see what happens:

> inA = data.frame("name"="pt", lat=35, lon=137.2)
> inA = st_as_sf(inA, coords=c("lon","lat"),crs="WGS84")

Then try the st_join function.

Often I find using the predicate function directly is more readable, so:

> st_within(some_point.sf, JPN)
Sparse geometry binary predicate list of length 1, where the predicate
was `within'
 1: (empty)

and then work on the returned list structure. But yes, often a spatial join is better, depends on where you are going with this.

Also note if you create points in a coordinate system you don't need to transform them as well - the second line here is redundant:

some_point.sf <- st_as_sf(tokyo.df, coords = c("lon","lat"), crs="WGS84")
some_point.sf <- st_transform(some_point.sf, crs="WGS84")
1
  • Thanks. This is perfect and the predicate function is exactly what I need to move the metadata. And having failed once to re-project spatial data, I do it somewhat redundantly!
    – Mark R
    Feb 4 at 20:21

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.