Created
November 20, 2025 14:54
-
-
Save cszang/62ef84a4fba6fbadedc8666ca9405598 to your computer and use it in GitHub Desktop.
Extract lon/lat pairs from the monthly DWD raster files.
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| library(R.utils) | |
| library(terra) | |
| library(sp) | |
| # unzip a folder structure directly downloaded from the DWD CDC; all archives end with .gz | |
| recursive_unzip <- function(dir) { | |
| f <- list.files(dir, pattern = "\\.asc.gz", recursive = TRUE, full.names = TRUE) | |
| lapply(f, gunzip) | |
| } | |
| get_dwd_month <- function(index) { | |
| prepend0 <- c(rep("0", 9), rep("", 3)) | |
| paste(paste0(prepend0, 1:12), month.abb, sep = "_")[index] | |
| } | |
| get_file_pattern <- function() { | |
| "*\\.asc$" | |
| } | |
| get_parameter_base_path <- function(.base_path, .parameter_name, .month, | |
| flat = FALSE) { | |
| if (!flat) { | |
| file.path(.base_path, .parameter_name, .month) | |
| } else { | |
| file.path(.base_path, .parameter_name) | |
| } | |
| } | |
| get_gridfiles <- function(.base_path, .parameter_name, .month_index, flat = FALSE) { | |
| list.files( | |
| get_parameter_base_path( | |
| .base_path, .parameter_name, | |
| get_dwd_month(.month_index), | |
| flat = flat), | |
| pattern = get_file_pattern(), | |
| full.names = TRUE) | |
| } | |
| read_dwd_raster <- function(.raster_path) { | |
| .raster <- terra::rast(.raster_path) | |
| crs(.raster) <- crs("EPSG:31467") | |
| .raster | |
| } | |
| coords_from_lat_lon_old <- function(.lat, .lon) { | |
| coords <- data.frame(y = .lat, x = .lon) | |
| coordinates(coords) <- ~ x + y | |
| proj4string(coords) <- CRS("+init=epsg:4326") | |
| coords <- spTransform(coords, CRS("+init=epsg:31467")) | |
| } | |
| coords_from_lat_lon <- function(.lat, .lon, .crs) { | |
| .crs <- crs("EPSG:4326") | |
| coords <- terra::vect(matrix(c(.lon, .lat), ncol = 2, | |
| byrow = FALSE), crs = .crs) | |
| coords <- terra::project(coords, "EPSG:31467") | |
| coords | |
| } | |
| get_year_from_gridname <- function(.name) { | |
| as.numeric(substr(regmatches(.name, | |
| regexpr("([0-9]{6})\\.asc$", .name)), | |
| 1, 4)) | |
| } | |
| extract_lat_lon <- function(.lat, .lon, .parameter_name, .base_path, flat = FALSE) { | |
| if (length(.lat) != length(.lon)) { stop("no way") } | |
| m <- length(.lat) | |
| .coords <- coords_from_lat_lon(.lat, .lon) | |
| out <- data.frame() | |
| for (j in 1:12) { | |
| cat("month:", j, "\n", "year:") | |
| grids <- get_gridfiles(.base_path, .parameter_name, j, flat = flat) | |
| for (i in seq_along(grids)) { | |
| .year <- get_year_from_gridname(grids[i]) | |
| cat(.year, ".") | |
| .raster <- read_dwd_raster(grids[i]) | |
| outslice <- data.frame( | |
| lon = .lon, | |
| lat = .lat, | |
| year = .year, | |
| month = j, | |
| value = terra::extract(.raster, .coords)[,2] | |
| ) | |
| out <- rbind(out, outslice) | |
| } | |
| cat("\n\n") | |
| } | |
| out <- out[order(out$year, out$month), ] | |
| out | |
| } | |
| extract_lat_lon_flat <- function(.lat, .lon, .parameter_name, .base_path) { | |
| .coords <- coords_from_lat_lon(.lat, .lon) | |
| .j <- 1 | |
| grids <- get_gridfiles(.base_path, .parameter_name, .j, flat = TRUE) | |
| n <- length(grids) | |
| out <- data.frame( | |
| year = integer(n), | |
| month = integer(n), | |
| value = numeric(n) | |
| ) | |
| for (i in 1:n) { | |
| date <- stringr::str_extract(grids[i], "[0-9]{6}(?=\\.asc$)") | |
| out$year[i] <- as.numeric(substr(date, 1, 4)) | |
| out$month[i] <- as.numeric(substr(date, 5, 6)) | |
| .raster <- read_dwd_raster(grids[i]) | |
| out$value[i] <- terra::extract(.raster, .coords)[1,2] | |
| } | |
| out %>% arrange(year, month) | |
| } |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment