I am new to geodata so maybe this is a silly question. I have two multipolygons and what I want to do is to test in which of the polygon of "constituencies" each polygon of "wards" is contained. Basically, to see for each ward in which consituency it is in.

However, I am having issues with the function st_within() from sf because it only returns TRUE if the first geometry is completely within the second geometry. That is a problem when the smaller units touch the borders of the larger ones. Here is an example

I have one multipolygon called "consituencies" and and one called "wards". You can find them here: https://www.dropbox.com/sh/xu08wm79rym00zz/AADFDpyPe0EuDSDSY-6SvqpYa?dl=0

There's 2 objects in constituencies:


enter image description here I then try to join them on whether the wards are within the constituencies

st_crs(wards) == st_crs(constituencies)#they have different CRS
test98 <- st_join(wards, constituencies, join=st_within)

And the outcome spits out many NAs on the constituencies' name and label a ward is within of.

Simple feature collection with 6 features and 4 fields
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: 531936.6 ymin: 180569.5 xmax: 533595.1 ymax: 182082.5
projected CRS:  OSGB 1936 / British National Grid
  WD98CD       WD98NM                             name label                       geometry
1 00AAFA   Aldersgate                             <NA>  <NA> MULTIPOLYGON (((532104.9 18...
2 00AAFB      Aldgate                             <NA>  <NA> MULTIPOLYGON (((533319.2 18...
3 00AAFC    Bassishaw Cities of London and Westminster   107 MULTIPOLYGON (((532587.3 18...
4 00AAFD Billingsgate Cities of London and Westminster   107 MULTIPOLYGON (((533167.9 18...
5 00AAFE  Bishopsgate                             <NA>  <NA> MULTIPOLYGON (((533410.7 18...
6 00AAFF Bread Street Cities of London and Westminster   107 MULTIPOLYGON (((532300.3 18...

When I plot them to see which ones have NA in name (in red), you can see it is the ones touching the consituency's borders:

plot(test98$geometry[!is.na(test98$name),], col = "green", add = T)
plot(test98$geometry[is.na(test98$name),], col = "red", add = T)

enter image description here

I assume one of two things is happening: 1) touching borders dos not count as being completely within a geometry, or 2) the borders of the wards and constituencies may not be perfectly matched so the wards aren't entirely within the constituency.

Is there a way of perfectly match the borders and then test which wards are within which constituencies OR is there a way of testing within which constituencies is not the entirety but the majority of a ward's boundaries?

  • Have answered in your cross-post on stackoverflow. R GIS q's are probably better here though.
    – Spacedman
    Jun 4 at 11:34
  • The "within" relationship is formally defined as "union of the interiors is the containing shape and the intersection of the exteriors is empty" so it looks like the function is behaving correctly.
    – Vince
    Jun 4 at 11:48
  • @Spacedman just seen it, thanks! Good to know they're better here. I was unsure but next time I'll know
    – AntVal
    Jun 4 at 11:50
  • @Vince, thanks! yes. I know it is. My issue is what should I do instead in order to get what I am looking for. Because neither st_within or other commans like intersects get right at what I am trying to get
    – AntVal
    Jun 4 at 11:51
  • @Spacedman answered your question --- you need to fix your data.
    – Vince
    Jun 4 at 11:56

As discussed on the StackOverflow question, this is because the geometries do not quite match exactly. If they did match exactly (down to the vertex level), ST_Within should produce the desired answer.

Another way to do this is to test whether an interior point of a ward is within a constituency. Interior points can be generated using ST_PointOnSurface. This is likely faster than testing with the full polygon, as well.

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.