Meteorological Review
Victor Munoz S.
2023-06-30
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= -73.962799
Latitude = 46.632005
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 = 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 | -73.96 | 46.63 | 545.5 |
2 SOURCES
2.1 DAYMET
library(daymetr)
library(magrittr)
library(terra)
daymet_info<-daymetr::download_daymet(site = "Daymet",
lat = Latitude,lon = Longitude,
start = 1980,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")
daymet_ras<-rast(here::here("data","netcdf","daymet_prcp_mean.nc"))
terra::plot(daymet_ras,type="interval",main="mean annual precipitation",breaks=seq(800,2000,200))
points(Longitude,Latitude)
2.2 ECCC
library(weathercan)
eccc_files<-here::here("data","rds","eccc_info.rdata")
if(file.exists(eccc_files)){
load(eccc_files)
}else{
weathercan::stations_dl()
eccc_station_tbl<-weathercan::stations() %>%
dplyr::mutate(distance=geosphere::distVincentyEllipsoid(p1=cbind(lon,lat),
p2=cbind(Longitude,Latitude))/1000) %>%
dplyr::select(prov,station_name,station_id,lat,lon,elev,start,end,interval,distance) %>%
dplyr::arrange(distance) %>%
dplyr::filter(interval=="day",distance<150,
end>2000)
eccc_station_db<-weathercan::weather_dl(station_ids=eccc_station_tbl$station_id,
interval = "day")
save(list = c("eccc_station_tbl","eccc_station_db"),
file = here::here("data","rds","eccc_info.rdata"))
}
eccc_prcp_z<-eccc_station_db %>%
dplyr::select(station_id,date,total_precip) %>%
reshape2::dcast(id.vars="date",formula = date~station_id,value.var = "total_precip") %>%
tk_zoo(.,silent=T) %>%
window(start=as.Date("1900-01-01"))
eccc_tavg_z<-eccc_station_db %>%
dplyr::select(station_id,date,mean_temp) %>%
reshape2::dcast(id.vars="date",formula = date~station_id,value.var = "mean_temp") %>%
tk_zoo(.,silent=T) %>%
window(start=as.Date("1900-01-01"))
eccc_prcp_sel1<-eccc_prcp_z %>%
coredata() %>%
colSums2(na.rm=T) %>%
magrittr::is_greater_than(5000)
eccc_prcp_z<-eccc_prcp_z[,eccc_prcp_sel1]
eccc_prcp_station_tbl<-eccc_station_tbl %>%
dplyr::filter(station_id%in%names(eccc_prcp_z))
eccc_prcp_z<-eccc_prcp_z[,eccc_prcp_station_tbl$station_id %>% as.character()]
lat.min<-46
lat.max<-47.5
lon.max<- -72.5
lon.min<- -75.5
eccc_prcp_sel2<-eccc_prcp_station_tbl$lat>=lat.min&
eccc_prcp_station_tbl$lat<=lat.max&
eccc_prcp_station_tbl$lon<=lon.max&
eccc_prcp_station_tbl$lon>=lon.min
eccc_prcp_z<-eccc_prcp_z[,eccc_prcp_sel2]
eccc_prcp_station_tbl<-eccc_station_tbl %>%
dplyr::filter(station_id%in%names(eccc_prcp_z))
eccc_prcp_z<-eccc_prcp_z[,eccc_prcp_station_tbl$station_id %>% as.character()]
eccc_prcp_station_tbl$ANN_all<-daily2annual.avg(eccc_prcp_z,FUN=sum)[,14]
eccc_prcp_station_tbl$ANN_1980<-daily2annual.avg(eccc_prcp_z %>%
window(start=as.Date("1980-01-01")),FUN=sum)[,14]
dwi.graph(eccc_prcp_z,
station.name = eccc_prcp_station_tbl$station_name,start.year = 1900
)
pandoc.table(eccc_prcp_station_tbl %>%
dplyr::select(-prov,-station_id,-interval),split.table=Inf,
caption="Mean annual precipitation / all records / from 1980 in mm")| station_name | lat | lon | elev | start | end | distance | ANN_all | ANN_1980 |
|---|---|---|---|---|---|---|---|---|
| ST MICHEL DES SAINTS | 46.68 | -73.92 | 350.5 | 1966 | 2022 | 6.261 | 927.8 | 920.9 |
| SAINT-MICHEL-DES-SAINTS | 46.82 | -74.09 | 429.9 | 2009 | 2023 | 23.05 | 910.8 | 910.8 |
| ST DONAT | 46.32 | -74.2 | 388.6 | 1964 | 2014 | 39.18 | 1098 | 1101 |
| ST COME | 46.28 | -73.75 | 244 | 1950 | 2021 | 42.41 | 1052 | 1087 |
| STE BEATRIX | 46.2 | -73.6 | 198 | 1974 | 2022 | 55.53 | 1074 | 1073 |
| ST CHARLES DE MANDEVILLE | 46.35 | -73.35 | 167.6 | 1921 | 2022 | 56.53 | 962.4 | 1076 |
| ST ALEXIS DES MONTS | 46.53 | -73.15 | 183 | 1963 | 2022 | 63.32 | 1045 | 1057 |
| LA MACAZA | 46.37 | -74.77 | 243.8 | 1976 | 2022 | 68.46 | 1020 | 1031 |
| ST-JOVITE | 46.08 | -74.56 | 238.5 | 1994 | 2023 | 76.67 | 1123 | 1123 |
| JOLIETTE VILLE | 46.02 | -73.43 | 56 | 1967 | 2011 | 79.44 | 1010 | 1001 |
| LOUISEVILLE | 46.27 | -73.02 | 45 | 1967 | 2022 | 82.86 | 966.6 | 983.3 |
| NOMININGUE | 46.4 | -75.08 | 274 | 1913 | 2013 | 89.53 | 969.9 | 1072 |
| SHAWINIGAN | 46.57 | -72.75 | 121.9 | 1902 | 2004 | 93.18 | 970.8 | 986.7 |
| SOREL | 46.03 | -73.12 | 14.6 | 1914 | 2022 | 93.22 | 957.3 | 1017 |
| GRANDE ANSE | 47.1 | -72.93 | 119 | 1982 | 2015 | 94.38 | 930.7 | 930.7 |
| SHAWINIGAN | 46.56 | -72.73 | 110 | 1998 | 2023 | 94.8 | NA | NA |
| TROIS-RIVIERES | 46.37 | -72.68 | 62.2 | 2004 | 2023 | 102.7 | 1002 | 1002 |
| HEROUXVILLE | 46.67 | -72.6 | 145 | 1966 | 2021 | 104.4 | 1026 | 1014 |
| PIERREVILLE | 46.08 | -72.83 | 15.2 | 1980 | 2022 | 106.6 | 982.2 | 982.2 |
| TROIS RIVIERES AQUEDUC | 46.38 | -72.62 | 54.9 | 1974 | 2009 | 106.8 | 1106 | 1119 |
| STE ANNE DU LAC | 46.85 | -75.33 | 262.1 | 1963 | 2022 | 107.3 | 1039 | 1036 |
| NICOLET | 46.23 | -72.66 | 8 | 1993 | 2023 | 109.6 | 949.6 | 949.6 |
| NICOLET | 46.2 | -72.62 | 30.4 | 1913 | 2018 | 113.9 | 973 | 923 |
| TROIS-RIVIERES | 46.35 | -72.52 | 6 | 1993 | 2023 | 115.1 | 918.4 | 918.4 |
| LA TUQUE | 47.4 | -72.78 | 152 | 1911 | 2004 | 124 | 908.2 | 912.7 |
| LA TUQUE | 47.41 | -72.79 | 168.9 | 1994 | 2023 | 124.2 | 890.7 | 890.7 |
2.3 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.
xy<-ggmap::make_bbox(lon = c(lon.min,lon.max),lat = c(lat.min,lat.max))
map<-ggmap::get_map(xy)
#Download SRTM
# SRTM_ras<-SRTM_capture(bbox = xy,file = "SRTM",
# trim = T,overwrite = T,
# maxcell = 1e5)
SRTM_ras<-raster(here::here("data","dem","SRTM.bil"))
daymet_ras<-crop(daymet_ras,rast(xmin=lon.min,
xmax=lon.max,
ymin=lat.min,
ymax=lat.max))
meteo_ras<-daymet_ras
meteo_ras$z<-resample(rast(SRTM_ras),daymet_ras)
names(meteo_ras)<-c("prcp","z")
plot(rast(SRTM_ras),type="interval",main="elevation",breaks=seq(00,900,100))
points(Longitude,Latitude)
3 Total Precipitation
3.1 Correlation matrix
tidyterra::fortify(meteo_ras) %>%
dplyr::filter(complete.cases(.)) %>%
data.frame() %>%
mutate(dist=geosphere::distVincentyEllipsoid(p1=cbind(Longitude,Latitude),
p2=cbind(x,y))/1000) %>%
dplyr::filter(dist<50) %>%
hydropairs(main="correlation matrix - point 50 km from site")
3.2 Local Geographical Distribution
The main issue is that the site is located 200 m over the elevation of the closest station St Michel Des Saints.
d_y=0.4
d_x=0.5
library(tidyterra)
#Map
ggmap(map, darken = c(0.3, "white"))+
#TOPOGRAHY
geom_spatraster(data=daymet_ras,alpha=0.7,na.rm = T)+
# geom_label_repel(data=eccc_prcp_station_tbl,
# aes(x=lon,y=lat,
# label=paste0(station_name,"\n",round(ANN_1980,1)) %>%
# str_wrap(20)),colour="darkblue",
# fontface="bold",alpha=0.6,
# box.padding =unit(0.5, "lines"),segment.color = "black",
# size=2 )+
#
geom_hline(yintercept =seq(40,50,d_y),linetype="dashed",col="red")+
geom_vline(xintercept =seq(-80.5,-70.5,d_x),linetype="dashed",col="red")+
geom_point(data = eccc_prcp_station_tbl,
aes(x=lon,y=lat,
fill=ANN_1980),size=3,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_binned(type = "viridis",breaks=seq(0,2000,50))+
# scale_fill_brewer(palette="Spectral",direction = -1)+
labs(size="Elevation [masl]",x="Longitude [deg]",y="Latitude [deg]",
shape="Station Name",
fill = "Background: \n * DAYMET\nPoints:\n * ECCC",
caption=
"Note: Site is presented as a red diamond",
title="Mean Annual Precipitation",
subtitle="(1980 - 2022)")+
theme(legend.position="right")+
ggsn::scalebar(x.min=-75.5,x.max=-72.5,y.min=46,y.max=47.5,dist = 20,
transform=T,
model = "WGS84",dist_unit = "km")
ggmap(map, darken = c(0.3, "white"))+
#TOPOGRAHY
geom_spatraster(data=rast(SRTM_ras),alpha=0.7,na.rm = T)+
# geom_label_repel(data=eccc_prcp_station_tbl,
# aes(x=lon,y=lat,
# label=paste0(station_name,"\n",round(ANN_1980,1)) %>%
# str_wrap(20)),colour="darkblue",
# fontface="bold",alpha=0.6,
# box.padding =unit(0.5, "lines"),segment.color = "black",
# size=2 )+
#
geom_hline(yintercept =seq(40,50,d_y),linetype="dashed",col="red")+
geom_vline(xintercept =seq(-80.5,-70.5,d_x),linetype="dashed",col="red")+
geom_point(data = eccc_prcp_station_tbl,
aes(x=lon,y=lat,
fill=elev),size=3,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_binned(type = "viridis",breaks=seq(0,2000,100))+
# scale_fill_brewer(palette="Spectral",direction = -1)+
labs(size="Elevation [masl]",x="Longitude [deg]",y="Latitude [deg]",
shape="Station Name",
fill = "Elevation [masl]:",
caption=
"Note: Site is presented as a red diamond",
title="Elevation - SRTM")+
theme(legend.position="right")+ ggsn::scalebar(x.min=-75.5,x.max=-72.5,y.min=46,y.max=47.5,dist = 20,
transform=T,
model = "WGS84",dist_unit = "km")
3.3 Relationships with Elevation / Longitude / Latitude
site_info_tbl<-data.frame(lon=rep(Longitude,2),
lat=rep(Latitude,2),
elev=rep(Elevation,2)) %>%
mutate(y_bin=cut(lat,seq(40,50,d_y))) %>%
mutate(x_bin=cut(lon,seq(-80.5,-70.5,d_x))) %>%
mutate(MAP=c(daymet_z$prcp %>%
daily2annual.avg(FUN=sum) %>% .[,13],
949),
source=c("Daymet","Previous work"))
fortify(meteo_ras) %>%
dplyr::filter(complete.cases(.)) %>%
data.frame() %>%
mutate(dist=geosphere::distVincentyEllipsoid(p1=cbind(Longitude,Latitude),
p2=cbind(x,y))/1000) %>%
mutate(y_bin=cut(y,seq(40,50,d_y))) %>%
mutate(x_bin=cut(x,seq(-80.5,-70.5,d_x))) %>%
ggplot(data=.,aes(x=z,y=prcp))+
geom_point(col="grey")+
geom_point(data=site_info_tbl,
aes(x=elev,y=MAP,col=source),size=4)+
geom_smooth(method="lm")+
ggpmisc::stat_poly_eq(use_label(c("eq", "adj.R2")),
formula = y ~ poly(x, 1, raw = TRUE),size=2.5)+
geom_point(data=eccc_prcp_station_tbl %>%
mutate(y_bin=cut(lat,seq(40,50,d_y))) %>%
mutate(x_bin=cut(lon,seq(-80.5,-70.5,d_x))),
aes(x=elev,y=ANN_1980),col="red")+
geom_text_repel(data=eccc_prcp_station_tbl %>%
mutate(y_bin=cut(lat,seq(40,50,d_y))) %>%
mutate(x_bin=cut(lon,seq(-80.5,-70.5,d_x))),
aes(x=elev,y=ANN_1980,label=station_name %>% str_wrap(10)),col="red",size=2.5)+
geom_vline(xintercept = Elevation,linetype="dashed")+
facet_grid(y_bin~x_bin)+
labs(col="site source for\nprecipitation:",
title="Mean annual precipitation versus elevation (1980 - 2022)",
subtitle="Based on grid 0.5 x 0.5 (deg) - grey points: Daymet / red points: ECCC",
x="Elevation [masl]",y="Mean annual precipitation [mm/yr]",
caption="note: dotted line represent site elevation")+
theme_light()+
theme(legend.position="bottom")
fortify(meteo_ras) %>%
dplyr::filter(complete.cases(.)) %>%
data.frame() %>%
mutate(dist=geosphere::distVincentyEllipsoid(p1=cbind(Longitude,Latitude),
p2=cbind(x,y))/1000) %>%
mutate(y_bin=cut(y,seq(40,50,d_y))) %>%
mutate(x_bin=cut(x,seq(-80.5,-70.5,d_x))) %>%
ggplot(data=.,aes(x=z,y=prcp))+
geom_point(col="grey")+
geom_point(data=site_info_tbl,
aes(x=elev,y=MAP,col=source),size=4)+
geom_vline(xintercept = Elevation,linetype="dashed")+
geom_smooth(method="loess",se=F)+
ggpmisc::stat_poly_eq(use_label(c("eq", "adj.R2")),
formula = y ~ poly(x, 1, raw = TRUE),size=3)+
geom_point(data=eccc_prcp_station_tbl %>%
mutate(y_bin=cut(lat,seq(40,50,d_y))) %>%
mutate(x_bin=cut(lon,seq(-80.5,-70.5,d_x))),
aes(x=elev,y=ANN_1980),col="red")+
geom_text_repel(data=eccc_prcp_station_tbl %>%
mutate(y_bin=cut(lat,seq(40,50,d_y))) %>%
mutate(x_bin=cut(lon,seq(-80.5,-70.5,d_x))),
aes(x=elev,y=ANN_1980,label=station_name),col="red",size=2)+
# facet_grid(y_bin~x_bin)+
labs(col="site source for\nprecipitation:",
title="Mean annual precipitation versus elevation (1980 - 2022)",
# subtitle="Based on grid 0.5 x 0.5 (deg) - grey points: Daymet / red points: ECCC",
x="Elevation [masl]",y="Mean annual precipitation [mm/yr]",
caption="note: dotted line represents site elevation")+
theme_light()+
theme(legend.position="bottom")
4 Climograph
Information at the site based on Daymet. As a surface water analysis, you have to consider this work as a simple overview; which assumes Daymet is representative of this location.
This analysis does not include or considers snow capture efficiency; which it should change or modified local snowfall estimations.
# source_gp<-"Daymet"
# 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 | 78.2 | 64.9 | 76.3 | 86.2 | 100.6 | 123.2 | 119.4 | 114.8 | 117.1 | 106.2 | 91.3 | 91.5 | 1170 |
| rain | 8 | 5.3 | 21.1 | 62.7 | 98.9 | 123.2 | 119.4 | 114.8 | 116.7 | 97.7 | 44.1 | 13.4 | 825.2 |
| snow | 70.2 | 59.7 | 55.1 | 23.5 | 1.7 | 0 | 0 | 0 | 0.4 | 8.5 | 47.3 | 78.2 | 344.6 |
| tmax | -8.2 | -5.6 | 0.2 | 7.5 | 16 | 20.7 | 23.1 | 21.9 | 17 | 9.4 | 1.7 | -5.1 | 8.2 |
| tmin | -20.7 | -19.7 | -13.4 | -4.5 | 2.4 | 7.7 | 10.4 | 9.6 | 5.2 | -0.4 | -6.8 | -15.6 | -3.8 |
| tavg | -14.4 | -12.7 | -6.6 | 1.5 | 9.2 | 14.2 | 16.7 | 15.7 | 11.1 | 4.5 | -2.6 | -10.4 | 2.2 |