22 - use case: map climate data to Landkreise

Shapefile of Landkreis districts:
https://public.opendatasoft.com/explore/dataset/landkreise-in-germany/export/ (file size 4 MB, unzipped 10 MB)

22.1 find available meteo stations for each district

# Select monthly climate data:
data("metaIndex") ; m <- metaIndex
m <- m[m$res=="monthly" & m$var=="kl" & m$per=="recent" & m$hasfile, ]
# Transform into spatial object:
msf <- sf::st_as_sf(m, coords=c("geoLaenge", "geoBreite"), crs=4326)

# Read district shapefile, see link above:
lk <- sf::st_read("landkreise-in-germany.shp", quiet=TRUE)

# intersections: list with msf rownumbers for each district:
int <- sf::st_intersects(lk, msf)

https://gis.stackexchange.com/a/318629/36710

# plot to check projection:
plot(lk[,"id_2"], reset=FALSE)
colPoints("geoLaenge", "geoBreite", "Stationshoehe", data=m, add=T, legend=F)
# berryFunctions::colPointsLegend + sf plots = set margins, see note there!
axis(1, line=-1); axis(2, line=-1, las=1)
points(m[int[[2]], c("geoLaenge", "geoBreite")], pch=16, col=2, cex=1.8)

22.2 Average data per district

Running analysis for a few selected districts only to reduce computation time.
Monthly rainfall average per Landkreis.

landkreis_rain <- function(lki) # LandKreisIndex (row number in lk)
{
rnr <- int[[lki]] # msf row number
if(length(rnr)<1)
  {
  warning("No rainfall data available for Landkreis ", lki, ": ", lk$name_2[lki], call.=FALSE)
  out <- data.frame(NA,NA)[FALSE,]
  colnames(out) <- c("MESS_DATUM", as.character(lk$name_2[lki]))
  return(out)
  }
urls <- selectDWD(id=m[rnr, "Stations_id"], res="monthly", var="kl", per="r")
clims <- dataDWD(urls, varnames=FALSE)
if(length(urls)==1) 
  {rainmean <- clims$MO_RR 
  monthlyrain <- clims[c("MESS_DATUM", "MO_RR")]
  } else
{
monthlyrain <- lapply(seq_along(clims), function(n) 
 {
 out <- clims[[n]][c("MESS_DATUM", "MO_RR")]
 colnames(out)[2] <- names(clims)[n] # no duplicate names
 out
 })
monthlyrain <- Reduce(function(...) merge(..., by="MESS_DATUM",all=TRUE), monthlyrain)
rainmean <- rowMeans(monthlyrain[,-1], na.rm=TRUE) # check also with median, variation is huge!
}
out <- data.frame(monthlyrain[,1], rainmean)
colnames(out) <- c("MESS_DATUM", as.character(lk$name_2[lki]))
return(out)
}

rainLK <- pbapply::pblapply(c(133,277,300,389), landkreis_rain)

## Warning: No rainfall data available for Landkreis 300: Offenbach

rainLK <- Reduce(function(...) merge(..., by="MESS_DATUM",all=TRUE), rainLK)
head(rainLK)

MESS_DATUM	Hohenlohekreis	Bodenseekreis	Berlin
2021-04-15			29.2
2021-09-15	27.9	25.3	30
2021-10-15	56.2	30.4	20.8
2021-11-15	36.9	35.5	63.5
2021-12-15	93.5	94.5	44
2022-01-15	69.6	33.7	32.9