TwoRivers_Climate
Mehmetcan Ozkadioglu & Victor Munoz
2022-11-30
1 Precipitation
This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown.
# Locations of the Projects
locs <- read.xlsx(here::here("data","csv","Two_Rivers_Locations.xlsx"))
project_loc <- locs
project_loc <- project_loc %>%
filter(Name == 'De Grooteboom TSF')# Setting Coordinates
Left_P1 <- 28.5 # East
Right_P1 <- 33.5 # West
Top_P1 <- -23 # North
Bottom_P1 <- (Top_P1 - (Right_P1 - Left_P1-0.5)) # To make it Square (Check Southern Hemisphere)
# Downloading the basemap tiles
base1 <- get_map(c(left = Left_P1,
right = Right_P1,
top = Top_P1,
bottom = Bottom_P1),
zoom = 8,
maptype="terrain",
#scale = 1,
)
# Creating the Map
map1 <- ggmap(base1,
#extent = "normal",
extent = "device",
#extent = "panel",
)
map1 +
geom_point(data = project_loc, aes(x = Lon, y = Lat, fill = Type, shape = Type), color="white", cex=5) +
scale_fill_manual(values = c("red"), labels=c("Project Location"), name=NULL) +
scale_shape_manual(values = c(22), labels=c("Project Location"), name=NULL) +
labs(title="Two Rivers Project",
subtitle = "Locality Map",
x="Latitude (X-Coords)",
y="Longitude (Y-Coords)") +
theme_bw() +
theme(plot.title = element_text(size = 14, face="bold"),
legend.position="right",
axis.text = element_text(size = 12),
axis.title=element_text(size=12.5,face="bold")) +
annotate('rect', xmin=28.95, ymin=-24.67, xmax=30.05, ymax=-24.55, col="navy", fill="white", alpha = 0.7) +
annotate('text', x=29.5, y=-24.6, label = 'Two Rivers Project', colour = I('red'), size = 5) +
annotate('segment', x=30, xend=project_loc$Lon, y=-24.67, yend=project_loc$Lat, colour=I('navy'), arrow = arrow(length=unit(0.1,"cm")), size = 1)
###############################################################################
# Mapping the TSF and SAWS
###############################################################################
# Setting Coordinates
Left_P2 <- 29.5 # East
Right_P2 <- 30.8 # West
Top_P2 <- -24.5 # North
Bottom_P2 <- (Top_P2 - (Right_P2 - Left_P2-0.5)) # To make it Square (Check Southern Hemisphere)
# Downloading the basemap tiles
base2 <- get_map(c(left = Left_P2,
right = Right_P2,
top = Top_P2,
bottom = Bottom_P2),
zoom = 11,
maptype="terrain",
#scale = 1,
)
# Creating the Map
map2 <- ggmap(base2,
#extent = "normal",
extent = "device",
#extent = "panel",
)
df <- read_csv(here::here("data","csv","Two_Rivers_SAWS_Table.csv"))
df
TSF <- locs %>%
filter(Type == "TSF")
De_Grooteboom_TSF <- TSF %>%
filter(Name == "De Grooteboom TSF")
Old_TSF <- TSF %>%
filter(Name == "Old TSF")
#Calculating SAWS Distance to Project TSF's
df <- df %>%
mutate(Dist_De_Grooteboom_TSF = distVincentyEllipsoid(cbind(Lon, Lat),
c(De_Grooteboom_TSF$Lon, De_Grooteboom_TSF$Lat),
a=6378137, b=6356752.3142, f=1/298.257223563)/1000,
Dist_Old_TSF = distVincentyEllipsoid(cbind(Lon, Lat), c(Old_TSF$Lon, Old_TSF$Lat),
a=6378137, b=6356752.3142, f=1/298.257223563)/1000)
# clipr::write_clip(df)
MAP_Plot <- map2 +
geom_point(data = df, mapping = aes(x = Lon,
y = Lat,
fill = MAP,
#label = Name,
),
shape = 21,
size = 5) +
geom_text_repel(data = df, aes(x = Lon,
y = Lat,
label = Name,
box.padding = 0.6,
)) +
scale_fill_viridis_c(option = "G",
direction = -1,
aesthetics = "fill",
begin = 0.35,
end = 1) +
geom_point(data = TSF, aes(x = Lon, y = Lat), shape = 22, size = 4.5, fill = "red") +
geom_label_repel(data = TSF, aes(x = Lon,
y = Lat,
label = Name,
box.padding = 1)) +
theme_bw() +
labs(x="Longitude (X-Coords)",
y="Latitude (Y-Coords)",
fill = "MAP (mm/yr)") +
theme(#plot.title = element_text(size = 14, face="bold"),
#legend.position="bottom",
axis.text = element_text(size = 10.5),
axis.title=element_text(size=11,face="bold"))
MAP_Plot
ggsave(here::here("outputs","TwoRivers_MAP_Stations_V2.png"),
width = 16,
height = 10,
units = c("cm"),
plot = MAP_Plot,
dpi = 600)2 Victor’s review
library(Hydro)
Longitude<-locs$Lon[10]
Latitude<-locs$Lat[10]
#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)
#minimum year of information used, to be defined after seen the available info
min.ywi<-10 2.1 NOAA Review
NOAA records were captured and downloaded
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,
#Get your own key from https://www.ncdc.noaa.gov//cdo-web/token
#Key from vmunoz@srk.co.uk
noaa.key = key_noaa)
#work just with the first 10 stations. Ranges can change based on the project
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"),
stringr::str_remove_all(as.Date(lubridate::now()),"-") %>% as.integer(),
selected.NOAA.GSOD$END)
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 |
|---|---|---|---|---|---|---|---|
| SF005934190 | MARTENSHOOP (POL) | -24.98 | 30.23 | 1408 | 1909-01-01 | 1997-12-31 | 12.63 |
| SF005937780 | KLIPFONTEIN | -24.97 | 30.43 | 1320 | 1936-01-01 | 1997-12-31 | 32.56 |
| SF005924740 | NEBO | -24.92 | 29.77 | 1570 | 1906-01-01 | 1997-12-31 | 34.32 |
| SF005547860 | LYDENBURG | -25.1 | 30.47 | 1425 | 1904-01-01 | 2023-03-14 | 39.62 |
| 681840_99999 | LYDENBURG | -25.1 | 30.47 | 1425 | 1962-08-01 | 2023-03-14 | 39.62 |
| 681850_99999 | LYDENBURG | -25.1 | 30.47 | 1434 | 1998-03-01 | 2023-03-14 | 39.62 |
| SF005931260 | MAANDAGSHOEK | -24.6 | 30.08 | 975 | 1924-01-01 | 1993-12-31 | 39.75 |
| SF005538590 | TONTELDOOS | -25.32 | 30 | 1771 | 1904-01-01 | 1980-12-31 | 41.93 |
| SF005941410 | RUSTPLAATS | -24.85 | 30.58 | 1341 | 1904-01-01 | 1997-12-31 | 49.13 |
| SF005541750 | DULLSTROOM (SAR) | -25.42 | 30.1 | 2030 | 1906-01-01 | 1993-12-31 | 51.64 |
| SF005536510 | LAERSDRIFT (POL) | -25.37 | 29.83 | 1450 | 1923-01-01 | 1997-12-31 | 53.9 |
| SF005554550 | CEYLON (BOS) | -25.1 | 30.75 | 1030 | 1931-01-01 | 1997-12-31 | 66.74 |
| SF005944440 | PILGRIMS REST (POL) | -24.92 | 30.77 | 1240 | 1907-01-01 | 1992-12-31 | 66.98 |
| SF005945390 | MAC MAC BOS | -24.98 | 30.82 | 1360 | 1914-01-01 | 1992-12-31 | 71.93 |
| SF005945900 | BLYDE (BOS) | -24.83 | 30.83 | 1433 | 1926-01-01 | 1997-12-31 | 74.28 |
| SF005554410 | RIETVALLEI | -25.33 | 30.72 | 1030 | 1915-01-01 | 1997-12-31 | 74.43 |
| 682870_99999 | GRASKOP | -24.93 | 30.85 | 1436 | 1998-03-01 | 2023-03-14 | 74.97 |
| SF005944570 | ELANDSFONTEIN | -24.62 | 30.77 | 1219 | 1906-01-01 | 1995-12-31 | 76.69 |
| SF005174300 | MACHADODORP | -25.67 | 30.25 | 1536 | 1903-01-01 | 1997-12-31 | 80.74 |
| SF005178160 | ELANDSHOEK | -25.6 | 30.48 | 1075 | 1914-01-01 | 1997-12-31 | 80.88 |
| 682720_99999 | OUDESTAD | -25.18 | 29.33 | 949 | 1988-11-30 | 2023-03-14 | 82.09 |
| 682850_99999 | BELFAST | -25.69 | 30.03 | 1879 | 2013-12-10 | 2019-07-16 | 82.25 |
| SF005556730 | WITKLIP (BOS) | -25.23 | 30.88 | 1060 | 1937-01-01 | 1997-12-31 | 83.61 |
| SF005948060 | WILGEBOOM (BOS) | -24.93 | 30.95 | 945 | 1935-01-01 | 1997-12-31 | 85.07 |
| SF005181860 | VLAKPLAATS | -25.6 | 30.6 | 961 | 1922-01-01 | 1997-12-31 | 87.13 |
| SF005177620 | WELTEVREDEN | -25.7 | 30.43 | 1844 | 1919-01-01 | 1997-12-31 | 88.93 |
| 681844_99999 | MAIEPSKOP(SAAFB) | -24.55 | 30.87 | 0 | 1988-01-06 | 1988-09-13 | 89.09 |
| SF005946960 | SALIQUE (BOS) | -24.6 | 30.9 | 808 | 1935-01-01 | 1997-12-31 | 89.36 |
| SF005165540 | ROODEPOORT | -25.73 | 29.82 | 1707 | 1905-01-01 | 1997-12-31 | 90.87 |
| SF005555670 | ALKMAAR | -25.45 | 30.83 | 706 | 1905-01-01 | 1997-12-31 | 91.19 |
| SF005558780 | WITWATER (BOS) | -25.13 | 31 | 1042 | 1929-01-01 | 1993-12-31 | 92.05 |
| SF006363080 | THE DOWNS | -24.13 | 30.18 | 1350 | 1913-01-01 | 1975-12-31 | 92.26 |
| SF006365180 | SCHELM | -24.13 | 30.3 | 840 | 1920-01-01 | 1993-12-31 | 94.01 |
| SF005167080 | WONDERFONTEIN (SKL) | -25.8 | 29.9 | 1789 | 1904-01-01 | 1981-12-31 | 96.23 |
| SF005951100 | BOSBOKRAND (POL) | -24.83 | 31.07 | 880 | 1905-01-01 | 1996-12-31 | 98.18 |
| SF005162850 | GEMSBOKFONTEIN | -25.75 | 29.67 | 1638 | 1913-01-01 | 1974-12-31 | 98.73 |
| 682880_99999 | NELSPRUIT | -25.43 | 30.98 | 674 | 1973-01-01 | 2006-02-19 | 102.9 |
| 682890_99999 | NELSPRUIT | -25.43 | 30.98 | 883 | 1998-03-01 | 2023-03-14 | 102.9 |
| 682920_99999 | FLEUR DE LYS | -24.53 | 31.03 | 620 | 1974-10-08 | 2004-05-17 | 104.7 |
| SF005951610 | CHAMPAGNE (NAF) | -24.68 | 31.1 | 549 | 1924-01-01 | 1997-12-31 | 104.9 |
| SF006344170 | UITZICHT | -24.45 | 29.23 | 945 | 1922-01-01 | 1989-12-31 | 105.1 |
| 686380_99999 | MIDDELBURG | -25.68 | 29.44 | 1489 | 1957-01-01 | 2003-04-07 | 105.3 |
| 682770_99999 | MIDDELBURG DAM | -25.77 | 29.55 | 1527 | 1998-03-01 | 2005-10-09 | 106.3 |
| SF005560880 | MAYFERN | -25.48 | 31.03 | 579 | 1930-01-01 | 1997-12-31 | 109.7 |
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("2005-01-01"),
end=as.Date("2022-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 Topography
The topographical information was captured from SRTM source, the original information is at 30 x 30 m, but the topography as upscaled up to 620 x 620 m.
TSF is presented as point.
library(terra)
library(Hydro)
library(ggmap)
# srtm_ras<-SRTM_capture(bbox = bbox2,overwrite = T,trim = T,maxcell = 1e5)
srtm_ras<-rast(here::here("data","dem","topo_SRTM.bil"))
plot(srtm_ras)
points(x = Longitude,y=Latitude,add=TRUE )
2.3 Precipitation
2.3.1 Summary
chirps_files<-dir(here::here("data","netcdf","chirps"),full.names = T)
chirps_prcp<-mean(terra::tapp(rast(chirps_files),index="years",fun=sum)) %>%
tidyterra::fortify() %>%
dplyr::filter(x>=29.8&x<=30.4,
y<=-24.6)
bbox<-ggmap::make_bbox(lon = c(chirps_prcp$x),lat=c(chirps_prcp$y))
base5 <- get_map(location = bbox,
zoom = 10,
maptype="terrain",
#scale = 1,
)
prcp_comp <- read_excel(here::here("data","Copy of 588618_Two_Rivers_MET_Data 02Mar2023.xlsx"),
sheet = "Precipitation", col_types = c("date",
"numeric", "numeric", "numeric",
"numeric", "numeric", "skip", "skip",
"skip", "skip", "skip", "skip", "skip",
"skip"))
map_prcp_loc<-data.frame(Glencore=prcp_comp$`RAIN OBS`,
Two_river=prcp_comp$`RAIN MODELLED`,
time=prcp_comp$DATES) %>%timetk::tk_zoo(silent=T) %>%
daily2annual.avg()
pandoc.table(map_prcp_loc,split.table=Inf,round=1,
caption="Comparison of records")| Station | Jan | Feb | Mar | Apr | May | Jun | Jul | Aug | Sep | Oct | Nov | Dec | Ann |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Glencore | 74.8 | 62 | 50.5 | 30.2 | 10.1 | 0.3 | 0.9 | 1.3 | 12.4 | 42.1 | 106.8 | 108.9 | 500.5 |
| Two_river | 128.3 | 141.3 | 127.3 | 125.9 | 67.7 | 34.8 | 12.8 | 16.5 | 4.6 | 6.4 | 26.1 | 67.8 | 759.5 |
2.3.2 Figures
ggmap(base5) +
geom_tile(data=chirps_prcp,aes(x=x,y=y,fill=mean),alpha=0.8)+
geom_point(data = df, mapping = aes(x = Lon,
y = Lat,
fill = MAP,
#label = Name,
),shape = 21,size = 5) +
geom_text_repel(data = df, aes(x = Lon,
y = Lat,
label = Name,
box.padding = 0.6,
)) +
# geom_point(data=NOAA.PRCP.idx,aes(x=LONGITUDE,y=LATITUDE))+
geom_point(data = TSF, aes(x = Lon, y = Lat,
fill=map_prcp_loc$Ann[2]),
shape = 22, size = 4.5) +
geom_label_repel(data = TSF, aes(x = Lon,
y = Lat,
label = Name,
box.padding = 1)) +
geom_point(aes(x=30.12228333,y=-24.81582777,
fill=map_prcp_loc$Ann[1]),size=5,shape=22)+
geom_label(aes(x=30.12228333,y=-24.81582777,label="Glencore Lion Weather"),
nudge_x = 0.02,nudge_y = 0.02)+
theme_bw() +
labs(x="Longitude (X-Coords)",y="Latitude (Y-Coords)",
fill = "MAP (mm/yr)",
title="Mean annual precipitation",
subtitle="Background: CHIRPS, Points: South African DB / point of interest")+
scale_fill_stepsn(n.breaks=6,#breaks=c(0,10,20,50,100,200,500,1000,1500),
colours = scales::brewer_pal(palette = "Spectral")(5))+
ggsn::scalebar(x.min = 29.8,x.max=30.4,y.min=-25.25,y.max=-24.6,dist = 5,
dist_unit = "km",transform = T,model = "WGS84")
ggsave(filename = here::here("outputs","mean annual precipitation.png"),
width = 10,height = 10,units = "in",dpi = 200)2.4 Temperature
2.4.1 Reading from ERA5-Land
era5_land_temp_files<-dir(here::here("data","netcdf","era5-land"),full.names = T)
#
# era5_land_dly_avg<-terra::tapp(rast(era5_land_temp_files),index="days",fun=mean)
# era5_land_dly_max<-terra::tapp(rast(era5_land_temp_files),index="days",fun=max)
# era5_land_dly_min<-terra::tapp(rast(era5_land_temp_files),index="days",fun=min)
# #
# time(era5_land_dly_avg)<-str_remove_all(names(era5_land_dly_avg),"X") %>% ymd()
# time(era5_land_dly_max)<-str_remove_all(names(era5_land_dly_max),"X") %>% ymd()
# time(era5_land_dly_min)<-str_remove_all(names(era5_land_dly_min),"X") %>% ymd()
#
# terra::writeCDF(x = era5_land_dly_max,filename = here::here("data","netcdf","era5-land_t2m_max.nc"),
# overwrite=TRUE)
#
# terra::writeCDF(x = era5_land_dly_min,filename = here::here("data","netcdf","era5-land_t2m_min.nc"),
# overwrite=TRUE)
#
# terra::writeCDF(x = era5_land_dly_avg,filename = here::here("data","netcdf","era5-land_t2m_avg.nc"),
# overwrite=TRUE)
#Load from Kelvin to degrees
era5_land_dly_ras_avg<-rast(here::here("data","netcdf","era5-land_t2m_avg.nc"))-273.15
era5_land_dly_ras_max<-rast(here::here("data","netcdf","era5-land_t2m_max.nc"))-273.15
era5_land_dly_ras_min<-rast(here::here("data","netcdf","era5-land_t2m_min.nc"))-273.15
era5_land_mon_ras_t2m<-c(mean(tapp(era5_land_dly_ras_max,index="years",fun=mean,na.rm=T)),
mean(tapp(era5_land_dly_ras_avg,index="years",fun=mean,na.rm=T)),
mean(tapp(era5_land_dly_ras_min,index="years",fun=mean,na.rm=T)))
era5_land_mon_ras_t2m<-c(terra::resample(srtm_ras,era5_land_mon_ras_t2m),
era5_land_mon_ras_t2m)
names(era5_land_mon_ras_t2m)<-c("z","max","avg","min")
era5_land_t2m<-era5_land_mon_ras_t2m %>%
tidyterra::fortify()
plot(era5_land_mon_ras_t2m)
t2m_comp <- read_excel(here::here("data","Copy of 588618_Two_Rivers_MET_Data 02Mar2023.xlsx"),
sheet = "Temperature", col_types = c("date",
"numeric", "numeric", "numeric",
"numeric", "numeric", "numeric",
"numeric", "numeric", "numeric",
"numeric", "numeric", "text", "text",
"text", "text", "text", "text", "text"))[,1:10]
t2m_comp_z1<-zoo(t2m_comp[,5:10],as.Date(t2m_comp$DATES)) %>%
# window(end=as.Date("2021-12-31")) %>%
subdaily2daily(FUN=mean)
plot(t2m_comp_z1,main="Daily temperature comparison - all records")
t2m_comp_z<-zoo(t2m_comp[,5:10],as.Date(t2m_comp$DATES)) %>%
window(end=as.Date("2021-12-31")) %>%
subdaily2daily(FUN=mean)
plot(t2m_comp_z,main="Daily temperature comparison - up to 2021-12-31")
# plot(zoo(t2m_comp[,5:10],as.Date(t2m_comp$DATES)))
t2m_comp_ann<-daily2annual.avg(t2m_comp_z,FUN=mean)%>%
separate(col = Station,sep = " ",into = c("par","location"))
pandoc.table(t2m_comp_ann,split.table=Inf,round=1,
caption="Comparison of monthly records. OBS is Glencore / MERRA2 is Two Rivers")| par | location | Jan | Feb | Mar | Apr | May | Jun | Jul | Aug | Sep | Oct | Nov | Dec | Ann |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Tmean | OBS | 25.7 | 25.4 | 24.7 | 21.7 | 18.8 | 15.8 | 15.8 | 18.7 | 21.8 | 23.3 | 24.6 | 25.2 | 21.8 |
| Tmax | OBS | 32.4 | 32.3 | 31.4 | 28.5 | 26.8 | 24.2 | 23.9 | 26.8 | 29.8 | 30.2 | 31.3 | 31.5 | 29.1 |
| Tmin | OBS | 20 | 19.9 | 18.2 | 14.7 | 10.4 | 7.2 | 7.2 | 10.3 | 14 | 16.4 | 18.2 | 19.7 | 14.7 |
| Tmean | MERRA2 | 22 | 22.3 | 21.1 | 19.3 | 17.4 | 14.4 | 12.7 | 13.2 | 16.3 | 19 | 20.1 | 21.5 | 18.3 |
| Tmax | MERRA2 | 29.2 | 29.8 | 28.2 | 26.9 | 25.7 | 23.3 | 21.8 | 22.7 | 25.7 | 28.4 | 28.7 | 29.3 | 26.6 |
| Tmin | MERRA2 | 15.6 | 15.7 | 15.1 | 13 | 10.9 | 7.6 | 5.9 | 6 | 8.7 | 11.3 | 12.7 | 14.7 | 11.4 |
TSF$Elev<-mapply(GM_elev,TSF$Lon,TSF$Lat,key_google = key_google)
Glencore_tbl<-tibble(Lon=30.12228333,Lat=-24.81582777,
Name="Glencore Lion Weather") %>%
mutate(Elev=GM_elev(Lon,Lat,key_google = key_google))
Glencore_tbl<-cbind(Glencore_tbl,t2m_comp_ann %>%
dplyr::filter(location=="OBS") %>%
dplyr::select(par,Ann) %>%
reshape2::dcast(formula = .~par,value.var = "Ann")) %>%
dplyr::select(-.)
TSF_tbl<-cbind(TSF,t2m_comp_ann %>%
dplyr::filter(location=="MERRA2") %>%
dplyr::select(par,Ann) %>%
reshape2::dcast(formula = .~par,value.var = "Ann")) %>%
dplyr::select(-.)
sel_loc<-rbind(Glencore_tbl,
TSF_tbl[names(Glencore_tbl)]) %>%
reshape2::melt(id.vars=c("Lon","Lat","Name","Elev")) %>%
mutate(variable=str_remove_all(variable,"T"),
variable=if_else(variable=="mean","avg",variable))
noaa_loc<-rbind(
NOAA.TAVG.idx %>%
mutate(par="avg"),
NOAA.TMIN.idx %>%
mutate(par="min"),
NOAA.TMAX.idx %>%
mutate(par="max"))
review_loc<-rbind(
sel_loc %>%
mutate(source="reviewed\nstations"),
noaa_loc %>%
dplyr::select(Lon=LONGITUDE,Lat=LATITUDE,Name=NAME,Elev=ELEVATION,
variable=par,value=ANN) %>%
mutate(source="noaa")) %>%
mutate(Name=toupper(Name))
pandoc.table(review_loc,split.table=Inf)| Lon | Lat | Name | Elev | variable | value | source |
|---|---|---|---|---|---|---|
| 30.12 | -24.82 | GLENCORE LION WEATHER | 810.1 | max | 29.1 | reviewed stations |
| 30.14 | -24.93 | DE GROOTEBOOM TSF | 1000 | max | 26.63 | reviewed stations |
| 30.11 | -24.96 | OLD TSF | 981 | max | 26.63 | reviewed stations |
| 30.12 | -24.82 | GLENCORE LION WEATHER | 810.1 | avg | 21.8 | reviewed stations |
| 30.14 | -24.93 | DE GROOTEBOOM TSF | 1000 | avg | 18.26 | reviewed stations |
| 30.11 | -24.96 | OLD TSF | 981 | avg | 18.26 | reviewed stations |
| 30.12 | -24.82 | GLENCORE LION WEATHER | 810.1 | min | 14.68 | reviewed stations |
| 30.14 | -24.93 | DE GROOTEBOOM TSF | 1000 | min | 11.44 | reviewed stations |
| 30.11 | -24.96 | OLD TSF | 981 | min | 11.44 | reviewed stations |
| 30.47 | -25.1 | LYDENBURG | 1434 | avg | 16.74 | noaa |
| 30.85 | -24.93 | GRASKOP | 1436 | avg | 15.77 | noaa |
| 29.33 | -25.18 | OUDESTAD | 949 | avg | 20.47 | noaa |
| 30.98 | -25.43 | NELSPRUIT | 883 | avg | 18.98 | noaa |
| 30.47 | -25.1 | LYDENBURG | 1434 | min | 10.49 | noaa |
| 30.85 | -24.93 | GRASKOP | 1436 | min | 11.07 | noaa |
| 29.33 | -25.18 | OUDESTAD | 949 | min | 12.64 | noaa |
| 30.98 | -25.43 | NELSPRUIT | 883 | min | 13.84 | noaa |
| 30.47 | -25.1 | LYDENBURG | 1434 | max | 25.37 | noaa |
| 30.85 | -24.93 | GRASKOP | 1436 | max | 23.24 | noaa |
| 29.33 | -25.18 | OUDESTAD | 949 | max | 30.62 | noaa |
| 30.98 | -25.43 | NELSPRUIT | 883 | max | 26.37 | noaa |
bbox2<-ggmap::make_bbox(lon = c(era5_land_t2m$x),lat=c(era5_land_t2m$y),f = 0)
base6 <- get_map(location = bbox2,
zoom = 8,
maptype="terrain",
#scale = 1,
)
# sel_loc %>% dplyr::select(Lon,Lat,Name,Elev) %>%
# unique() %>% write.excel()2.4.2 Machine Learning
ml_file<-here::here("outputs","air temperature raster_ML.nc")
if(file.exists(ml_file)){
era5_land_inter_ras<-rast(here::here("outputs",
"air temperature raster_ML.nc"))
names(era5_land_inter_ras)<-c("z","avg","max","min")
}else{
library(MACHISPLIN)
#Topographical information
srtm_ras<-rast(here::here("data","dem","topo_SRTM.bil"))
#ground parameters estimated
topo_terrain_crop.ras<-terra::terrain(srtm_ras,
v = c( "TRI","TPI","roughness","slope","aspect"),
unit = "degrees")
#information implemented as raster (instead of terra!)
raster_stack<-stack(raster(srtm_ras),
brick(topo_terrain_crop.ras)) %>%
raster::crop(x=.,y=raster(era5_land_mon_ras_t2m))
map_utm.df<-era5_land_mon_ras_t2m %>%
tidyterra::fortify() %>%
dplyr::select(x,y,max,avg,min) %>%
data.frame()
#Machine learning model
interp.rast<-machisplin.mltps(int.values=map_utm.df, tps = T,
covar.ras=raster_stack,
smooth.outputs.only=TRUE, n.cores=16)
#rename the layers as were defined on the original dataset
names(interp.rast)<-c("max","avg","min")
era5_land_inter_ras<-c(rast(interp.rast$avg$final),
rast(interp.rast$max$final),
rast(interp.rast$min$final))
era5_land_inter_ras<-c(terra::resample(srtm_ras,era5_land_inter_ras),
era5_land_inter_ras)
names(era5_land_inter_ras)[1]<-"z"
terra::writeCDF(x = era5_land_inter_ras,
filename = here::here("outputs",
"air temperature raster_ML.nc"))
# NOAA.TAVG.syn$Daily.Comp %>% dwi.graph(start.year = 1900)
# tidyterra::fortify(era5_land_inter_ras) %>%
# dplyr::filter(complete.cases(.)) %>%
# reshape2::melt(id.var=c("x","y","topo_SRTM")) %>%
# ggplot(data=.,aes(x=topo_SRTM,y=value,col=variable))+
# # geom_point(alpha=0.1,shape=22)+
# geom_smooth(method="lm",aes(col=variable))+
# geom_point(data=noaa_loc,aes(x=ELEVATION,y=ANN,col=par))+
# geom_point(data=sel_loc,aes(x=Elev,y=value,col=variable),size=5,
# shape=22,
# fill="black")
}
plot(era5_land_inter_ras)
2.4.3 Elevation Gradients
era5_land_t2m %>%
dplyr::filter(complete.cases(.)) %>%
reshape2::melt(id.var=c("x","y","z")) %>%
ggplot(data=,aes(x=z,y=value,col=variable))+
geom_point(alpha=0.5,shape=22,col="black")+
geom_smooth(data=review_loc %>% dplyr::filter(source=="noaa"),
aes(x=Elev,y=value,col=variable),
method="lm",se=F,col="red")+
geom_smooth(method="lm",col="black",se=F)+
# geom_point(data=review_loc,aes(x=ELEVATION,y=ANN,col=par))+
geom_point(data=review_loc,aes(x=Elev,y=value,
col=source),size=5)+
geom_label_repel(data=review_loc,aes(x=Elev,y=value,
label=Name %>% str_wrap(12)),alpha=0.8,
col="black",
size=2.5)+
facet_wrap(~variable,scale="free_y")+
stat_poly_eq(aes(label = paste("ERA5:",..eq.label.., ..rr.label.., sep = "~~~")),
formula = y ~ poly(x, 1, raw = TRUE),
label.x.npc = 0.1,label.y.npc = 0.1,col="black")+
stat_poly_eq(data=review_loc %>% dplyr::filter(source=="noaa"),aes(x=Elev,y=value,
label = paste("NOAA:",
..eq.label.., ..rr.label.., sep = "~~~")),
formula = y ~ poly(x, 1, raw = TRUE),
label.x.npc = 0.1,label.y.npc = 0.05,col="red"
)+
theme_light()+
theme(legend.position = "bottom")+
scale_x_continuous(breaks=seq(0,5000,200))+
labs(x="Elevation [masl]",y="Mean Annual Air Temperature [degC]",col="air temperature:",
shape="source:",
caption="Note: NOAA records from 2005 to 2022\nreviewed station records from 2012 to 2021\nERA5-Land records from 2000 to 2022",
title="Mean annual air temperature gradients",
subtitle="black : ERA5-Land \nred : NOAA DB\nblue : points of interest")+
scale_color_manual(values=c("red","blue"))
ggsave(filename = here::here("outputs","mean annual precipitation.png"),
width = 10,height = 10,units = "in",dpi = 200)2.4.4 Figures
era5_land_inter_tavg<-raster(era5_land_inter_ras$avg) %>%
rasterToPoints() %>%
data.frame()
era5_land_t2m_gp_HD<-tidyterra::fortify(era5_land_inter_ras) %>%
dplyr::filter(complete.cases(.)) %>%
# dplyr::rename(z="topo_SRTM") %>%
reshape2::melt(id.vars=c("x","y","z"))
era5_land_t2m_gp<-era5_land_t2m %>%
reshape2::melt(id.vars=c("x","y","z"))
ggmap(base6) +
# geom_tile(data=era5_land_t2m_gp_HD,aes(x=x,y=y,fill=value),alpha=0.8)+
geom_tile(data=era5_land_t2m_gp,aes(x=x,y=y,fill=value),alpha=0.8)+
geom_point(data = review_loc,
mapping = aes(x = Lon,y = Lat,fill = value,
shape=source),size = 5,
alpha=0.8) +
geom_text_repel(data = review_loc,
mapping=aes(x = Lon,y = Lat,label = Name)) +
theme_bw() +
labs(x="Longitude (X-Coords)",y="Latitude (Y-Coords)",
fill = "mean annual\nair temperature (degC)",
caption="Note: NOAA records from 2005 to 2022\nreviewed station records from 2012 to 2021\nERA5-Land records from 2000 to 2022",
title="Mean annual temperature",
subtitle="Background: Era5-Land / Point: NOAA and point of interest"
)+
scale_fill_stepsn(n.breaks=8,#breaks=c(0,10,20,50,100,200,500,1000,1500),
colours = scales::brewer_pal(palette = "Spectral")(5))+
coord_quickmap()+
scale_y_continuous(limits=c(-25.8,-24.5))+
facet_grid(variable~.)+
theme(legend.position="bottom")+
ggsn::scalebar(x.min = 29,x.max=31,y.min=-25.65,y.max=-24.8,dist = 25,
dist_unit = "km",transform = T,model = "WGS84",height = 0.05)+
scale_shape_manual(values=c(21,22))
ggsave(filename = here::here("outputs","mean annual air temperature - geographical distribution.png"),
width = 8,height = 12,units = "in",dpi = 200)
ggmap(base5) +
geom_tile(data=era5_land_t2m_gp_HD %>%
dplyr::filter(x>=29.8,x<=30.4,
y>=-25.2,y<=-24.6,
variable=="avg"),
aes(x=x,y=y,fill=value),alpha=0.8)+
# geom_tile(data=era5_land_t2m_gp,aes(x=x,y=y,fill=value),alpha=0.8)+
geom_point(data = review_loc %>% dplyr::filter(variable=="avg"),
mapping = aes(x = Lon,y = Lat,fill = value,
shape=source),shape = 21,size = 5,
alpha=0.8) +
geom_text(data = review_loc,
mapping=aes(x = Lon,y = Lat,label = Name),
nudge_x = 0.02,nudge_y = 0.02) +
theme_bw() +
labs(x="Longitude (X-Coords)",y="Latitude (Y-Coords)",
fill = "mean annual\nair temperature (degC)",
caption="Note: NOAA records from 2005 to 2022\nreviewed station records from 2012 to 2021\nERA5-Land records from 2000 to 2022",
title="Mean annual air temperature based on ERA5-Land",
subtitle="values downscaled with machine learning model")+
scale_fill_stepsn(n.breaks=10,#breaks=c(0,10,20,50,100,200,500,1000,1500),
colours = scales::brewer_pal(palette = "Spectral")(5))+
coord_quickmap()+
scale_y_continuous(limits=c(-25.2,-24.6))+
scale_x_continuous(limits=c(29.8,30.4))+
theme(legend.position="bottom")+
ggsn::scalebar(x.min = 29.8,x.max=30.4,y.min=-25.15,y.max=-24.6,dist = 5,
dist_unit = "km",transform = T,model = "WGS84")
ggsave(filename = here::here("outputs","mean annual air temperature - geographical distribution_HD_AVG.png"),
width = 12,height = 12,units = "in",dpi = 200)3 Precipitation work
# Daily Rainfall Dataset Import
pr <- read.xlsx(here::here("data","csv","MARTENSHOOP (POL)_Merged_RAW.xlsx"))
pr$Datetime <- convertToDate(pr$Datetime)
pr <- pr %>%
filter(Datetime >= "1906-01-01")
pr$Month <- month(pr$Datetime)
pr$Year <- year(pr$Datetime)
pr <- pr %>%
relocate(Datetime, Month, Year);
# Monthly Precipitation
prM <- pr %>%
group_by(Month, Year) %>%
summarize(MARTENSHOOP = sum(MARTENSHOOP),
CHIRPS = sum(CHIRPS),
ERA5 = sum(ERA5),
TRMM = sum(TRMM),
GPM = sum(GPM),
MERRA2 = sum(MERRA2),
CDR = sum(CDR),
)
prM$Period <- make_date(year = prM$Year, month = prM$Month, day = 15)
prM$Period <-as.Date(prM$Period, format(prM$Period,'%Y-%m'))
prM <- prM %>% relocate(Period)
# Monthly Count of Data
prcountM <- pr
prcountM$Period <- make_date(year = prcountM$Year, month = prcountM$Month, day = 15)
prcountM$Period <-as.Date(prcountM$Period, format(prcountM$Period,'%Y-%m'))
prcountM <- prcountM %>%
group_by(Period, Month, Year) %>%
summarise(MARTENSHOOP = sum(!is.na(MARTENSHOOP)),
CHIRPS = sum(!is.na(CHIRPS)),
ERA5 = sum(!is.na(ERA5)),
TRMM = sum(!is.na(TRMM)),
GPM = sum(!is.na(GPM)),
MERRA2 = sum(!is.na(MERRA2)),
CDR = sum(!is.na(CDR )),
)
prcountM
prcountM$DaysInMonth <- days_in_month(prcountM$Period)
prcountM$MARTENSHOOP_Pct <- round((prcountM$MARTENSHOOP / prcountM$DaysInMonth) * 100, digits = 2)
prcountM$CHIRPS_Pct <- round((prcountM$CHIRPS / prcountM$DaysInMonth) * 100, digits = 2)
prcountM$ERA5_Pct <- round((prcountM$ERA5 / prcountM$DaysInMonth) * 100, digits = 2)
prcountM$TRMM_Pct <- round((prcountM$TRMM / prcountM$DaysInMonth) * 100, digits = 2)
prcountM$GPM_Pct <- round((prcountM$GPM / prcountM$DaysInMonth) * 100, digits = 2)
prcountM$MERRA2_Pct <- round((prcountM$MERRA2 / prcountM$DaysInMonth) * 100, digits = 2)
prcountM$CDR_Pct <- round((prcountM$CDR / prcountM$DaysInMonth) * 100, digits = 2)
#prcountM <- subset(prcountM, select = -c(MARTENSHOOP, CHIRPS, ERA5, TRMM, GPM, MERRA2, CDR))
#colnames(prcountM) <- c("Period", "Month", "Year", "DaysInMonth", "MARTENSHOOP", "CHIRPS", "ERA5", "TRMM", "GPM", "MERRA2", "CDR")
prcountM
# Annual Count of Data
prcountY <- prcountM
prcountY <- prcountY %>%
group_by(Year) %>%
summarise(MARTENSHOOP = sum(MARTENSHOOP),
CHIRPS = sum(CHIRPS),
ERA5 = sum(ERA5),
TRMM = sum(TRMM),
GPM = sum(GPM),
MERRA2 = sum(MERRA2),
CDR = sum(CDR),
DaysInYear = sum(DaysInMonth)
)
prcountY$MARTENSHOOP_Pct <- round((prcountY$MARTENSHOOP / prcountY$DaysInYear) * 100, digits = 2)
prcountY$CHIRPS_Pct <- round((prcountY$CHIRPS / prcountY$DaysInYear) * 100, digits = 2)
prcountY$ERA5_Pct <- round((prcountY$ERA5 / prcountY$DaysInYear) * 100, digits = 2)
prcountY$TRMM_Pct <- round((prcountY$TRMM / prcountY$DaysInYear) * 100, digits = 2)
prcountY$GPM_Pct <- round((prcountY$GPM / prcountY$DaysInYear) * 100, digits = 2)
prcountY$MERRA2_Pct <- round((prcountY$MERRA2 / prcountY$DaysInYear) * 100, digits = 2)
prcountY$CDR_Pct <- round((prcountY$CDR / prcountY$DaysInYear) * 100, digits = 2)
prcountY <- subset(prcountY, select = -c(MARTENSHOOP, CHIRPS, ERA5, TRMM, GPM, MERRA2, CDR))
colnames(prcountY) <- c("Year", "DaysInYear", "MARTENSHOOP", "CHIRPS", "ERA5", "TRMM", "GPM", "MERRA2", "CDR")
prcountY
# Annual Long Data for Vis
prcountY_L <- prcountY
prcountY_L <- pivot_longer(prcountY_L,
cols = 3:9,
names_to = "Station",
values_to = "PCT")
prcountY_L <- prcountY_L %>% arrange(Station, Year, .by_group6 = TRUE)
prcountY_L # .height=3.93701, fig.width=6.29921
# Annual Data Coverage Graph
p1 <- ggplot(data = prcountY_L, aes(x = Year, y = Station, fill = PCT)) +
geom_tile(color = "black",
size = 0.005,
alpha = 0.9) +
#scale_fill_gradient(low = "white", high = "blue") +
scale_fill_distiller(type = "seq",
#palette = "Greys",
palette = "YlGnBu",
direction = 1,
na.value = "Gray") +
theme_bw() +
theme(plot.title = element_text(size = 13),
axis.text = element_text(size = 10),
axis.title = element_text(size = 11))+
labs(#title = "Meteorological Dataset Annual Data Coverage (%)",
#title = "OMAS OKSUT Project",
#subtitle = "Local Meteorological Stations' Annual Data Coverage (%)",
fill = "Percentage (%)",
y = NULL) +
guides(fill = guide_colourbar(barwidth = 0.75,
barheight = 7))+
coord_fixed(ratio = 10,
xlim = c(1900,2023),
ylim = NULL,
expand = FALSE,
clip = "on") +
scale_y_discrete(limits = c("CDR",
"MERRA2",
"GPM",
"TRMM",
"ERA5",
"CHIRPS",
"MARTENSHOOP"))
p1
#pr2 <- read_excel("MARTENSHOOP (POL)_Merged_RAW.xlsx")
#pr2 <- pr2 %>%
# filter(Datetime >= "1906-01-01")
pr2 <- pr
#pr2 <- select(pr2, c("Month","Year","MARTENSHOOP", "CHIRPS", "ERA5", "TRMM", "GPM","MERRA2", "CDR"))
pr2 <- dplyr::select(pr2, c("Month","Year","MARTENSHOOP", "CHIRPS", "ERA5","MERRA2", "CDR"))
pr2 <- pr2 %>%
drop_na()
pr2 <- pr2 %>%
group_by(Month, Year) %>%
summarise(MARTENSHOOP = sum(MARTENSHOOP),
CHIRPS = sum(CHIRPS),
ERA5 = sum(ERA5),
MERRA2 = sum(MERRA2),
CDR = sum(CDR),
) %>%
arrange(Year)
pr2 <- subset(pr2, select = -c(Month, Year)) # Delete Month and Year Cols.
pr2
p2 <- GGally::ggpairs(pr2,
#title = "Develi OMGI - Gridded Models Precipitation Correlation",
#columns = 1:4,
progress = FALSE, # To suppress ggpairs progress notes.
proportions = "auto",
upper = list(continuous = GGally::wrap("cor", size = 3, color = "red",
display_grid = FALSE)),
#diag = list(continuous = "densityDiag"),
lower = list(continuous = GGally::wrap("smooth", alpha = 0.4,
size = 1, color = "blue"))
) +
theme_bw(base_size = 8.75) +
theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1)) +
theme(strip.text.y.right = element_text(angle = 0, vjust = 0.5, hjust=0))
p2
ggsave(here::here("outputs","MARTENSHOOP_Pr_ScatterPlot_Matrix.png"),
width = 16,
height = 11,
units = c("cm"),
plot = p2,
dpi = 600)saws <- read_excel(here::here("data","csv","TwoRivers_CHIRPS_Raster_Elv_Precipitation.xlsx"),
sheet = "saws" )
CHIRPS2000 <- read_excel(here::here("data","csv","TwoRivers_CHIRPS_Raster_Elv_Precipitation.xlsx"),
sheet = "df_CHIRPS2000_DEM")
CHIRPS <- read_excel(here::here("data","csv","TwoRivers_CHIRPS_Raster_Elv_Precipitation.xlsx"),
sheet = "df_CHIRPS_DEM")
df_CHIRPS <- read_excel(here::here("data","csv","TwoRivers_CHIRPS_Raster_Elv_Precipitation.xlsx"),
sheet = "df_CHIRPS")
df_CHIRPS_1 <- df_CHIRPS %>%
dplyr::filter(Type != "CHIRPS" & Type != "SAWS SRTM (1981-2000)")
df_CHIRPS_SRTM <- df_CHIRPS %>%
dplyr::filter(Type != "CHIRPS" & Type != "SAWS (1981-2000)")
saws; CHIRPS2000; CHIRPS; df_CHIRPS; df_CHIRPS_1; df_CHIRPS_SRTM
p3 <- ggplot(data = df_CHIRPS_1, mapping = aes(y = MAP, x = Elevation, color = Type)) +
geom_point(#shape = 1,
size = 2,
#color = "black",
#fill = Type,
alpha = 1) +
#(method = lm, fullrange = TRUE) +
scale_colour_brewer(palette = "Set1") +
stat_poly_line() +
stat_poly_eq(mapping = use_label(c("eq", "R2"))) +
geom_segment(aes(x = 1000, xend = 1000, y = -Inf, yend = +Inf), linetype = 2, color = "black") +
geom_segment(aes(x = 970, xend = 970, y = -Inf, yend = +Inf), linetype = 2, color = "black") +
theme_bw() +
theme(legend.position = "bottom") +
labs(x = "Elevation (mamsl)",
y = "MAP (mm)",
color = "Dataset") +
geom_text(aes(label="De Grooteboom TSF",
x=1020,
y = 750,
angle = 90),size=4,color = "black") +
geom_text(aes(label="Old TSF",
x=940,
y = 750,
angle = 90),size=4,color = "black") +
geom_text_repel(data = df_CHIRPS_1,
aes(x = Elevation,
y = MAP,
label = NAME,
), size = 3)
p4 <- ggplot(data = df_CHIRPS_SRTM, mapping = aes(y = MAP, x = Elevation, color = Type)) +
geom_point(#shape = 1,
size = 2,
#color = "black",
#fill = Type,
alpha = 1) +
#(method = lm, fullrange = TRUE) +
scale_colour_brewer(palette = "Set1") +
stat_poly_line() +
stat_poly_eq(mapping = use_label(c("eq", "R2"))) +
geom_segment(aes(x = 1000, xend = 1000, y = -Inf, yend = +Inf), linetype = 2, color = "black") +
geom_segment(aes(x = 970, xend = 970, y = -Inf, yend = +Inf), linetype = 2, color = "black") +
theme_bw() +
theme(legend.position = "bottom") +
labs(x = "Elevation (mamsl)",
y = "MAP (mm)",
color = "Dataset") +
geom_text(aes(label="De Grooteboom TSF",
x=1020,
y = 750,
angle = 90),size=4,color = "black") +
geom_text(aes(label="Old TSF",
x=940,
y = 750,
angle = 90),size=4,color = "black") +
geom_text_repel(data = df_CHIRPS_SRTM,
aes(x = Elevation,
y = MAP,
label = NAME,
), size = 3)
p3; p4
GPM <- read_excel(here::here("data","csv","TwoRivers_GPM_Raster_Elv_Precipitation.xlsx"),
sheet = "df_GPM_DEM")
df_GPM <- read_excel(here::here("data","csv","TwoRivers_GPM_Raster_Elv_Precipitation.xlsx"),
sheet = "df_GPM")
df_GPM_1 <- df_GPM %>%
dplyr::filter(Type != "SAWS SRTM (1981-2000)")
df_GPM_SRTM <- df_GPM %>%
dplyr::filter(Type != "SAWS (1981-2000)")
saws; GPM; df_GPM; df_GPM_1; df_GPM_SRTM
p5 <- ggplot(data = df_GPM_1, mapping = aes(y = MAP, x = Elevation, color = Type)) +
geom_point(#shape = 1,
size = 2,
#color = "black",
#fill = Type,
alpha = 1) +
#(method = lm, fullrange = TRUE) +
scale_colour_brewer(palette = "Set1") +
stat_poly_line() +
stat_poly_eq(mapping = use_label(c("eq", "R2"))) +
geom_segment(aes(x = 1000, xend = 1000, y = -Inf, yend = +Inf), linetype = 2, color = "black") +
geom_segment(aes(x = 970, xend = 970, y = -Inf, yend = +Inf), linetype = 2, color = "black") +
theme_bw() +
theme(legend.position = "bottom") +
labs(x = "Elevation (mamsl)",
y = "MAP (mm)",
color = "Dataset") +
geom_text(aes(label="De Grooteboom TSF",
x=1020,
y = 650,
angle = 90),size=4,color = "black") +
geom_text(aes(label="Old TSF",
x=940,
y = 650,
angle = 90),size=4,color = "black") +
geom_text_repel(data = df_GPM_1,
aes(x = Elevation,
y = MAP,
label = NAME,
), size = 3)
p6 <- ggplot(data = df_GPM_SRTM, mapping = aes(y = MAP, x = Elevation, color = Type)) +
geom_point(#shape = 1,
size = 2,
#color = "black",
#fill = Type,
alpha = 1) +
#(method = lm, fullrange = TRUE) +
scale_colour_brewer(palette = "Set1") +
stat_poly_line() +
stat_poly_eq(mapping = use_label(c("eq", "R2"))) +
geom_segment(aes(x = 1000, xend = 1000, y = -Inf, yend = +Inf), linetype = 2, color = "black") +
geom_segment(aes(x = 970, xend = 970, y = -Inf, yend = +Inf), linetype = 2, color = "black") +
theme_bw() +
theme(legend.position = "bottom") +
labs(x = "Elevation (mamsl)",
y = "MAP (mm)",
color = "Dataset") +
geom_text(aes(label="De Grooteboom TSF",
x=1020,
y = 650,
angle = 90),size=4,color = "black") +
geom_text(aes(label="Old TSF",
x=940,
y = 650,
angle = 90),size=4,color = "black") +
geom_text_repel(data = df_GPM_SRTM,
aes(x = Elevation,
y = MAP,
label = NAME,
), size = 3)
p5; p6
p7 <- p3 / p5
p7
ggsave(here::here("outputs","MAP-Elevation_Relationship_V2.png"),
width = 16,
height = 19,
units = c("cm"),
plot = p7,
dpi = 600)4 Output
file.image<-here::here("data","rds",paste0("Backup Meteorological review_",as.Date(now()),".rdata"))
save.image(file = file.image)