Maximiliano Micheli
Maximiliano Micheli

Reputation: 55

Intersection keeping non-intersecting polygons in sf - R

I'm looking to intersect 2 spatial layers, keeping all the non-intersecting features as well.

My first layer is the SA2 from NSW, Australia, which look like

enter image description here

My second layer is the Areas of Regional Koala Significance (ARKS):

enter image description here

When I intersect them, I get part of my desired result, which is subdividing the SA2 by the ARKS.

enter image description here

The thing is that I'd like to have also the rest of the SA2 polygons that don't intersect. The desired result would be a map of the SA2, where the intersecting ones would be subdivided by where they intersect to the ARKS layer, and the ones that don't intersect would contain NA. Something like in the next picture but in a single dataset instead of two: enter image description here

I post my code below:

library(sf)
SA2 <- st_read('C:/Users/Maxi/Desktop/NSW Planning/Koala/Data/SA2_2021_GDA94/SA2_2021_AUST_GDA94.shp')
ARKS <- st_read('C:/Users/Maxi/Desktop/ARKS_data/ARKS_renamedv6_PriorityPopulations_GDA94geog.shp')
ARKS <- st_transform(ARKS, 3577)
SA2 <- st_transform(SA2, 3577)
SA2 <- SA2[SA2$STE_CODE21=='1',]

intersection <- st_intersection(SA2, ARKS)

The data can be obtained from:

SA2: https://www.abs.gov.au/statistics/standards/australian-statistical-geography-standard-asgs-edition-3/jul2021-jun2026/access-and-downloads/digital-boundary-files#downloads-for-gda94-digital-boundary-files

ARKS: https://datasets.seed.nsw.gov.au/dataset/areas-of-regional-koala-significance-arks

Hope this is enough information!

Upvotes: 5

Views: 2217

Answers (1)

dieghernan
dieghernan

Reputation: 3402

Please consider this approach: Once that you have your intersection, you can remove the intersecting parts with st_difference. That would effectively split the intersecting SA2 in zones based on ARKS, and leave the rest as they are originally. After that, you can rejoin the dataset with dplyr::bind_rows, having the ARKS layer, the SA2 intersected split and the SA2 non-intersected as they are originally:

library(sf)
library(ggplot2)
library(dplyr)

SA2 <- st_read("./SA2_2021_AUST_SHP_GDA94/SA2_2021_AUST_GDA94.shp")
ARKS <-
  st_read(
    "./Biodiversity_KoalaPrioritisationProjectNSW_ARKS/KoalaPrioritisationProjectNSW_ARKS.shp"
  )
ARKS <- st_transform(ARKS, 3577)
SA2 <- st_transform(SA2, 3577)
SA2 <- SA2[SA2$STE_CODE21 == "1", ]



intersection <- st_intersection(SA2, ARKS)

# Control var, you can remove this
intersection$koala <- "Yes"

# Difference
SA2_diff <-
  st_difference(SA2, st_union(st_geometry(intersection)))

# Check visually, we should see the gaps
ggplot(SA2_diff) +
  geom_sf(fill = "green")

enter image description here


# Join all back
SA2_end <- dplyr::bind_rows(
  SA2_diff,
  intersection
)


ggplot(SA2_end) +
  geom_sf(aes(fill = koala))

enter image description here

Upvotes: 5

Related Questions