I have the same city map (BH city) for two different shapefiles with different locations. The first one is the "shared management", this one have bigger polygons, as we can see

shared_management <- read_sf(dsn = 'TGC_POP_DENGUE_00_11_12_TAXA.shp')

enter image description here

And I have the smallers ones ("sectors")

sectors <- read_sf(dsn = 'SETOR_CENSITARIO_2010.shp')

enter image description here

Then I need a list of polygons of "sectors" which is inside of each polygon of "shared_management".

I know the number of polygons of "sectors" inside of which polygon of "shared_management". This number is given by

geo.int <- st_intersection(sector, shared_management)
geo.agg <- aggregate( ID_CS_2010~ID_TGC, data = geo.int, FUN = length)
geo.agg1 <- aggregate( ID_CS_2010~IDENT_TERR, data = geo.int, FUN = length)
1          B1        112
2          B2         98
3          B3         85
4          B4        123
5          B5        127

Any idea of how can I obtain which polygons of "sectors" are inside of each polygon of "shared_management"?


You can use the function st_join in the package sf to add the field you want from 'shared management' to the polygons within them that belong to 'sectors'

The code will look something like


geo_join <- st_join(shared_management, sectors, left = FALSE, largest = TRUE)

Setting the argument 'largest' to true will deal with overlap by assigning the sector the value from the shared_management zone that it overlaps with the most.

Once you have done this, you can filter the sector polygons by the field that indicates which of the larger polygons it fell within.

  • 1
    This is not the most stable solution. Please take a look at DE-9IM (Dimensionally-Extended 9 Intersection Model) behavior to understand how GEOS deals with overlay.intersect operations. Oct 30 '19 at 18:05
  • Thanks @JeffreyEvans, this is helpful, as I'd not thought about the math behind this. When you say that it's not the most stable, do you mean that it can produce differing results?
    – Kartograaf
    Oct 30 '19 at 18:20

Not the answer you're looking for? Browse other questions tagged or ask your own question.