I'm using the sf package in R to understand how land cover changed from forest fires in BC, by overlapping burned areas (NBAC) with data from Dynamic World (DW). I have the DW data in the form of 5km by 5km grids with land cover proportions for each grid. When I use st_intersect(sf_2018, burn_2018), it duplicated four of the rows, with different geometries. Using st_intersection(sf_2018, st_difference(burn_2018)) does not resolve the problem, as suggested in this question. How can I avoid this?

Here is my code:

#Open Dynamic World Data
dw_2018 = read.csv("landcover_2018.csv")

##We can make this into a spatial object using the geoJSON information stored in .geo
sf_2018 = st_as_sf(data.frame(dw_2018, geom=geojson_sf(dw_2018$.geo)))

#Remove .geo column
sf_2018 = subset(sf_2018, select = -c(.geo))

#Open 2018 burned area data.
burn_2018 = st_read("nbac_2018.shp")

#Put a Thompson-Nicola Regional District (TNRD) shapefile into same CRS as burned areas
TNRD_conal = st_transform(TNRD, crs = crs(burn_2018))

#Subset burned areas to TNRD
burn_2018 = st_intersection(burn_2018, TNRD_conal)

#Now put the burn areas into UTM to match with the exported data
burn_2018 = st_transform(burn_2018, crs = crs(sf_2018))

#Find area of sample tiles
sf_2018$area = as.numeric(st_area(sf_2018))

#Subset grid data to only include grids that intersect polygons and the intersection area
##Note: problem still occurs when I use subset_18 = st_intersection(sf_2018, st_difference(burn_2018)) instead
subset_18 = st_intersection(sf_2018, burn_2018)

#Subset data so we only have the ID ("system.index")
subset_18 = subset_18[c(1)]

#check which rows of subset_18 are duplicated 
duplicated_rows <- subset_18[duplicated(subset_18$system.index), ]

#this number should be zero (there are four duplicated rows, they have different geometries but the same system.index)
duplicated_rows <- subset_18 %>%
  filter(system.index %in% c(duplicated_rows$system.index))

Here is the table of duplicated rows:

> duplicated_rows
Simple feature collection with 8 features and 1 field
Geometry type: GEOMETRY
Dimension:     XY
Bounding box:  xmin: -121.589 ymin: 50.24116 xmax: -119.7244 ymax: 51.99503
Geodetic CRS:  WGS 84
       system.index                       geometry
1670         2_3921 MULTIPOLYGON (((-120.2922 5...
2416         2_4770 POLYGON ((-119.7521 50.7734...
456          2_1766 POLYGON ((-121.5678 50.2530...
1738         2_3994 POLYGON ((-120.2048 51.9569...
1670.1       2_3921 POLYGON ((-120.2472 51.9927...
1738.1       2_3994 POLYGON ((-120.2454 51.9930...
2416.1       2_4770 POLYGON ((-119.7256 50.7660...
456.1        2_1766 POLYGON ((-121.5838 50.2439...
  • 1
    If you could make this code reproducible, it would greatly help others trying to fix the problem. Feb 27 at 15:14
  • Where does TNRD come from? The object suddenly appears without any prior assignments, and you did not reference the source. But most of all, we would need some subset of dw_2018, which you can e.g. create via head(dw_2018, 20) |> dput(). Right now, this is a black box and there is no starting point to help you.
    – dimfalk
    Mar 30 at 13:03


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Browse other questions tagged or ask your own question.