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?