Malawi_Project/01_Raster_Mosaicking_Buildings_2024.R at main · wpgp/Malawi_Project · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
library(terra)
library(sf)
library(tictoc)

# Specify data path
drive_path <- "//Working/MALAWI/Ortis/"
base_path <- paste0(drive_path, "Input_Data/") ## Base path where the folders are located
shp_path <- paste0(drive_path, "Input_Data/Shapefiles/") ## Shapefile path
result_path <- paste0(drive_path, "Input_Data/Mosaic_Buildings_2024/") # Result path
building_path <- paste0(drive_path, "Input_Data/Malawi_Covs/2024_Buildings/")

#Load data
boundary <- st_read(paste0(shp_path, "Country_Shapefile_Buffer_10km.shp"))
r1 <- rast(paste0(building_path, "mwi_buildings_count_2023_glv2_5_t0_5_C_100m_v1.tif"))

#Reproject boundary to r1
boundary <- st_transform(boundary, crs = st_crs(r1))

# Define folder names
folders <- c("Malawi_Covs/2024_Buildings","Mozambique_Covs/2024_Buildings",
             "Tanzania_Covs/2024_Buildings", "Zambia_Covs/2024_Buildings")


# Initialize a list to store raster file names
raster_files <- list()

tic()

# Loop through each folder and read all raster files  (.tif files)
for (folder in folders) {
  folder_path <- file.path(base_path, folder)
  files <- list.files(folder_path, pattern = "\\.tif$", full.names = TRUE)
  raster_files[[folder]] <- files
}

# Extract unique raster names (except the first three alphabets)
unique_raster_names <- unique(sapply(basename(unlist(raster_files)), function(x) substr(x, 4, nchar(x))))

# Function to process each unique raster name
process_raster <- function(raster_name) {
  rasters <- list()

  # Collect rasters with same name from each folder
  for (folder in folders) {
    matching_files <- raster_files[[folder]][sapply(basename(raster_files[[folder]]), function(x) substr(x, 4, nchar(x)) == raster_name)]
    rasters <- c(rasters, lapply(matching_files, rast))
  }

  # Get the CRS of the first raster (CMR)
  ref_crs <- crs(rasters[[1]])

  # Reproject all rasters to the CRS of the first raster
  rasters_reprojected <- lapply(rasters, function(r) project(r, ref_crs))

  # Mosaic the rasters together, prioritizing the first raster in case of overlap
  mosaic_raster <- do.call(mosaic, c(rasters_reprojected, fun = "first"))

  # Crop the mosaicked raster using the shapefile boundary
  cropped_raster <- crop(mosaic_raster, boundary)

  # Mask the cropped raster to the boundary
  masked_raster <- mask(cropped_raster, boundary)

  # Save the masked raster to a file with a name based on the original raster file name
  output_name <- paste0("MOS_MLW", raster_name)
  writeRaster(masked_raster, file.path(result_path, output_name), overwrite = TRUE)

  # Display a message after saving the raster
  message("Saved ", output_name)

  # Remove the rasters from memory
  rm(list = ls(pattern = "raster"))
  gc()
}

# Loop through each unique raster name and process
for (raster_name in unique_raster_names) {
  process_raster(raster_name)
}
toc()