Merge branch 'dev' of https://github.com/KWB-R/kwb.rabimo into dev

Author: hsonne

R/natural-scenario.R

@@ -32,6 +32,9 @@ data_to_natural <- function(data, type = "undeveloped")
   nat_data <- data
   nat_data[urban_columns] <- 0
 
+  # vegetation class 50
+  nat_data["veg_class"] <- 50L
+
   if (type == "undeveloped") {
     return(nat_data)
   }
@@ -143,7 +146,7 @@ calculate_delta_w_1 <- function(
 calculate_delta_w_2 <- function(
     natural,
     urban,
-    columns_water_balance = c("surface_runoff", "infiltration", "evaporation"),
+    columns_water_balance = c("runoff", "infiltr", "evapor"),
     column_code = "code"
 )
 {
@@ -173,7 +176,7 @@ calculate_delta_w_2 <- function(
 calculate_delta_w_3 <- function(
     natural,
     urban,
-    columns_water_balance = c("surface_runoff", "infiltration", "evaporation"),
+    columns_water_balance = c("runoff", "infiltr", "evapor"),
     column_code = "code"
 )
 {

inst/extdata/column-names.csv

@@ -1,27 +1,27 @@
 rabimo_berlin,abimo_berlin,by_100,meaning,unit,type,data_type,default
 code,CODE,,Unique identifier of block,-,required,character,
 prec_yr,REGENJA,,Long-term average precipitation per year,mm,required,integer,600
-prec_s,REGENSO,,Long-term average precipitation per summer half year (May to October!),mm,required,integer,300
+prec_s,REGENSO,,Long-term average precipitation per summer half year (May to October),mm,required,integer,300
 epot_yr,,,Long-term average potential evapotranspiration per year,mm,required,integer,700
-epot_s,,,Long-term average potential evapotranspiration per summer half year (May to October!),mm,required,integer,500
-district,BEZIRK,,Specific to Berlin: name of city district,-,required?,character,0
+epot_s,,,Long-term average potential evapotranspiration per summer half year (May to October),mm,required,integer,500
+district,BEZIRK,,Specific to Berlin: identifier of city district,-,required?,character,0
 total_area,,,Total block area,m2,required,numeric,100
 area_main,FLGES,,Total non-road area within the block,m2,required?,numeric,100
 area_road,STR_FLGES,,Total road area within the block,m2,required?,numeric,0.0
 main_frac,,,Non-road fraction of total block area (area_main/total_area),0..1,required,numeric,
 roof,PROBAU,x,Roof fraction of non-road built area,0..1,required,numeric,0.2
 green_roof,,,Green roof fraction of roof area,0..1,required,numeric,0.0
-swg_roof,KAN_BEB,x,Fraction of ... ,0..1,required,numeric,1.0
-pvd,PROVGU,x,Fraction of ...,0..1,required,numeric,0.6
-swg_pvd,KAN_VGU,x,Fraction of ...,0..1,required,numeric,0.7
+swg_roof,KAN_BEB,x,Fraction of roof area connected to the sewer ,0..1,required,numeric,1.0
+pvd,PROVGU,x, Paved fraction of non-road area,0..1,required,numeric,0.6
+swg_pvd,KAN_VGU,x,Fraction of paved area connected to the sewer,0..1,required,numeric,0.7
 srf1_pvd,BELAG1,x,Fraction of paved area belonging to surface class 1,0..1,required,numeric,0.5
 srf2_pvd,BELAG2,x,Fraction of paved area belonging to surface class 2,0..1,required,numeric,0.2
 srf3_pvd,BELAG3,x,Fraction of paved area belonging to surface class 3,0..1,required,numeric,0.1
 srf4_pvd,BELAG4,x,Fraction of paved area belonging to surface class 4,0..1,required,numeric,0.1
 srf5_pvd,BELAG5,x,Fraction of paved area belonging to surface class 5,0..1,required,numeric,0.1
 road_frac,,,Road fraction of total block area (area_road/total_area),0..1,required,numeric,0.0
-pvd_r,VGSTRASSE,x,Fraction of ...,0..1,required,numeric,0.9
-swg_pvd_r,KAN_STR,x,Fraction of ...,0..1,required,numeric,1.0
+pvd_r,VGSTRASSE,x,Paved fraction of road area,0..1,required,numeric,0.9
+swg_pvd_r,KAN_STR,x,Fraction of paved road area connected to the sewer,0..1,required,numeric,1.0
 srf1_pvd_r,STR_BELAG1,x,Fraction of road area belonging to surface class 1,0..1,required,numeric,0.9
 srf2_pvd_r,STR_BELAG2,x,Fraction of road area belonging to surface class 2,0..1,required,numeric,0.1
 srf3_pvd_r,STR_BELAG3,x,Fraction of road area belonging to surface class 3,0..1,required,numeric,0.0

inst/scripts/epot-raster-to-district.R

@@ -0,0 +1,157 @@
+library(terra)
+library(stringr)
+library(R.utils)
+
+# Set the directory where files are stored
+data_dir <- "C:/Users/fdpunt/Downloads/dwd/grids_germany/monthly/evapo_p"
+output_dir <- "C:/Users/fdpunt/Documents/Projekte/AMAREX/Daten/Raster/epot"
+unzipped_dir <- file.path(data_dir, "unzipped")
+
+# download and save data
+if(FALSE) {
+  from <- "199101"
+  to <- "202312"
+  pot_evaporation_monthly <- kwb.dwd::load_potential_evaporation_berlin(from, to)
+}
+
+# List all .asc.gz files
+files <- list.files(data_dir, pattern = "grids_germany_monthly_evapo_p_\\d{6}\\.asc\\.gz$", full.names = TRUE)
+file_names <- basename(files)
+
+grid <- kwb.dwd::read_asc_gz_file(files[1])
+
+# Extract year and month information
+timestamp <- str_extract(file_names, "\\d{6}")
+year <- as.numeric(substr(timestamp, 1, 4))
+month <- as.numeric(substr(timestamp, 5, 6))
+
+# Create the data frame
+file_info <- data.frame(
+  file = files,
+  year = year,
+  month = month
+)
+
+districts <- sf::read_sf("C:/Users/fdpunt/Documents/Projekte/AMAREX/Daten/Bezirke/Bezirke/ISU5_UA_FL_2020_Bezirke.shp")
+
+# Function to check and handle raster values before saving
+process_raster <- function(raster, scale = 0.1) {
+
+  raster <- raster * scale
+  return(raster)
+}
+
+# Process yearly total for all months
+years <- 1991:2020
+raster_list_year <- list()
+raster_list_summer <- list()
+
+# Loop over each year
+for (i in seq_along(years)) {
+  # Get files for the current year
+  yr <- years[i]
+  year_files <- file_info$file[file_info$year == yr]
+  summer_files <- file_info$file[file_info$year == yr &
+                                   file_info$month %in% 5:10]
+
+  # Read all monthly rasters for the current year
+  rasters_year <- lapply(year_files, kwb.dwd::read_asc_gz_file)
+  rasters_year_stack <- raster::stack(rasters_year)
+
+  # Read all summer rasters for the current year
+  rasters_summer <- lapply(summer_files, kwb.dwd::read_asc_gz_file)
+  rasters_summer_stack <- raster::stack(rasters_summer)
+
+  # Sum rasters to get the yearly total
+  yearly_sum <- sum(rasters_year_stack, na.rm = TRUE)
+  summer_sum <- sum(rasters_summer_stack, na.rm = TRUE)
+
+  # Process raster (data scaling)
+  yearly_sum <- process_raster(yearly_sum)
+  summer_sum <- process_raster(summer_sum)
+
+  # Save raster to list
+  raster_list_year[[as.character(yr)]] <- yearly_sum
+  raster_list_summer[[as.character(yr)]] <- summer_sum
+  print(paste("Stored rasters for year:", yr))
+
+  # Define the output file path
+  output_file_year <- file.path(output_dir, paste0("evapo_p_total_", yr, ".tif"))
+  output_file_summer <- file.path(output_dir, paste0("evapo_p_summer_", yr, ".tif"))
+
+  # Save the raster to a .tif file
+  terra::writeRaster(yearly_sum, output_file_year, overwrite = TRUE)
+  terra::writeRaster(summer_sum, output_file_summer, overwrite = TRUE)
+
+  # visualize the result (for verification)
+  plot(yearly_sum, main = paste("Yearly Sum of Potential Evaporation", yr))
+  plot(summer_sum, main = paste("Summer Sum of Potential Evaporation", yr))
+}
+
+# Calculate average raster
+years_stack <- raster::stack(raster_list_year)
+average_epot_year <- mean(rast(years_stack), na.rm = TRUE)
+
+summers_stack <- raster::stack(raster_list_summer)
+average_epot_summer <- mean(rast(summers_stack), na.rm = TRUE)
+
+# Reproject the shapefile to match the raster's CRS
+average_epot_crs <- sf::st_crs(average_epot_year)$proj4string
+berlin_reprojected <- sf::st_transform(districts, average_epot_crs)
+
+# Save the rasters
+output_file <- file.path(output_dir, paste0("yearly-average_", min(years), "-", max(years), ".tif"))
+terra::writeRaster(average_epot_year, output_file, overwrite = TRUE)
+
+output_file <- file.path(output_dir, paste0("summer-average_", min(years), "-", max(years), ".tif"))
+terra::writeRaster(average_epot_summer, output_file, overwrite = TRUE)
+
+# Visualize the raster with the reprojected shapefile
+plot(average_epot_year, , main = paste0("average epot ", min(years), "-", max(years)))
+plot(average_epot_summer, , main = paste0("average summer epot ", min(years), "-", max(years)))
+
+# mask and crop Berlin
+berlin_epot_avg_year <- (raster::mask(average_epot_year, berlin_reprojected[,1]))
+berlin_epot_avg_year_crop <- raster::crop(berlin_epot_avg_year, berlin_reprojected[,1])
+plot(berlin_epot_avg_year_crop,
+     main = "mean yearly potential evaporation berlin", axes = NA)
+plot(berlin_reprojected[,1], col = NA, border = "white", add = TRUE)
+district_centroids <- sf::st_centroid(berlin_reprojected)
+text(sf::st_coordinates(district_centroids),
+     labels = berlin_reprojected$bezneu, col = 'white', cex = 0.4)
+
+berlin_epot_avg_summer <- (raster::mask(average_epot_summer, berlin_reprojected[,1]))
+berlin_epot_avg_summer_crop <- raster::crop(berlin_epot_avg_summer, berlin_reprojected[,1])
+plot(berlin_epot_avg_summer_crop,
+     main = "mean summer potential evaporation berlin", axes = NA)
+plot(berlin_reprojected[,1], col = NA, border = "white", add = TRUE)
+district_centroids <- sf::st_centroid(berlin_reprojected)
+text(sf::st_coordinates(district_centroids),
+     labels = berlin_reprojected$bezneu, col = 'white', cex = 0.4)
+
+# Convert the sf object to a SpatVector
+districts_spatvector <- terra::vect(berlin_reprojected)
+
+# Extract mean raster values for each district
+district_means_year <- terra::extract(berlin_epot_avg_year_crop,
+                                 districts_spatvector,
+                                 fun = mean,
+                                 na.rm = TRUE)
+
+district_means_summer <- terra::extract(berlin_epot_avg_summer_crop,
+                                      districts_spatvector,
+                                      fun = mean,
+                                      na.rm = TRUE)
+
+epot_berlin_1991_2020 <- berlin_reprojected %>%
+  dplyr::select(ID = bezneu, district = bezirk) %>%
+  dplyr::mutate(ID = as.integer(ID)) %>%
+  sf::st_drop_geometry() %>%
+  dplyr::left_join(district_means_year, by = "ID") %>%
+  dplyr::mutate(total_mean = as.integer(round(mean))) %>%
+  dplyr::select(-mean) %>%
+  dplyr::left_join(district_means_summer, by = "ID") %>%
+  dplyr::mutate(summer_mean = as.integer(round(mean))) %>%
+  dplyr::select(-mean)
+
+

inst/scripts/test-rabimo.R

@@ -55,7 +55,7 @@ if (FALSE)
 
   rabimo_inputs_2020 <- kwb.rabimo::prepare_berlin_inputs(
     data = berlin_2020_data,
-    config = kwb.abimo::read_config()
+    config = kwb.abimo:::read_config()
   )
 
   # Fehler: Column 'gw_dist' must not contain missing values (NA, found 4