Meteorological Review
Victor Munoz S.
2023-09-08
1 Introduction
SRK Consulting was retained by xxxx to prepare a review of the future meteorological condition to be expected at the site xxxx.
SRK has already prepared a meteorological evaluation at the site including site analogue daily records for total precipitation and mean temperature.
Longitude= -80.6448
Latitude = 8.8520
source(here::here("scripts","hydro_support.r"))
#For more information
#https://developers.google.com/maps/gmp-get-started\
key_google<-rstudioapi::askForSecret("Google Key")
#For more information
# https://www.ncdc.noaa.gov//cdo-web/token
key_noaa<-rstudioapi::askForSecret("NOAA Key")
ggmap::register_google(key = key_google)
#Get elevation from google maps
Elevation=GM_elev(Longitude,Latitude,key_google)
site.info<-data.frame(Location="Site",
Longitude=Longitude,
Latitude=Latitude,
Elevation=Elevation)
#Site in looped with variables
site.info_xy<-site.info %>%
dplyr::rename(x="Longitude",y="Latitude",z="Elevation") %>%
reshape2::melt(id.vars="Location")
pandoc.table(site.info)| Location | Longitude | Latitude | Elevation |
|---|---|---|---|
| Site | -80.64 | 8.852 | 141.9 |
2 SOURCES
2.1 NOAA
Local meteorological station records from publicly available National Oceanic and Atmospheric Administration (NOAA) databases: Global Historical Climatology Network Database (GHCND) and Global Surface Summary of the Day (GSOD) were obtained to implement the site data. These stations typically compile precipitation and temperature daily records. A radius of about ~ 100 km around the site was used as an area of interest to compile these meteorological records.
From NOAA, the amount of information in the region of interest is limited. Most of the relevant records present a gap of information after 1990. Meteorological stations with records less than 10years were removed. Records from NOAA have been complemented with Land precipitation and temperature records.
NOAA.file<-here::here("data","NOAA","NOAA Info.rdata")
if(file.exists(NOAA.file)){
load(NOAA.file)
}else{
NOAA.find(longitude = Longitude,latitude = Latitude,radius_deg = 1.5,
#Get your own key from https://www.ncdc.noaa.gov//cdo-web/token
#Key from vmunoz@srk.co.uk
noaa.key = "PacdvUWVYvHOPwXrHRILLOJpvDnORxEB")
#correction from DB
selected.NOAA.GHCND$maxdate<-if_else(selected.NOAA.GHCND$maxdate>ymd("2020-01-01"),
as.Date(lubridate::now()),selected.NOAA.GHCND$maxdate)
selected.NOAA.GSOD$END<-if_else(ymd(selected.NOAA.GSOD$END)>ymd("20200101"),
str_remove_all(as.Date(lubridate::now()),"-") %>% as.integer(),
selected.NOAA.GSOD$END)
#work just with the first 10 stations. Ranges can change based on the project
# selected.NOAA.GHCND<-selected.NOAA.GHCND[1:10,]
# selected.NOAA.GSOD<-selected.NOAA.GSOD[1:10,]
NOAA.capture.GHCND()
NOAA.capture.GSOD()
NOAA.prepare()
save(file = NOAA.file,list = c("NOAA.db","NOAA.idx"))
}
NOAA.idx %>%
dplyr::select(-WMO_ID,-ICAO) %>%
pandoc.table(style="simple",split.table=Inf,
caption="NOAA STATIONS")| ID | NAME | LATITUDE | LONGITUDE | ELEVATION | MINDATE | MAXDATE | DISTANCE |
|---|---|---|---|---|---|---|---|
| PM000105001 | BOCA DEL TOABRE | 8.91 | -80.55 | 170 | 1979-01-01 | 1993-12-31 | 12.26 |
| PM000132001 | EL PALMAR | 8.53 | -81.06 | 1000 | 1979-01-01 | 1993-12-31 | 58.07 |
| 749035_99999 | AGUADULCE | 8.25 | -80.57 | 46 | 1942-03-30 | 1946-01-25 | 67.56 |
| 787957_99999 | COCLE/CAP SCARLETT MARTINEZ ARPT/RIO HATO | 8.383 | -80.13 | 858.9 | 2005-10-02 | 2023-09-08 | 76.8 |
| 788078_99999 | RIO HATO | 8.367 | -80.12 | 32 | 1939-01-17 | 1989-12-19 | 79.31 |
| 749025_99999 | CHAME | 8.6 | -79.92 | 37 | 1942-07-14 | 1946-10-24 | 84.85 |
| 699404_99999 | FT SHERMAN (ROCOB) | 9.333 | -79.98 | NA | 1980-04-23 | 1981-06-12 | 90.33 |
| 788010_99999 | FT SHERMAN (ROCOB) | 9.333 | -79.98 | 52 | 1973-01-02 | 2006-01-02 | 90.33 |
| 787950_99999 | RUBEN CANTU | 8.086 | -80.94 | 82.9 | 1974-01-01 | 2023-09-08 | 91.45 |
| 749027_99999 | LA MITRA | 8.85 | -79.78 | 94 | 1942-04-08 | 1945-10-19 | 94.79 |
| PMW00010719 | FT SHERMAN | 9.367 | -79.95 | 9.1 | 1965-01-01 | 1971-12-31 | 95.47 |
| 999999_10719 | FORT SHERMAN | 9.367 | -79.95 | 9.1 | 1965-04-01 | 1971-07-01 | 95.49 |
| 749032_99999 | LAKE GATUN | 9.117 | -79.78 | 79 | 1945-05-01 | 1946-01-31 | 99.24 |
| PMW00010715 | COCO SOLO | 9.367 | -79.9 | 4.3 | 1945-01-01 | 1957-12-31 | 99.92 |
| 999999_10715 | COCO SOLO | 9.367 | -79.9 | 4 | 1945-05-01 | 1957-07-18 | 99.94 |
| 787830_99999 | ENRIQUE ADOLFO JIMENEZ | 9.357 | -79.87 | 7.6 | 1951-01-01 | 2019-04-14 | 102.3 |
| PMW00010705 | HOWARD | 8.917 | -79.6 | 10.1 | 1949-01-01 | 1949-12-31 | 115.1 |
| PMW00010718 | FT KOBBE | 8.917 | -79.6 | 16.2 | 1954-01-01 | 1967-12-31 | 115.1 |
| 787955_99999 | PANAMA PACIFICO | 8.917 | -79.6 | 13 | 2006-10-10 | 2023-09-08 | 115.1 |
| 788060_99999 | HOWARD | 8.917 | -79.6 | 13 | 1941-10-01 | 2004-03-15 | 115.1 |
| 788064_99999 | HOWARD AFB | 8.917 | -79.6 | 16 | 1980-04-20 | 1980-08-09 | 115.1 |
| 999999_10718 | FORT KOBBE | 8.917 | -79.6 | 15.8 | 1954-12-13 | 1968-01-01 | 115.1 |
| 788065_99999 | CHIVA-CHIVA | 9.033 | -79.6 | 27 | 1983-05-16 | 1985-10-08 | 116.7 |
| PMW00010716 | BALBOA | 8.933 | -79.57 | 13.1 | 1953-01-01 | 1957-12-31 | 118.9 |
| PMW00010701 | BALBOA ALBROOK | 8.967 | -79.55 | 9.1 | 1949-01-01 | 1967-12-31 | 121.1 |
| 783842_99999 | MARCOS A GELABERT I | 8.967 | -79.55 | 10 | 2005-02-08 | 2023-09-08 | 121.1 |
| 788071_99999 | ALBROOK AFB/BALBOA | 8.967 | -79.55 | 9 | 1936-10-01 | 1971-03-14 | 121.1 |
| 999999_10701 | BALBOA ALBROOK | 8.967 | -79.55 | 9.1 | 1949-01-01 | 1967-12-30 | 121.1 |
| 788067_99999 | MARCOS A. GELABERT | 8.983 | -79.52 | 13 | 1973-01-02 | 2001-09-21 | 124.9 |
| 749029_99999 | CALZADA LARGA | 9.166 | -79.55 | 120.1 | 1942-07-08 | 1944-11-01 | 125.9 |
| PM000128005 | LA MESA DE MACARACAS | 7.63 | -80.61 | 180 | 1979-01-01 | 1993-12-31 | 136.1 |
| 749034_99999 | POCRI AIRFIELD | 7.733 | -80.13 | 14 | 1943-06-15 | 1947-12-23 | 136.7 |
| 787920_99999 | TOCUMEN INTL | 9.071 | -79.38 | 41.1 | 1973-01-01 | 2023-09-08 | 140.9 |
| 788075_99999 | TOCUMEN/GEN. OMAR | 9.083 | -79.37 | 41 | 1973-06-01 | 1973-06-30 | 142.8 |
| 749028_99999 | LA JOYA | 9.133 | -79.23 | 48 | 1942-07-16 | 1945-09-27 | 158.3 |
#used for MERRA2 - small area coordinates +/- 2.5 degrees
xy <-ggmap::make_bbox(lon = c(0.5*round(Longitude/0.5,1)-2.5,0.5*round(Longitude/0.5,1)+2.5),
lat = c(0.5*round(Latitude/0.5,1)-2.5,0.5*round(Latitude/0.5,1)+2.5),f = 0)
#used for ggsn graphical scale, same ranges but 8% smaller
xy_min <-ggmap::make_bbox(lon = c(xy[1],xy[3]),
lat = c(xy[2],xy[4]),f=-0.08)
#ggmap capture for goggle maps
map <- get_map(xy,scale=2,zoom=7,maptype="hybrid",source="google",
messaging=FALSE)NOAA.file.par<-here::here("data","NOAA","NOAA Info par.rdata")
if(file.exists(NOAA.file.par)){
load(NOAA.file.par)
}else{
par.review<-data.frame(var=c("PRCP"),
op=c("sum"),
names=c("TOTAL PRECIPITATION"),
stringsAsFactors = F)
apply(par.review,1,FUN=function(par){
# par<-par.review[2,]
#variable definition
par.var<<-par[1]
par.fun<<-eval(parse(text = par[2]))
par.title<<-par[3]
#Title (two spaces before \n )
pander::pandoc.header(par.title,level=2)
# files
db.name<-paste0("NOAA.",par.var,".db")
idx.name<-paste0("NOAA.",par.var,".idx")
syn.name<-paste0("NOAA.",par.var,".syn")
#Prepare of DB
NOAA.par.db<-NOAA.db %>%
dplyr::filter(variable%in%par.var) %>%
dplyr::select(-NAME,-variable) %>%
as.data.frame() %>%
split(x=.[c("DATE","value")],f=.$ID)
#Prepare index
NOAA.par.idx<-NOAA.idx[NOAA.idx$ID%in%names(NOAA.par.db),]
#Sort db
NOAA.par.db<-NOAA.par.db[NOAA.par.idx$ID]
#days with information
st.dwi<-sapply(NOAA.par.db,FUN=function(x){dwi(x)%>%sum})
#1.selection of the stations with information more than 'n' years
st.sel<-st.dwi>=min.ywi*365.25
#clean the DB based on dwi
NOAA.par.db<-NOAA.par.db[st.sel]
NOAA.par.idx<-NOAA.par.idx[st.sel,]
Meteo.analysis(db = NOAA.par.db,
FUN.month =par.fun,
Min.days.annual = 300,
Max.miss.days.month = 5,output = syn.name,EC = F,
start.year = 1900) %>% invisible()
#Mean annual based on a range
NOAA.par.idx %<>%mutate(ANN=eval(parse(text=syn.name))$"Daily.Comp" %>%
window(start=as.Date("1980-01-01"),
end=as.Date("2019-12-31")) %>%
daily2annual.avg(FUN = par.fun) %>%
.$Ann)
#2.Remove is ANN is not available
st.sel<-!is.na(NOAA.par.idx$ANN)
#clean the DB based on dwi
NOAA.par.db<-NOAA.par.db[st.sel]
NOAA.par.idx<-NOAA.par.idx[st.sel,]
Meteo.analysis(db = NOAA.par.db,
FUN.month =par.fun,
Min.days.annual = 300,
Max.miss.days.month = 5,output = syn.name,EC = F,
start.year = 1900) %>% invisible()
#assign the variables
assign(x = db.name,value=NOAA.par.db,envir = .GlobalEnv)
assign(x = idx.name,value=NOAA.par.idx,envir = .GlobalEnv)
#remove the auxiliar variables
rm(par.var,par.fun,par.title,
envir = .GlobalEnv)
pander::pandoc.p('')
pander::pandoc.p('')
})
save(file = NOAA.file.par,
list =c(glue::glue('NOAA.{par.review$var}.db'),
glue::glue('NOAA.{par.review$var}.syn'),
glue::glue('NOAA.{par.review$var}.idx')))
}2.2 SRTM
Information obtained indirectly as a CGIAR-SRTM (90 m resolution), this topography covers an extension from 60 deg Lat North / South. For area beyond these limitations, it is recommend to use ArticDEM or another topographical source.
#Download SRTM
SRTM_ras<-SRTM_capture(bbox = xy,file = "SRTM",
trim = T,overwrite = T,
maxcell = 1e5)
# to a data.frame
SRTM_tbl<-raster::rasterToPoints(SRTM_ras) %>%
data.frame()
names(SRTM_tbl)<-c("x","y","value")
SRTM_tbl%<>%
dplyr::filter(complete.cases(value))
plot(SRTM_ras)
points(Longitude,Latitude)
2.3 MERRA2
MERRA2 (Modern-Era Retrospective analysis from Research and Applications, Version 2) from the National Aeronautics and Space Administration (NASA). MERRA includes daily data from 1981 to present (2023) for the entire world, based on a 0.5 degree latitude by 0.5 degree longitude grid, these records were obtained daily.
This analysis upscale the topographical information to match MERRA2 grid size (~50 x 50 km); also total precipitation is divided in rainfall and snowfall as a daily timeseries for every MERRA2 grid. This support some understanding of local / regional rainfall and snowfall.
It is recommended to use MERRA2 as a general overview; however, based on the analysis, these records can be complemented with. ERA5, ERA5-Land, PERSIANN, PERSIANN-CDR, PERSIANN-CCS, CHIRPS, and other local reanalysis sources depending on the geographical location.
merra2_file<-here::here("data","rds","MERRA2.par")
if(file.exists(merra2_file)){
merra2_db <- arrow::read_parquet(file = merra2_file)
}else{
#Multi core capabilities as a way to speed up downloading processes
cl<-parallel::makeCluster(parallel::detectCores())
clusterExport(cl,varlist = c("xy"))
merra2_db<-pbapply::pblapply(cl=cl,1981:year(now()),FUN=function(date_per){
library(magrittr)
merra2_tbl <- nasapower::get_power(
community = "ag",
lonlat = as.vector(xy),
pars = c("PRECTOTCORR","T2M"
,"T2M_MAX","T2M_MIN",
"RH2M","T2MDEW","WS2M",
"ALLSKY_SFC_SW_DWN"
),
dates = c(as.Date(paste0(date_per,"-01-01")), #from january first
#to the min between now and end of year
min(as.Date(lubridate::now())-lubridate::days(5),
as.Date(paste0(date_per,"-12-31")))),
temporal_api = "daily")
merra2_db<-merra2_tbl[,-c(3,4,5,6)] %>%
dplyr::rename(dates="YYYYMMDD") %>%
reshape2::melt(id.vars=c("LAT","LON","dates")) %>%
data.table::as.data.table()
}) %>% data.table::rbindlist()
stopCluster(cl)
arrow::write_parquet(x = merra2_db,
sink = merra2_file)
}
#Mean Annual values region
merra2_MA_tbl<-merra2_db %>%
dplyr::filter(variable%in%c("PRECTOTCORR","T2M","T2M_MAX","T2M_MIN","WS2M")) %>%
dplyr::mutate(years=lubridate::year(dates)) %>%
data.table::as.data.table() %>%
# dplyr::rename(PRCP="PRECTOTCORR") %>%
reshape2::dcast(data = .,formula = LON+LAT+dates+years~variable,
value.var = "value") %>%
mutate(RAIN=if_else(T2M>0,PRECTOTCORR,0), #rainfall, snowfall at daily scale
SNOW=if_else(T2M<=0,PRECTOTCORR,0)) %>%
data.table::as.data.table() %>%
dplyr::group_by(LON,LAT,years) %>%
dplyr::summarize(prcp=sum(PRECTOTCORR,na.rm=T),
rain=sum(RAIN,na.rm=T),
snow=sum(SNOW,na.rm=T),
tavg=mean(T2M,na.rm=T),
tmin=mean(T2M_MIN,na.rm=T),
tmax=mean(T2M_MAX,na.rm=T),
ws=mean(WS2M,na.rm=T)) %>%
dplyr::group_by(LON,LAT) %>%
dplyr::summarize(prcp=mean(prcp),
rain=mean(rain),
snow=mean(snow),
tavg=mean(tavg),
tmin=mean(tmin),
tmax=mean(tmax),
ws=mean(ws)) %>%
as.data.frame()
#WGS84
merra2_ras<-raster::rasterFromXYZ(merra2_MA_tbl)
crs(merra2_ras)<-crs("+init=epsg:4326")
merra2_ras$z<-resample(SRTM_ras,merra2_ras,method="bilinear")
#Merra table is overwritten to include elevation
merra2_MA_tbl<-rasterToPoints(merra2_ras) %>%
data.frame()
plot(merra2_ras)
rasterToPoints(merra2_ras) %>%
data.frame() %>%
mutate(coast=coast.dist(x,y)) %>%
dplyr::select(x,y,z,coast,tavg,tmin,tmax,prcp,rain,snow,ws) %>%
hydropairs(dec=2,
main="Correlation matrix for mean annual values based on MERRA2")
merra2_x<-merra2_db$LON %>% unique()
merra2_y<-merra2_db$LAT %>% unique()
#MERRA TS
sel_LON<-
merra2_x[(Longitude-merra2_x)^2%in%min((Longitude-merra2_x)^2)][1]
sel_LAT<-
merra2_y[(Latitude-merra2_y)^2%in%min((Latitude-merra2_y)^2)][1]
merra2_z<-merra2_db %>%
dplyr::filter(LON==sel_LON, LAT==sel_LAT) %>%
dplyr::select(-LON,-LAT) %>%
as.data.frame() %>%
#If the values is < - 50 then it is a NAN
# mutate(value=ifelse(value< -50,NaN,value)) %>%
reshape2::dcast(data=.,formula=dates~variable,value.var = "value",
fun.aggregate = mean) %>%
rename(prcp=PRECTOTCORR,
tavg=T2M,
tmin=T2M_MIN,
tmax=T2M_MAX) %>%
mutate(rain=if_else(tavg>0,prcp,0), #rainfall, snowfall at daily scale
snow=if_else(tavg<=0,prcp,0)) %>%
tk_zoo(silent=T)
plot.zoo(merra2_z,yax.flip = T,
main="Local time series based on MERRA2")
2.4 DAYMET (Just available in Mexico, US and Canada)
# library(daymetr)
#
# daymet_info<-daymetr::download_daymet(site = "Daymet",lat = Latitude,lon = Longitude,
# start = 2000,end=lubridate::year(now())-1,silent = T)
#
#
# daymet_z<-daymet_info$data %>%
# dplyr::mutate(dates=as.Date(paste0(year,"-01-01"))+lubridate::days(yday)-1) %>%
# dplyr::select(dates,prcp="prcp..mm.day.",tmax="tmax..deg.c.",tmin="tmin..deg.c.") %>%
# dplyr::mutate(tavg=0.5*tmax+0.5*tmin,
# rain=if_else(tavg>0,prcp,0),
# snow=if_else(tavg<=0,prcp,0)) %>%
# tk_zoo(silent=T)
#
# plot(daymet_z, yax.flip = T,
# main="Local time series based on Daymet")3 Total Precipitation
NOAA.PRCP.idx$NAME<-NOAA.PRCP.idx$NAME %>%make.names(unique=T) %>%
str_replace_all("\\."," ") %>%
str_wrap(string = .,width = 20)
st.names<-glue('{NOAA.PRCP.idx$NAME}\n{round(NOAA.PRCP.idx$ANN,1)} mm/yr')
dwi.graph(x=NOAA.PRCP.syn$"Daily.Comp",
station.name = st.names,
start.year = 1900)
NOAA_prcp_gp<-NOAA.PRCP.idx %>%
dplyr::select(x=LONGITUDE,y=LATITUDE,z=ELEVATION,ANN,NAME) %>%
reshape2::melt(id.vars=c("NAME","ANN"))
rbind(
merra2_MA_tbl %>%
dplyr::select(x,y,z,ANN=prcp) %>%
mutate(source="MERRA"),
NOAA.PRCP.idx %>%
dplyr::select(x=LONGITUDE,y=LATITUDE,z=ELEVATION,ANN) %>%
mutate(source="NOAA")
) %>%
reshape2::melt(id.vars=c("ANN","source")) %>%
ggplot(data=.,aes(x=value,y=ANN))+
geom_point(size=4,aes(col=source),alpha=0.4)+
geom_smooth(method="lm",se=F,show.legend = F,linetype=2,aes(col=source))+
geom_text_repel(data=NOAA_prcp_gp,aes(x=value,label=NAME),size=2.5)+
stat_poly_eq(aes(label = ifelse(after_stat(adj.r.squared) > 0.6,
paste(after_stat(adj.rr.label),
after_stat(eq.label),
sep = "*\", \"*"),
after_stat(adj.rr.label)),col=source),
formula = y ~ poly(x, 1, raw = TRUE),label.x.npc = 0.1,vstep = 0.1,
size=4)+
geom_vline(data=site.info_xy,aes(xintercept = value),linetype=2)+
theme_bw()+
labs(col="source:",x="",title="Mean annual comparison values MERRA and NOAA",
subtitle="x as Longitude, y as Latitude and z as Elevation",
caption="Note: site as dotted line")+
theme(legend.position = "bottom")+
facet_wrap(~variable,scales = "free_x",ncol=1)
NOAA.PRCP.idx %>%
dplyr::select(LONGITUDE,LATITUDE,ELEVATION,ANN) %>%
data.frame %>%
psych::pairs.panels(smooth = T,scale = T,cor = T,
smoother = F,ci=F,cex.cor = 0.8,lm=T,
main="Correlation matrix for parameters - NOAA")
prcp.sel<-1:4
prcp.comp<-merge(NOAA.PRCP.syn$Daily.Comp[,prcp.sel],
MERRA2=merra2_z$prcp)
names(prcp.comp)[prcp.sel]<-NOAA.PRCP.idx$NAME[prcp.sel]
# prcp.comp %>%
# data.frame %>%
# psych::pairs.panels(smooth = T,scale = T,cor = T,
# smoother = T,ci=F,cex.cor = 0.8,lm=T,
# main="Correlation matrix for daily precipitation")
prcp.comp %>%
daily2monthly(FUN=sum) %>%
data.frame %>%
psych::pairs.panels(smooth = T,scale = F,cor = T,
smoother = F,ci=F,cex.cor = 0.8,lm=T,
main="Correlation matrix for monthly precipitation")
#Definition of the ranges
prcp_ranges<-c(merra2_MA_tbl$prcp,NOAA.PRCP.idx$ANN) %>% na.omit()
prcp_bin<-seq(min(prcp_ranges) %>% floor,
max(prcp_ranges) %>% ceiling,length.out=11) %>%
signif(2) %>% unique()
prcp_lbl<-paste(prcp_bin[1:(length(prcp_bin)-1)] %>%
str_pad(string = .,width = 4,side = "left",pad = "0"),
"–",
prcp_bin[2:(length(prcp_bin))] %>%
str_pad(string = .,width = 4,side = "left",pad = "0"))
#Map
ggmap(map, darken = c(0.3, "white"))+
#TOPOGRAHY
# geom_tile(data = SRTM_tbl, aes(x = x, y = y,
# fill=cut(value,seq(0,8000,500),
# include.lowest=T,dig.lab=10)),
# alpha=0.4)+
#MERRA
geom_tile(data = merra2_MA_tbl,aes(x=x,y=y,
fill=cut(prcp,prcp_bin,prcp_lbl,dig.lab=5,ordered_result=T)),alpha=0.6)+
geom_label_repel(data=NOAA.PRCP.idx,
aes(x=LONGITUDE,y=LATITUDE,
label=paste0(NAME,"\n",round(ANN,1))),colour="darkblue",
fontface="bold",alpha=0.6,
box.padding =unit(0.5, "lines"),segment.color = "black",
size=3 )+
geom_point(data = NOAA.PRCP.idx,
aes(x=LONGITUDE,y=LATITUDE,
fill=cut(ANN,prcp_bin,prcp_lbl,dig.lab=5,ordered_result=T)),size=4,shape=22,alpha=0.6)+
geom_point(data=site.info,
aes(x=Longitude,y=Latitude),shape=5,col="red")+
geom_label_repel(data=site.info,
aes(x=Longitude,y=Latitude,label=Location))+
scale_fill_brewer(palette="Spectral",direction = -1)+
labs(size="Elevation [masl]",x="Longitude [deg]",y="Latitude [deg]",
shape="Station Name",
fill = "Mean Annual\nPrecipitation\nBackground from MERRA2\nPoints from NOAA",
# fill = "Topography\nfrom SRTM",
# col="MAAT from\nNOAA[degC]\n(1980-2022)",
caption=
"Note: Site is presented as a red diamond")+
theme(legend.position="bottom")+
ggsn::scalebar(x.min=xy_min[1],x.max=xy_min[3],
y.min=xy_min[2],y.max=xy_min[4],dist = 50,dist_unit = "km",transform = T,model = "WGS84")+
coord_quickmap()
4 Site review
prcp_2021 <- read_excel(here::here("data","csv","Precipitation 2015 to 2023.xls"),
sheet = "2021", skip = 2)
prcp_2022 <- read_excel(here::here("data","csv","Precipitation 2015 to 2023.xls"),
sheet = "2022", skip = 2)
prcp_2023 <- read_excel(here::here("data","csv","Precipitation 2015 to 2023.xls"),
sheet = "2023",col_names = F)
names(prcp_2023)<-c("Date","mm")
prcp_z<-rbind(
prcp_2021[,1:2],
prcp_2022[,1:2],
prcp_2023[,1:2]) %>%
tk_zoo()
time(prcp_z)<-as.Date(time(prcp_z))
site_max<-IDF(prcp_z,1) %>% coredata()
reg_sta_max<-NOAA.PRCP.syn$Daily.Comp$PM000105001 %>%
IDF(.,1)
reg_sta_max<-coredata(reg_sta_max)[!is.infinite(coredata(reg_sta_max))]data.frame(prob=FA.val,
RP=RP.text(FA.val,dry.wet=F),
site=FA(1.13*site_max,distrib = "GUMBEL")$results[,2],
"boca del toabre"=FA(1.13*reg_sta_max,distrib = "GUMBEL")$results[,2]) %>%
dplyr::filter(prob>=0.5,prob<=0.99) %>%
reshape2::melt(id.vars=c("prob","RP")) %>%
ggplot(data=.,aes(x=RP,y=value,col=variable))+
geom_point()+
geom_line()+
scale_x_log10(limits=c(NA,100),breaks=c(2,5,10,25,50,100))+
scale_y_log10(breaks=c(150,200,300,400),limits=c(100,NA))+
annotation_logticks(sides = "bl")+
labs(x="return period [yrs]",y="maximum precipitation 24 hrs [mm]",
col="source:",title="Maximum precipitation based on Gumbel Distribution",
caption="note: maximum values corrected with 1.13 factor from daily to 24 hours")+
theme_light()+
theme(legend.position="bottom")
5 Synthesis
Information at the site based on MERRA2. As a surface water analysis, you have to consider this work as a simple overview; which assumes MERRA2 is representative of this location.
This analysis does not include or considers snow capture efficiency; which it should change or modified local snowfall estimations.
BTW## Climograph
source_gp<-"MERRA2"
reanalysis_z<-merra2_z
#In the case of DAYMET "unlock" these lines
# source_gp<-"DAYMET"
#
# reanalysis_z<-daymet_z
year_start_z<-year(start(reanalysis_z))
year_end_z<-year(last(reanalysis_z))
#Mean monthly
reanalysis_mean_mon<-rbind(
daily2annual.avg(reanalysis_z,FUN = sum) %>%
dplyr::filter(Station%in%c("prcp","rain","snow")),
daily2annual.avg(reanalysis_z,FUN = mean)%>%
dplyr::filter(!Station%in%c("prcp","rain","snow"))
)
reanalysis_mean_mon_gp<-reanalysis_mean_mon %>%
reshape2::melt(data = .,id.vars="Station") %>%
dplyr::filter(!variable=="Ann") %>%
reshape2::dcast(data =.,formula=variable~Station)
prcp_range_max<-max(reanalysis_mean_mon_gp$prcp)
prcp_range_min<-min(reanalysis_mean_mon_gp$prcp,0)
tavg_range_max<-max(reanalysis_mean_mon_gp$tavg)
tavg_range_min<-min(reanalysis_mean_mon_gp$tavg)
#ranges in Kelvin
tavg_range_max<-(tavg_range_max+273.15)*1.003-273.15
tavg_range_min<-(tavg_range_min+273.15)*0.998-273.15
slope_temp<-(prcp_range_max-prcp_range_min)/(tavg_range_max-tavg_range_min)
inter_temp<-prcp_range_min-tavg_range_min*slope_temp
reanalysis_mean_mon_prcp<-reanalysis_mean_mon_gp %>%
rename(mon="variable") %>%
reshape::melt.data.frame(measure.vars=c("rain","snow"))
#Climograph
ggplot(data=reanalysis_mean_mon_gp,aes(x=variable))+
geom_bar(data=reanalysis_mean_mon_prcp,stat = "Identity", aes(x=mon,y=value,fill=variable) ,col="black") +
geom_point( aes(y=slope_temp*tavg+inter_temp), size=2, color="red",alpha=0.6,shape=21,fill="red")+
geom_line(aes(y=slope_temp*tavg+inter_temp,group="tavg"), size=1, color="red",alpha=0.6)+
scale_y_continuous(sec.axis=sec_axis(trans=~(.-inter_temp)/slope_temp,
name="Temperature [degC]"))+
geom_text(aes(y=prcp*1.02,label=round(prcp,1)),alpha=0.6,nudge_y = 0.6)+
geom_text(aes(y=slope_temp*tavg+inter_temp,label=round(tavg,1)),alpha=0.6,
nudge_y = -0.6,col="red",fontface="bold")+
theme_bw()+
theme(axis.line.y.right = element_line(color = "red"),
axis.ticks.y.right = element_line(color = "red"),
axis.text.y.right = element_text(color = "red"),
axis.title.y.right = element_text(color = "red")
) +
labs(x="",y="Total Precipitation [mm]",
title=glue::glue("Climatogram for Longitude: {round(Longitude,2)} / Latitude: {round(Latitude,2)}"),
subtitle=glue::glue("Based on {source_gp} ({year_start_z} to {year_end_z})"),
fill="")+
scale_fill_manual(values = c("rain"="dodgerblue3","snow"="snow2"))+
theme(legend.position="bottom")
pander::pandoc.table(reanalysis_mean_mon,split.table=Inf,style="simple",round=1,
caption=glue::glue("Mean monthly parameters based on {source_gp}"))| Station | Jan | Feb | Mar | Apr | May | Jun | Jul | Aug | Sep | Oct | Nov | Dec | Ann |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| prcp | 25.1 | 12.6 | 18.3 | 87.8 | 219.1 | 241.6 | 230.8 | 238.7 | 246.3 | 247.5 | 234.5 | 87.9 | 1890 |
| rain | 25.1 | 12.6 | 18.3 | 87.8 | 219.1 | 241.6 | 230.8 | 238.7 | 246.3 | 247.5 | 234.5 | 87.9 | 1890 |
| snow | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| tavg | 24.7 | 25.2 | 25.9 | 26.5 | 26.3 | 25.9 | 25.6 | 25.6 | 25.6 | 25.4 | 25.2 | 25 | 25.6 |
| tmax | 27.6 | 28.7 | 30 | 30.5 | 29.3 | 28.2 | 27.7 | 27.8 | 27.9 | 27.7 | 27.6 | 27.6 | 28.4 |
| tmin | 22.8 | 22.9 | 23.3 | 24 | 24.4 | 24.3 | 24.1 | 24 | 23.9 | 23.7 | 23.6 | 23.2 | 23.7 |
| RH2M | 83.4 | 78.7 | 75.2 | 76.8 | 84 | 87.5 | 88.9 | 88.7 | 87.5 | 87.9 | 88.6 | 86.8 | 84.5 |
| T2MDEW | 21.6 | 21 | 20.9 | 21.8 | 23.3 | 23.6 | 23.6 | 23.5 | 23.3 | 23.2 | 23.1 | 22.6 | 22.6 |
| WS2M | 2.5 | 2.6 | 2.6 | 2.3 | 1.6 | 1.4 | 1.5 | 1.4 | 1.3 | 1.4 | 1.7 | 2.2 | 1.9 |
| ALLSKY_SFC_SW_DWN | 17.7 | 19.6 | 20.6 | 19.3 | 16.4 | 15 | 15 | 15.3 | 15.7 | 15 | 14.2 | 15.5 | 16.6 |
5.1 Koppen- Geiger Climate Classification
clim_class_tbl<-data.frame(col=
c("#960000", "#FF0000", "#FF6E6E", "#FFCCCC", "#CC8D14", "#CCAA54", "#FFCC00", "#FFFF64", "#007800", "#005000", "#003200", "#96FF00", "#00D700", "#00AA00", "#BEBE00", "#8C8C00", "#5A5A00", "#550055", "#820082", "#C800C8", "#FF6EFF", "#646464", "#8C8C8C", "#BEBEBE", "#E6E6E6", "#6E28B4", "#B464FA", "#C89BFA", "#C8C8FF", "#6496FF", "#64FFFF", "#F5FFFF"),
class=c('Af', 'Am', 'As', 'Aw', 'BSh', 'BSk', 'BWh', 'BWk', 'Cfa', 'Cfb','Cfc', 'Csa', 'Csb', 'Csc', 'Cwa','Cwb', 'Cwc', 'Dfa', 'Dfb', 'Dfc','Dfd', 'Dsa', 'Dsb', 'Dsc', 'Dsd','Dwa', 'Dwb', 'Dwc', 'Dwd', 'EF','ET', 'Ocean'),
description=c("Tropical rainforest climate", "Tropical monsoon climate", "Tropical dry savanna climate", "Tropical savanna, wet", "Hot semi-arid (steppe) climate", "Cold semi-arid (steppe) climate", "Hot deserts climate", "Cold desert climate", "Humid subtropical climate", "Temperate oceanic climate", "Subpolar oceanic climate", "Hot-summer Mediterranean climate", "Warm-summer Mediterranean climate", "Cool-summer Mediterranean climate", "Monsoon-influenced humid subtropical climate", "Subtropical highland climate or temperate oceanic climate with dry winters", "Cold subtropical highland climate or subpolar oceanic climate with dry winters", "Hot-summer humid continental climate", "Warm-summer humid continental climate", "Subarctic climate", "Extremely cold subarctic climate", "Hot, dry-summer continental climate", "Warm, dry-summer continental climate", "Dry-summer subarctic climate", "Very Cold Subarctic Climate", "Monsoon-influenced hot-summer humid continental climate", "Monsoon-influenced warm-summer humid continental climate", "Monsoon-influenced subarctic climate", "Monsoon-influenced extremely cold subarctic climate", "Ice cap climate", "Tundra", "Ocean"))
clim_class<-data.table(lon=kgc::genCoords(latlong='lon',full='true'),
lat=kgc::genCoords(latlong='lat',full='true'),
kg=kgc::kmz) %>%
dplyr::filter(lon>=xy[1],lon<=xy[3],lat>=xy[2],lat<=xy[4]) %>%
as.data.frame() %>%
mutate(kg_clas=getZone(kg)) %>%
left_join(.,clim_class_tbl,by=c("kg_clas"="class")) %>%
mutate(description=paste0(kg_clas,": ",description) %>%
str_wrap(.,30))
clim_class_fill<-clim_class %>%
dplyr::select(kg_clas,col,description) %>% unique() %>%
arrange(kg_clas)
ggmap(map, darken = c(0.3, "white"))+
geom_tile(data=clim_class,aes(x=lon,y=lat,fill=description),alpha=0.4)+
geom_point(data=site.info,
aes(x=Longitude,y=Latitude),shape=5,col="red")+
geom_label_repel(data=site.info,
aes(x=Longitude,y=Latitude,label=Location))+
geom_label_repel(data=NOAA.PRCP.idx,
aes(x=LONGITUDE,y=LATITUDE,
label=paste0(NAME)),colour="darkblue",
fontface="bold",alpha=0.6,
box.padding =unit(0.5, "lines"),segment.color = "black",
size=3 )+
geom_point(data = NOAA.PRCP.idx,
aes(x=LONGITUDE,y=LATITUDE),size=4,shape=22,alpha=0.6)+
coord_quickmap()+
ggsn::scalebar(x.min=xy_min[1],x.max=xy_min[3],
y.min=xy_min[2],y.max=xy_min[4],dist = 50,dist_unit = "km",transform = T,model = "WGS84")+
scale_fill_manual(breaks=clim_class_fill$description,values = clim_class_fill$col)+
labs(caption="Rubel, F., Brugger, K., Haslinger, K., Auer, I., 2016. The climate of the European Alps: Shift
of very high resolution Köppen-Geiger climate zones 1800–2100. Meteorologische Zeitschrift.
doi:10.1127/metz/2016/0816 http://koeppen-geiger.vu-wien.ac.at/",
title="Koppen-Geiger Climate Classification",
fill="",
x="Longitude [deg]",y="Latitude[deg]")+
theme(legend.position="bottom")+
guides(fill=guide_legend(ncol=3))