17 - 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)

17.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")
## Reading layer `landkreise-in-germany' from data source `C:\Dropbox\Rpack\rdwd\misc\vign\landkreise-in-germany.shp' using driver `ESRI Shapefile'
## Simple feature collection with 403 features and 14 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: 5.866251 ymin: 47.27012 xmax: 15.04182 ymax: 55.05653
## geographic CRS: WGS 84
# intersections: list with msf rownumbers for each district:
int <- sf::st_intersects(lk, msf)
## although coordinates are longitude/latitude, st_intersects assumes that they are planar

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)

17.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"], # set dir if needed
                  res="monthly", var="kl", per="r", outvec=TRUE)
clims <- dataDWD(urls, varnames=FALSE, dir="../localdata")
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_DATUMHohenlohekreisBodenseekreisOffenbachBerlin
2017-12-1590.894.439.9 
2018-01-15125  115  67.8 
2018-02-1524.853.33.18
2018-03-1546.851.147.7 
2018-04-1537.811.735.6 
2018-05-1560.162.220.3