Meteorological Review
Victor Munoz S.
2023-01-25
1 Library Installation
#Tidyverse
library(plyr); library(dplyr)
library(tidyverse); library(lubridate);library(magrittr)
#Time series
library(zoo);library(xts);library(timetk);library(hydroTSM);
library(dygraphs)
#ggplot support
library(ggmap);library(ggpmisc);library(ggrepel);library(ggsn)
#spreadsheets
library(openxlsx);library(readxl)
#paralel computing
library(parallel);library(pbapply)
#general support
library(here);library(glue);library(udunits2);library(usethis);
library(matrixStats);library(pander);library(reshape2);library(Hydro)
library(geosphere);library(nasapower)
#data table support
library(data.table);library(dtplyr);library(arrow)
#raster support
library(raster);library(sp);library(kgc)
#These directories are implemented to organize the regional analysis
use_directory("data")
use_directory("data/csv")
use_directory("data/netcdf")
use_directory("data/rds")
use_directory("graphs")
use_directory("reports")
use_directory("scripts")
use_directory("outputs")
glue::glue("Work directory: {here()}")2 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= -152.946304
Latitude = 60.122490
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 | -152.9 | 60.12 | 705.2 |
#minimum year of information used, to be defined after seen the available info
min.ywi<-10 3 SOURCES
3.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 = 0.8,
noaa.key = key_noaa)
#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 |
|---|---|---|---|---|---|---|---|
| USC00501818 | CHISIK ISLAND | 60.1 | -152.6 | 9.1 | 1967-01-01 | 1969-12-31 | 20.29 |
| USC00508211 | SARGENT CREEK | 59.97 | -152.7 | 18 | 1964-01-01 | 1965-12-31 | 21.53 |
| USC00508477 | SILVER SALMON CREEK | 59.97 | -152.7 | 46 | 1966-01-01 | 1967-12-31 | 23.29 |
| USR0000ACHG | CHIGMIT MTNS ALASKA | 60.22 | -153.5 | 1372 | 2009-01-01 | 2023-01-26 | 30.91 |
| USC00501810 | CHINITNA BAY | 59.75 | -153.2 | 92 | 1936-01-01 | 1939-12-31 | 44.45 |
| USC00503925 | INISKIN | 59.75 | -153.2 | 92 | 1954-01-01 | 1962-12-31 | 44.45 |
| 994690_99999 | DRIFT RIVER TERMINAL | 60.55 | -152.1 | 15.9 | 2002-02-19 | 2019-04-07 | 65.1 |
| USC00506441 | NINILCHIK | 60.05 | -151.7 | 36.9 | 1940-01-01 | 1968-12-31 | 71.49 |
| USR0000ANIN | NINILCHIK ALASKA | 60.04 | -151.7 | 39.6 | 1992-01-01 | 2023-01-26 | 71.65 |
| USC00506450 | NINILCHIK 5 NE | 60.1 | -151.6 | 64 | 1968-01-01 | 1971-12-31 | 75.65 |
| 703627_99999 | PEDRO BAY | 59.78 | -154.1 | 14 | 2004-07-02 | 2004-09-28 | 76.18 |
| 703406_26546 | PORT ALSWORTH AIRPORT | 60.2 | -154.3 | 79.3 | 2005-01-01 | 2023-01-26 | 76.4 |
| 703406_99999 | PORT ALSWORTH AIR | 60.2 | -154.3 | 85 | 1973-01-01 | 2004-12-31 | 76.41 |
| USC00507570 | PORT ALSWORTH | 60.2 | -154.3 | 79.2 | 1960-01-01 | 2023-01-26 | 76.42 |
| USW00026562 | PORT ALSWORTH 1 SW | 60.2 | -154.3 | 97.8 | 2009-01-01 | 2023-01-26 | 76.5 |
| 999999_26562 | PORT ALSWORTH 1 SW | 60.2 | -154.3 | 97.8 | 2009-09-25 | 2023-01-26 | 76.52 |
| USR0000APAL | PORT ALSWORTH ALASKA | 60.2 | -154.3 | 79.2 | 1992-01-01 | 2023-01-26 | 76.53 |
| USC00509005 | TANALIAN POINT | 60.22 | -154.4 | 92 | 1940-01-01 | 1959-12-31 | 79.35 |
| USC00503672 | HOMER 8 NW | 59.74 | -151.6 | 329.2 | 1977-01-01 | 2023-01-26 | 84.39 |
| 702986_26557 | BIG RIVER LAKE | 60.81 | -152.3 | 18.3 | 2006-01-01 | 2019-08-08 | 84.82 |
| USC00500788 | BIG RIVER LAKES | 60.81 | -152.3 | 18.3 | 1978-01-01 | 2010-12-31 | 84.84 |
| 702986_99999 | BIG RIVER LAKE | 60.82 | -152.3 | 12 | 1978-08-23 | 2005-12-31 | 85.06 |
| 994700_99999 | AUGUSTINE ISLAND | 59.38 | -153.3 | 9.1 | 2002-07-03 | 2023-01-26 | 85.7 |
| USC00509144 | THE TICES | 59.68 | -151.7 | 274.9 | 1973-01-01 | 1975-12-31 | 87.33 |
#used for MERRA2 - small area coordinates +/- 2.5 degrees
xy <-ggmap::make_bbox(lon = c(0.5*round(Longitude/0.5,0)-2.5,
0.5*round(Longitude/0.5,0)+2.5),
lat = c(0.5*round(Latitude/0.5,0)-2.5,
0.5*round(Latitude/0.5,0)+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","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')))
}3.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)
#
# arctic_dem_100m<-raster("V:/Data/arctic_dem/arcticdem_mosaic_100m_v3.0.tif")
#
# bbox_ras<-raster(xmn=xy[1],xmx=xy[3],ymn=xy[2],ymx=xy[4]) %>%
# projectRaster(from = .,crs = crs(arctic_dem_100m))
#
# arctic_dem_cropped<-raster::crop(arctic_dem_100m,bbox_ras) %>%
# projectRaster(from=.,crs=crs("+init=epsg:4326"))
#
# arctic_dem_cropped[arctic_dem_cropped<0]<-NA
#
# raster::writeRaster(x = arctic_dem_cropped,filename = here::here("data","dem","arctic_dem_100m.grd"))
SRTM_ras<-raster(here::here("data","dem","arctic_dem_100m.grd")) %>%
aggregate(.,10)
# 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)
3.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")
3.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")
4 Mean Air Temperature
NOAA.TAVG.idx$NAME<-NOAA.TAVG.idx$NAME %>%make.names(unique=T) %>%
str_replace_all("\\."," ") %>%
str_wrap(string = .,width = 20)
names(NOAA.TAVG.syn$Daily.Comp)<-NOAA.TAVG.idx$NAME
st.names<-glue('{NOAA.TAVG.idx$NAME}\n{round(NOAA.TAVG.idx$ANN,1)} degC')
dwi.graph(x=NOAA.TAVG.syn$"Daily.Comp",
station.name = st.names,
start.year = 1900)
NOAA_tavg_gp<-NOAA.TAVG.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=tavg) %>%
mutate(source="MERRA"),
NOAA.TAVG.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_tavg_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.TAVG.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")
tavg.sel<-1:3
tavg.comp<-merge(NOAA.TAVG.syn$Daily.Comp[,tavg.sel],
MERRA2=merra2_z$tavg)
# tavg.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 temperature")
tavg.comp %>%
daily2monthly(FUN=mean) %>%
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 monthly temperature")
# Definition of the ranges
tavg_ranges<-c(merra2_MA_tbl$tavg,NOAA.TAVG.idx$ANN) %>% na.omit()
tavg_bin<-seq(min(tavg_ranges) %>% floor,
max(tavg_ranges) %>% ceiling,length.out=11) %>%
signif(2) %>% unique()
tavg_lbl<-paste(tavg_bin[1:(length(tavg_bin)-1)] %>%
str_pad(string = .,width = 2,side = "left",pad = "0"),
"–",
tavg_bin[2:(length(tavg_bin))] %>%
str_pad(string = .,width = 2,side = "left",pad = "0"))
#Map
ggmap(map, darken = c(0.3, "white"))+
#MERRA
geom_tile(data = merra2_MA_tbl,aes(x=x,y=y,
fill=cut(tavg,tavg_bin,tavg_lbl)),alpha=0.6)+
geom_label_repel(data=NOAA.TAVG.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.TAVG.idx,
aes(x=LONGITUDE,y=LATITUDE,
fill=cut(ANN,tavg_bin,tavg_lbl)),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\nAir Temperature\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()
5 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")
6 PMP
pmp_tbl<-lapply(NOAA.PRCP.syn$Max.Daily.data,FUN = function(x){
# x<-NOAA.PRCP.syn$Max.Daily.data[[1]]
pmp_res<-x$`Annual Peak` %>%
na.omit() %>%
PMP(val.max=.,input.hrs = 24,number.of.obs.unit = 1,area.km = 1,graphs = F)
return(pmp_res$PMP)
}) %>%
do.call(rbind,.)
NOAA.PRCP.idx %>%
mutate(PMP=pmp_tbl[,1]) %>%
dplyr::select(NAME,LATITUDE,LONGITUDE,DISTANCE,ANN,PMP) %>%
pandoc.table(.,caption="Regional PMP - Hershfield")| NAME | LATITUDE | LONGITUDE | DISTANCE | ANN | PMP |
|---|---|---|---|---|---|
| DRIFT RIVER TERMINAL | 60.55 | -152.1 | 65.1 | 6.758 | 268 |
| PORT ALSWORTH AIR | 60.2 | -154.3 | 76.41 | 0.9737 | 93 |
| PORT ALSWORTH | 60.2 | -154.3 | 76.42 | 370.2 | 383 |
| PORT ALSWORTH 1 SW | 60.2 | -154.3 | 76.5 | 551.9 | 222 |
| PORT ALSWORTH 1 SW 1 | 60.2 | -154.3 | 76.52 | 554.2 | 222 |
| HOMER 8 NW | 59.74 | -151.6 | 84.39 | 762.7 | 329 |
| BIG RIVER LAKES | 60.81 | -152.3 | 84.84 | 1438 | 458 |
| BIG RIVER LAKE | 60.82 | -152.3 | 85.06 | 544.3 | 1074 |
| AUGUSTINE ISLAND | 59.38 | -153.3 | 85.7 | 7.369 | 268 |
#Definition of the ranges
prcp_ranges<-c(merra2_MA_tbl$prcp,NOAA.PRCP.idx$ANN) %>% na.omit()
summary(prcp_ranges)
prcp_bin<-seq(min(prcp_ranges) %>% floor,
max(prcp_ranges) %>% ceiling+50,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 ERA5-land\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()
7 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.
7.1 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 | 72.3 | 60.2 | 39.2 | 38 | 32.7 | 56.1 | 85.7 | 105.7 | 150.1 | 127.9 | 118.9 | 111 | 997.8 |
| rain | 12 | 7.7 | 5.2 | 12.4 | 27.8 | 56.1 | 85.7 | 105.7 | 143.8 | 90.5 | 33.9 | 12.5 | 593.2 |
| snow | 60.3 | 52.5 | 34 | 25.6 | 4.9 | 0 | 0 | 0 | 6.4 | 37.4 | 84.9 | 98.5 | 404.6 |
| tmax | -6 | -2.8 | -2.3 | 3.4 | 9.2 | 13.5 | 14.9 | 13.5 | 8.7 | 2.3 | -3.6 | -4.7 | 3.8 |
| tmin | -13.3 | -10.9 | -12.3 | -5.8 | -0.7 | 3.4 | 6.1 | 5.2 | 1.2 | -3.9 | -10.3 | -11.3 | -4.4 |
| tavg | -9.6 | -6.8 | -7.3 | -1.2 | 4.2 | 8.4 | 10.5 | 9.4 | 5 | -0.8 | -6.9 | -8 | -0.3 |
7.2 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))
8 Output
file.image<-here::here("data","rds",paste0("Backup Meteorological review_",as.Date(now()),".rdata"))
save.image(file = file.image)