1

I have these two shapefiles in R:

I imported these files into R:

file_1 <- sf::st_read("C:/Users/me/OneDrive/Documents/folder1/lada000b21a_e.shp", options = "ENCODING=WINDOWS-1252")
file_2 <- sf::st_read("C:/Users/me/OneDrive/Documents/folder2/lpc_000b21a_e.shp", options = "ENCODING=WINDOWS-1252")

file_1 <- st_transform(file_1, crs = "+proj=longlat +datum=WGS84")
file_2 <- st_transform(file_2, crs = "+proj=longlat +datum=WGS84")


> file_1

Simple feature collection with 1507 features and 4 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -95.15386 ymin: 41.68132 xmax: -74.34347 ymax: 56.05945
CRS:           +proj=longlat +datum=WGS84
First 10 features:
       ADAUID             DGUID LANDAREA PRUID                       geometry
1567 35010001 2021S051635010001 643.4007    35 MULTIPOLYGON (((-74.48809 4...
1568 35010002 2021S051635010002 605.0164    35 MULTIPOLYGON (((-74.55843 4...
1569 35010003 2021S051635010003 515.4641    35 MULTIPOLYGON (((-74.90049 4...

>file_2
> pop_mini
Simple feature collection with 310 features and 9 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -116.9893 ymin: 41.98072 xmax: -64.09933 ymax: 53.53916
CRS:           +proj=longlat +datum=WGS84
First 10 features:
   PCUID PCPUID         DGUID DGUIDP                 PCNAME PCTYPE PCCLASS LANDAREA PRUID                       geometry
1   0001 350001 2021S05100001   <NA>                  Acton      2       2   7.8376    35 MULTIPOLYGON (((-80.00597 4...
5   0006 350006 2021S05100006   <NA>             Alexandria      4       2   2.0689    35 MULTIPOLYGON (((-74.63831 4...
6   0007 350007 2021S05100007   <NA>                 Alfred      4       2   0.8654    35 MULTIPOLYGON (((-74.87997 4...

My Question: I am trying to construct a matrix which shows what percent of each polygon in file_1 is covered by polygons in file_2 - and what percent of each polygon in file_2 is covered by polygons in file_1.

Based on some research I did (e.g. R: Checking Matrix Sums and invalid sf geometries ), I first repaired the geometries with both of these files:

library(lwgeom)
library(sf)
library(dplyr)

# Repair invalid geometries in file_1
file_1$geometry <- st_make_valid(file_1$geometry)

# Repair invalid geometries in file_2
file_2$geometry <- st_make_valid(file_2$geometry)

Then, I tried to write a matrix-loop procedure to calculate the percent coverage of pairwise polygons in both files:

# Calculate the number of polygons in each file
n_file_1 <- nrow(file_1)
n_file_2 <- nrow(file_2)

# Create a matrix to store coverage percentages
coverage_matrix <- matrix(0, n_file_1, n_file_2)

# Calculate the number of polygons in each file
n_ada <- nrow(file_1)
n_pop <- nrow(file_2)

# Create a matrix to store coverage percentages
coverage_matrix <- matrix(0, n_file_1, n_file_2)

# Calculate coverage percentages for each pair of polygons
for (i in seq_len(n_file_1)) {
    for (j in seq_len(n_file_2)) {
        intersection_area <- st_area(st_intersection(file_1[i,], file_2[j,]))
        if (length(intersection_area) > 0) {
            file_1_area <- st_area(file_1[i,])
            coverage_matrix[i,j] <- 100 * intersection_area / file_1_area
        }
        # Print intermediate results
        cat(paste0("i: ", i, ", j: ", j, ", coverage: ", coverage_matrix[i,j], "\n"))
    }
}

# Set row and column names for the coverage matrix
rownames(coverage_matrix) <- paste0("file_1 ", seq_len(n_file_1))
colnames(coverage_matrix) <- paste0("file_2 ", seq_len(n_file_2))

# Print the coverage matrix
print(coverage_matrix)

The code appears to be running :

i: 1, j: 1, coverage: 0
i: 1, j: 2, coverage: 0.349480438105992
i: 1, j: 3, coverage: 0
i: 1, j: 4, coverage: 0

But I am not sure if I am doing this correctly.

How can I do this?

4
  • Is this different to your previous question? gis.stackexchange.com/questions/456175/…
    – Spacedman
    Commented Mar 27, 2023 at 7:20
  • @ spacedman: thank you for your reply! I deleted that question and made a new one that is more reproducible
    – stats_noob
    Commented Mar 27, 2023 at 7:31
  • 1
    The way to progress here is to a) write your code as a function so you can run it easily on multiple data sets, and b) create some simple data examples where you can run and verify by hand it gets the right results. Try making some rectangular polygons with integer coordinates with different degrees of overlap and that should be a good test data set. Time spent preparing tests like that will save you time in the end.
    – Spacedman
    Commented Mar 27, 2023 at 12:38
  • @ spacedman: thank you for your reply! I tried to make a test over here : stackoverflow.com/questions/75849216/…
    – stats_noob
    Commented Mar 29, 2023 at 13:51

0

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.