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
MESS_DATUM | Hohenlohekreis | Bodenseekreis | Offenbach | Berlin |
---|---|---|---|---|
2020-10-15 | 48.8 | 102 | 64.2 | |
2020-11-15 | 19.8 | 24 | 17.1 | |
2020-12-15 | 65.2 | 86.2 | 21.1 | |
2021-01-15 | 71.2 | 150 | 44.3 | |
2021-02-15 | 53.9 | 53.8 | 30.4 | |
2021-03-15 | 36.8 | 39.4 | 33.6 |