Intersections Between Two Shapefiles in R

Question

I have these two shapefiles in R:

• file_1
• file_2

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?