Meteorological Review
Victor Munoz S.
2023-08-16
1 Introduction
Longitude = -49.39485
Latitude = -14.23154
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 | -49.39 | -14.23 | 384.6 |
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,]
# Output
NOAA.capture.GHCND()
NOAA.capture.GSOD()
NOAA.prepare()
save(file = NOAA.file,list = c("NOAA.db","NOAA.idx"))
}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","TAVG","TMAX","TMIN"),
op=c("sum","mean","mean","mean"),
names=c("TOTAL PRECIPITATION",
"MEAN AIR TEMPERATURE",
"MAXIMUM AIR TEMPERATURE",
"MINIMUM AIR TEMPERATURE"),
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')))
}
# NOAA.PRCP.idx %>%
# dplyr::select(-WMO_ID,-ICAO,-MINDATE,-MAXDATE) %>%
#
# pandoc.table(split.table=Inf)3 Create NETCDF files
As an overall review 3 reanalysis sources were presented as comparison: MERGE (local reanalysis), CHIRPS and ERA5-LAND. ERA5-LANd is typically considered an upgrade from ERA5 (used in Golder report); however, normally ERA5-LAND tends to perform similarly as ERA5.
3.1 MERGE
library(terra)
# merge_tbl<-readr::read_rds(here::here("data","rds","MERGE_db.rds"))
#
# merge_lst<-split(merge_tbl,merge_tbl$yyyymmdd)
#
# merge_ras<-pbapply::pblapply(merge_lst,FUN=function(x){
#
# rast_1<-terra::rast(x[,1:3],type="xyz")
#
# terra::time(rast_1)<-dplyr::pull(x[1,4])
#
# return(rast_1)
#
# }) %>% terra::rast()
#
# terra::writeCDF(x = merge_ras,filename = here::here("data","netcdf","merra_dly.nc"))
#
# mean(tapp(merge_ras,index="years",fun=sum)) %>% plot()
#
# points(Longitude,Latitude,add=T)
reanalysis_fn(file = here::here("data","netcdf","MERGE_dly_prcp.nc"),
par = "prcp",fun=sum,yrs=c(1980:2020))
points(Longitude,Latitude,add=T)
3.2 CHIRPS
#
# chirps_tbl<-readr::read_rds(here::here("data","rds","chirps_db.rds"))
#
# chirps_lst<-split(chirps_tbl,chirps_tbl$yyyymmdd)
#
# chirps_ras<-pbapply::pblapply(chirps_lst,FUN=function(x){
#
# rast_1<-terra::rast(x[,c(1,2,4)],type="xyz")
#
# terra::time(rast_1)<-dplyr::pull(x[1,3])
#
# return(rast_1)
#
# }) %>% terra::rast()
#
# terra::writeCDF(x = chirps_ras,filename = here::here("data","netcdf","chirps_dly.nc"))
#
# mean(tapp(chirps_ras,index="years",fun=sum)) %>% plot()
#
# points(Longitude,Latitude,add=T)
reanalysis_fn(file = here::here("data","netcdf","CHIRPS_dly_prcp.nc"),
par = "prcp",fun=sum,yrs=c(1980:2020))
points(Longitude,Latitude,add=T)
3.3 ERA5 - LAND
#
# era5l_files<-dir(here::here("data","netcdf","ERA5-Land"),pattern = ".nc",full.names = T)
#
# era5l_ras<-pbapply::pblapply(era5l_files,FUN=function(x){
#
# # x<-era5l_files[[1]]
#
# rast_1<-terra::rast(x)
#
# # plot(rast_1[[c(1,2,4,8,20,22,23,24,25,26,27)]]*1000)
#
# time_sel<-lubridate::hour(terra::time(rast_1))==0
#
# # lubridate::hour(terra::time(rast_1))[c(1,2,4,8,20,22,23,24,25,26,27)]
#
# rast_2<-rast_1[[time_sel]]*1000
#
# terra::time(rast_2)<-as.Date(terra::time(rast_2))-days(1)
#
# return(rast_2)
#
# }) %>% terra::rast()
#
# terra::writeCDF(x = era5l_ras,filename = here::here("data","netcdf","ERA5-LAND_dly_prcp.nc"),overwrite=TRUE)
# mean(tapp(era5l_ras,index="years",fun=sum)[[1:42]]) %>% plot()
reanalysis_fn(file = here::here("data","netcdf","ERA5-LAND_dly_prcp.nc"),
par = "prcp",fun=sum,yrs=c(1980:2020))
points(Longitude,Latitude,add=T)
4 NOAA Review
The regional information was compared with records from CHIRPS and ERA5-LAND.MERGE was not included because NOAA records tends to be before 1994; while MERGE is available from 2000.
lapply(1:5,FUN=function(x){
site_name<-NOAA.PRCP.idx$NAME[x]
noaa_lat<-NOAA.PRCP.idx$LATITUDE[x]
noaa_lon<-NOAA.PRCP.idx$LONGITUDE[x]
pander::pandoc.header(site_name,level=2)
prcp_comb_z<-merge(NOAA=NOAA.PRCP.syn$Daily.Comp[,x],
# MERGE=nc2z_fn(Longitude = noaa_lon,Latitude = noaa_lat,
# source = "MERGE"),
CHIRPS=nc2z_fn(Longitude = noaa_lon,Latitude = noaa_lat,
source = "CHIRPS"),
"ERA5-LAND"=nc2z_fn(Longitude = noaa_lon,Latitude = noaa_lat,
source = "ERA5-LAND") %>%
window(end=as.Date("2022-10-30")))
pander::pandoc.p('')
pander::pandoc.p('')
data.frame(prcp_comb_z %>% daily2monthly.V2(FUN=sum)) %>%
hydropairs(dec = 2)
pander::pandoc.p('')
pander::pandoc.p('')
daily2annual.avg(prcp_comb_z,FUN=sum) %>%
pandoc.table(split.table=Inf,caption=site_name)
pander::pandoc.p('')
pander::pandoc.p('')
prcp_comb_z%>%
fortify.zoo()%>%
reshape2::melt(id.vars=c("Index","NOAA")) %>%
mutate(mon=month(Index) %>% as.factor()) %>%
openair::TaylorDiagram(mydata = .,obs = "NOAA",
mod = "value",group=c("variable"),
main=paste0(site_name,"- Daily Records"),
key.title="source:")
pander::pandoc.p('')
pander::pandoc.p('')
prcp_comb_z%>%
daily2monthly.V2(FUN=sum) %>%
fortify.zoo()%>%
reshape2::melt(id.vars=c("Index","NOAA")) %>%
mutate(mon=month(Index) %>% as.factor()) %>%
openair::TaylorDiagram(mydata = .,obs = "NOAA",
mod = "value",group=c("variable"),
main=paste0(site_name,"- Monthly Records"),
key.title="source:")
pander::pandoc.p('')
pander::pandoc.p('')
}) %>% invisible()4.1 SANTA TEREZINHA DE GOIAS
| Station | Jan | Feb | Mar | Apr | May | Jun | Jul | Aug | Sep | Oct | Nov | Dec | Ann |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| NOAA | 350.6 | 251.1 | 206.6 | 83 | 27.3 | 4.295 | 3.414 | 8.8 | 45.04 | 150.7 | 204.5 | 304 | 1639 |
| CHIRPS | 305.8 | 247.5 | 220.5 | 103 | 17.03 | 3.558 | 1.436 | 4.147 | 35.6 | 128.5 | 220.2 | 290 | 1577 |
| ERA5-LAND | 283.8 | 213.1 | 204.1 | 90.44 | 13.57 | 2.205 | 1.564 | 3.438 | 24.61 | 104.8 | 208.2 | 282.6 | 1432 |
4.2 PORTO URUACU
| Station | Jan | Feb | Mar | Apr | May | Jun | Jul | Aug | Sep | Oct | Nov | Dec | Ann |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| NOAA | 280.3 | 213.1 | 207.1 | 88.5 | 24.23 | 4.177 | 1.858 | 8.05 | 41.64 | 130.2 | 199.3 | 288.5 | 1487 |
| CHIRPS | 273.1 | 245.1 | 218.4 | 89.92 | 18.84 | 3.526 | 1.439 | 4.263 | 31.97 | 131.2 | 223.1 | 283.1 | 1524 |
| ERA5-LAND | 282.1 | 208.3 | 213.7 | 92.51 | 16.67 | 2.998 | 1.484 | 4.516 | 26.54 | 121.9 | 233.2 | 295.8 | 1500 |
4.3 PILAR DE GOIAS
| Station | Jan | Feb | Mar | Apr | May | Jun | Jul | Aug | Sep | Oct | Nov | Dec | Ann |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| NOAA | 458.6 | 339.5 | 260.7 | 121.4 | 33.08 | 3.948 | 4.971 | 16.51 | 53.75 | 159.8 | 247.7 | 381.8 | 2082 |
| CHIRPS | 300.7 | 237.8 | 213.2 | 106.2 | 18.81 | 3.024 | 1.494 | 4.237 | 40.94 | 135.5 | 221.4 | 286.8 | 1570 |
| ERA5-LAND | 305.8 | 235.2 | 231.8 | 106.6 | 18 | 3.525 | 2.447 | 5.124 | 31.86 | 131.9 | 236.8 | 296.6 | 1606 |
4.4 ESTRELA DO NORTE
| Station | Jan | Feb | Mar | Apr | May | Jun | Jul | Aug | Sep | Oct | Nov | Dec | Ann |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| NOAA | 319.1 | 266.7 | 259.5 | 114.1 | 21.1 | 3.874 | 3.791 | 7.639 | 49.44 | 142.4 | 239.8 | 340.9 | 1768 |
| CHIRPS | 289.6 | 239.9 | 233.2 | 124 | 18.49 | 2.535 | 1.417 | 2.674 | 36.5 | 125.8 | 237.9 | 312.9 | 1625 |
| ERA5-LAND | 335.5 | 286.6 | 291.3 | 143 | 23.15 | 3.123 | 1.187 | 4.647 | 32.77 | 145.4 | 263.6 | 327.1 | 1857 |
4.5 CRIXAS
| Station | Jan | Feb | Mar | Apr | May | Jun | Jul | Aug | Sep | Oct | Nov | Dec | Ann |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| NOAA | 258.2 | 210.9 | 192.3 | 91.12 | 7.758 | 6.975 | 0.6364 | 12.24 | 37.5 | 117.4 | 180 | 286.1 | 1401 |
| CHIRPS | 326.6 | 229 | 210.3 | 107.5 | 18 | 3.024 | 1.313 | 4.17 | 36.91 | 132.8 | 224.4 | 292.4 | 1586 |
| ERA5-LAND | 299.3 | 230.3 | 219.4 | 97.62 | 16.54 | 2.897 | 2.139 | 4.244 | 28.59 | 122.6 | 231.4 | 309.7 | 1565 |
5 MAP Review
Mean annual precipitation is compared between climate reanalisis sources and NOAA records, the site is presented as a red dot. It is noted as ERA5-Land and CHIRPS tends to present “similar” regional trends.
Differences must be associated to the different periods as NOAA are available before 1994 and reanalysis sources tends to be presented from 1980 to 2020; which the exception of MERGE from 2000.
library(tidyterra)
NOAA.PRCP.idx<-NOAA.PRCP.idx %>%
dplyr::filter(ANN>1000&ANN<2500) %>%
dplyr::filter(DISTANCE<150)
NOAA.PRCP.syn$Daily.Comp<-NOAA.PRCP.syn$Daily.Comp[,NOAA.PRCP.idx$ID]
dwi.graph(NOAA.PRCP.syn$Daily.Comp,start.year = 1900,
station.name = NOAA.PRCP.idx$NAME)
| ID | NAME | LATITUDE | LONGITUDE | ELEVATION | DISTANCE | ANN |
|---|---|---|---|---|---|---|
| BR001449002 | SANTA TEREZINHA DE GOIAS | -14.43 | -49.71 | 400 | 35.75 | 1640 |
| BR001449001 | PORTO URUACU | -14.52 | -49.05 | 547 | 49.34 | 1580 |
| BR001449000 | PILAR DE GOIAS | -14.76 | -49.58 | 850 | 57.47 | 2127 |
| BR001349000 | ESTRELA DO NORTE | -13.87 | -49.07 | 0 | 58.42 | 1771 |
| BR001449003 | CRIXAS | -14.53 | -49.96 | 0 | 64.74 | 1401 |
| BR001448005 | PALMEIRINHA | -14.02 | -48.61 | 0 | 92.14 | 1271 |
| BR001349002 | PORANGATU (DESCOBERTO) | -13.45 | -49.14 | 600 | 95.92 | 1598 |
| BR001549004 | NOVA AMERICA | -15.02 | -49.89 | 800 | 97.74 | 1885 |
| BR001450002 | GOVERNADOR LEONINO | -14.1 | -50.33 | 0 | 99.41 | 1460 |
| BR001448001 | NIQUELANDIA | -14.48 | -48.46 | 0 | 106.8 | 1783 |
| BR001348003 | TROMBAS | -13.51 | -48.74 | 450 | 112.1 | 1556 |
| BR001448002 | PONTE QUEBRA LINHA | -14.98 | -48.68 | 0 | 113 | 1410 |
| BR001349003 | ENTRONCAMENTO SAO MIGUEL | -13.27 | -49.2 | 427 | 113.5 | 1693 |
| BR001349001 | NOVO PLANALTO | -13.24 | -49.51 | 250 | 114.6 | 1656 |
| BR001350001 | RIO PINTADO (FAZ.PONTAL) | -13.53 | -50.19 | 200 | 116.4 | 1339 |
| BR001549000 | CERES (POSTO BIQUINHA) | -15.31 | -49.55 | 550 | 117 | 1546 |
| BR001549001 | GOIANESIA | -15.33 | -49.12 | 0 | 122.9 | 1540 |
| BR001448003 | PORTO RIO BAGAGEM | -14.37 | -48.2 | 373 | 132.6 | 1400 |
| BR001450001 | MOZARLANDIA (CHAC.FOGUEIRA) | -14.74 | -50.58 | 400 | 135.2 | 1675 |
| BR001350002 | SAO MIGUEL DO ARAGUAIA | -13.27 | -50.16 | 331 | 136.4 | 1893 |
| BR001448004 | MOQUEM-FAZ.VAU DA ONCA | -14.55 | -48.17 | 0 | 139 | 1597 |
| BR001549009 | URUANA | -15.5 | -49.69 | 0 | 140.3 | 1624 |
| BR001450000 | LAGOA DA FLECHA (FAZENDA) | -14.33 | -50.73 | 200 | 140.9 | 1505 |
| BR001448000 | COLINAS DO SUL | -14.15 | -48.08 | 0 | 145.7 | 1526 |
ggplot()+
geom_spatraster(data = era5_land_ma_ras, alpha = 0.6,na.rm = T)+
geom_point(data=NOAA.PRCP.idx,aes(x=LONGITUDE,y=LATITUDE,fill=ANN),size=5,shape=22)+
geom_text_repel(data=NOAA.PRCP.idx,aes(x=LONGITUDE,y=LATITUDE,label=NAME),size=2.5)+
geom_point(data=site.info,aes(x = Longitude, y = Latitude), color = "black",
shape = 21, fill = "red",size=5)+
scale_fill_binned(type = "viridis",breaks=seq(1000,2100,200))+
labs(title="ERA5-Land and NOAA")
ggplot()+
geom_spatraster(data = merge_ma_ras, alpha = 0.6,na.rm = T)+
geom_point(data=NOAA.PRCP.idx,aes(x=LONGITUDE,y=LATITUDE,fill=ANN),size=5,shape=22)+
geom_text_repel(data=NOAA.PRCP.idx,aes(x=LONGITUDE,y=LATITUDE,label=NAME),size=2.5)+
geom_point(data=site.info,aes(x = Longitude, y = Latitude), color = "black",
shape = 21, fill = "red",size=5)+
scale_fill_binned(type = "viridis",breaks=seq(1000,2100,200))+
labs(title="MERGE and NOAA")
ggplot()+
geom_spatraster(data = chirps_ma_ras, alpha = 0.6,na.rm = T)+
geom_point(data=NOAA.PRCP.idx,aes(x=LONGITUDE,y=LATITUDE,fill=ANN),size=5,shape=22)+
geom_text_repel(data=NOAA.PRCP.idx,aes(x=LONGITUDE,y=LATITUDE,label=NAME),size=2.5)+
geom_point(data=site.info,aes(x = Longitude, y = Latitude), color = "black",
shape = 21, fill = "red",size=5)+
scale_fill_binned(type = "viridis",breaks=seq(1000,2100,200))+
labs(title="CHIRPS and NOAA")
6 PMP Review
topo_ras<-terra::rast(here::here("data","dem","SRTM.bil"))
NOAA.PRCP.idx$PMP<-sapply(1:nrow(NOAA.PRCP.idx),FUN=function(x){
noaa_prcp_z<-NOAA.PRCP.syn$Daily.Comp[,x]
year_sel<-dwi(noaa_prcp_z)>200
year_max<-IDF(noaa_prcp_z,1)
year_max_val<-coredata(year_max)[year_sel]
PMP_val<-PMP(val.max = year_max_val,input.hrs = 24,number.of.obs.unit = 1,area.km = 1,graphs = F)
return(PMP_val$PMP)
})
NOAA.PRCP.idx %>%
mutate(z=mapply(GM_elev,LONGITUDE,LATITUDE,key_google)) %>%
dplyr::select(LATITUDE,LONGITUDE,ANN,z,PMP) %>%
data.frame() %>%
hydropairs()
ggplot()+
geom_spatraster(data = topo_ras, alpha = 0.6,na.rm = T)+
geom_point(data=NOAA.PRCP.idx,aes(x=LONGITUDE,y=LATITUDE,col=PMP),size=5)+
geom_text_repel(data=NOAA.PRCP.idx,aes(x=LONGITUDE,y=LATITUDE,label=round(PMP,0)),size=2.5)+
geom_point(data=site.info,aes(x = Longitude, y = Latitude), color = "black",
shape = 21, fill = "white",size=5)+
# scale_fill_continuous()+
scale_color_viridis_c(option = "A")+
scale_fill_binned(type = "viridis",breaks=seq(0,1600,200))+
labs(title="Topography and PMP")
7 Max 1day precipitation Review
Golder’s precipitation seems in the higher range compared with NOAA records.
No relationship between maximum precipitation with the exception of PMP and small correlations with MAP.
prcp_max_tbl<-lapply(1:nrow(NOAA.PRCP.idx),FUN=function(x){
noaa_prcp_z<-NOAA.PRCP.syn$Daily.Comp[,x]
year_sel<-dwi(noaa_prcp_z)>200
year_max<-IDF(noaa_prcp_z,1)
year_max_val<-coredata(year_max)[year_sel]
fa_res<-FA(year_max_val,distrib = "GUMBEL") %>% .$results %>%
mutate(RP=RP.text(Prob))
names(fa_res)[2]<-"prcp_max"
fa_res %>%
mutate(site=NOAA.PRCP.idx$NAME[x]) %>%
reshape2::dcast(data=.,formula=site~Prob,value.var = "prcp_max")
}) %>% do.call(rbind,.)
NOAA.PRCP.idx$"prcp_100_yr"<-prcp_max_tbl$"0.99"
NOAA.PRCP.idx<-NOAA.PRCP.idx %>%
mutate(z=mapply(GM_elev,LONGITUDE,LATITUDE,key_google))
prcp_max_tbl %>%
reshape2::melt() %>%
mutate(variable=as.numeric(as.character(variable))) %>%
mutate(RP=RP.text(variable,FALSE)) %>%
dplyr::filter(variable>=0.5) %>%
left_join(.,NOAA.PRCP.idx,by=c("site"="NAME")) %>%
select(variable,value,RP,DISTANCE,ANN,PMP,z) %>%
reshape2::dcast(data=.,formula=DISTANCE+ANN+PMP+z~RP,value.var = "value") %>%
hydropairs(dec = 2)
site_max_prcp_tb<-data.frame(RP=c(2,5,10,25,50,100,200,500,1000,2000,5000),
day_1=c(92,129,153,183,206,228,250,280,302,324,376))
prcp_max_tbl %>%
reshape2::melt() %>%
mutate(variable=as.numeric(as.character(variable))) %>%
mutate(RP=RP.text(variable,FALSE)) %>%
dplyr::filter(variable>=0.5) %>%
left_join(.,NOAA.PRCP.idx,by=c("site"="NAME")) %>%
ggplot(data=.,aes(x=RP,y=value,col=site))+
geom_point(aes(shape=site))+
geom_point(data=site_max_prcp_tb,aes(x=RP,y=day_1),inherit.aes = F,col="black",size=5)+
geom_line()+
scale_x_log10()+
scale_y_log10()+
scale_shape_manual(values=LETTERS)+
labs(x="Return periods",y="maxumum precipitation 1day[mm]",
caption="Note: black dots from Golder report",
title="Maximum precipitation of 1 days - benchmark comparison")+
theme_light()+
theme(legend.position="bottom")
ggplot()+
geom_spatraster(data = topo_ras, alpha = 0.6,na.rm = T)+
geom_point(data=NOAA.PRCP.idx,aes(x=LONGITUDE,y=LATITUDE,col=prcp_100_yr),size=5)+
geom_text_repel(data=NOAA.PRCP.idx,aes(x=LONGITUDE,y=LATITUDE,label=round(prcp_100_yr,0)),size=2.5)+
geom_point(data=site.info,aes(x = Longitude, y = Latitude), color = "black",
shape = 21, fill = "white",size=5)+
geom_text_repel(data=site.info,aes(x = Longitude, y = Latitude,
label="Golder\n228mm"))+
# scale_fill_continuous()+
scale_color_viridis_c(option = "A")+
scale_fill_binned(type = "viridis",breaks=seq(0,1600,200))+
labs(title="Topography and Max 1 dat 1:100 yrs",
fill="topography\n[masl]",col="max prcp\n100yr [1day]")
8 Annual Precipitation
For frequency analysis information was compiled with a GUMBEL distribution.
All the values seems to be within a regional context.
map_group_z<-merge(
chirps=nc2z_fn(Longitude = Longitude,Latitude = Latitude,
source = "CHIRPS") %>%
window(start=as.Date("2001-01-01"),end=as.Date("2022-12-31")),
"era5_land"=nc2z_fn(Longitude = Longitude,Latitude = Latitude,
source = "ERA5-LAND") %>%
window(start=as.Date("1980-01-01"),end=as.Date("2022-10-31")),
merge=nc2z_fn(Longitude = Longitude,Latitude = Latitude,
source = "MERGE") %>%
window(start=as.Date("1980-01-01"),end=as.Date("2022-12-31")),
NOAA.PRCP.syn$Daily.Comp[,1:6]
)
#note this may change if more of less TS are included
names(map_group_z)[4:9]<-NOAA.PRCP.idx$NAME[1:6]
map_fa_tbl<-lapply(map_group_z,FUN=function(x){
# x<-map_group_z[,1]
sel_yrs<-coredata(dwi(x)>350)
total_yrs<-coredata(daily2annual(x,FUN=sum)[sel_yrs])
FA(total_yrs,distrib = "GUMBEL") %>%
.$results %>%
.[,2]
}) %>%
do.call(rbind,.) %>%
data.frame()
names(map_fa_tbl)<-FA.valcc_report_tbl<-data.frame(RP=c(5,10,20,50,100),Wet=c(1923,2011,2083,2164,2218),
Dry=c(1087,1009,944,871,822)) %>%
reshape2::melt(id.vars="RP") %>%
mutate(variable=paste0(RP," ",variable),
stations="Climate change Golder") %>%
left_join(.,data.frame(prob=FA.val,RP=RP.text(FA.val)),by=c("variable"="RP")) %>%
mutate(type=if_else(prob>=0.5,"wet","dry"))
map_fa_tbl %>%
rownames_to_column(var = "stations") %>%
reshape2::melt(id.vars="stations") %>%
mutate(RP=RP.text(as.numeric(as.character(variable)),FALSE)) %>%
mutate(variable=as.character(variable) %>% as.numeric()) %>%
mutate(type=if_else(variable>=0.5,"wet","dry")) %>%
ggplot(data=.,aes(x=RP,y=value,col=stations,group=stations))+
geom_point(aes(shape=stations),size=5)+
geom_point(data=cc_report_tbl,aes(x=RP,y=value),size=5,col="black")+
geom_line()+
# scale_x_discrete()+
scale_x_log10()+
# scale_x_discrete(breaks=FA.val,label=str_replace_all(RP.text(FA.val)," ","\n"))+
scale_y_log10(limits=c(500,NA),breaks=c(500,800,1000,1500,2000,3000,5000,10000))+
scale_shape_manual(values=LETTERS)+
annotation_logticks(sides="lb")+
theme_light()+
theme(legend.position="bottom")+
labs(x="Return Period",
title="Annual Precipitation - Return Period comparison",
y="annual precipitation [mm/yr]",col="sources:",shape="sources:",
caption="Note: black dots from Golder report")+
facet_wrap(~type)