Lab3_ET
John Gardner
2022-09-22
Load and organize data from DAYMET
First review the DAYMET data, https://daymet.ornl.gov/overview For reference on the R package and how to use it, https://cran.r-project.org/web/packages/daymetr/daymetr.pdf and https://cran.r-project.org/web/packages/daymetr/vignettes/daymetr-vignette.html For more info on SPEI, https://spei.csic.es/home.html
### Directions: look up the latitude and longitude of Phoenix, AZ and Miami, FL. Load the DAYMET data for these two ### city coordinates from 1990-2020 using the download_daymet function (see R package links above). If you have your own data of monthly precipitation, monthly mean, min, and max temperature for another site, feel free to use that instead. FINISH this code
phx <- download_daymet(
lat = 33.4484,
lon = -112.074,
start =1990,
end = 2020,
internal = TRUE,
simplify = TRUE) ## Downloading DAYMET data for: Daymet at 33.4484/-112.074 latitude/longitude !## Done !mia <- download_daymet(
lat = 25.7617,
lon = -80.1918,
start = 1990,
end = 2020,
internal = TRUE,
simplify = TRUE) ## Downloading DAYMET data for: Daymet at 25.7617/-80.1918 latitude/longitude !
##
## Done !# Take a look at the data structure
str(phx)## tibble [79,205 × 9] (S3: tbl_df/tbl/data.frame)
## $ site : chr [1:79205] "Daymet" "Daymet" "Daymet" "Daymet" ...
## $ tile : num [1:79205] 11014 11014 11014 11014 11014 ...
## $ latitude : num [1:79205] 33.4 33.4 33.4 33.4 33.4 ...
## $ longitude : num [1:79205] -112 -112 -112 -112 -112 ...
## $ altitude : num [1:79205] 334 334 334 334 334 334 334 334 334 334 ...
## $ year : int [1:79205] 1990 1990 1990 1990 1990 1990 1990 1990 1990 1990 ...
## $ yday : int [1:79205] 1 2 3 4 5 6 7 8 9 10 ...
## $ measurement: chr [1:79205] "dayl..s." "dayl..s." "dayl..s." "dayl..s." ...
## $ value : num [1:79205] 35386 35417 35451 35488 35527 ...# look at what the unique variables are. NOTE, this data is in long format.
unique(phx$measurement)## [1] "dayl..s." "prcp..mm.day." "srad..W.m.2." "swe..kg.m.2."
## [5] "tmax..deg.c." "tmin..deg.c." "vp..Pa."### First for the Phoenix data,
#1) use the pivot_wider function to convert to wide format (HINT: using measurement and value columns),
#2) calculate the mean temperature using tmax and tmin,
#3) do date conversions (code is provided, but try to understand how code works),
#4) convert shortwave radiation flux to total daily radiation (HINT: Look at the DAYMET link at top of script, go to "Description tab" then "model outputs". There is conversion equation shown in the table).
#5) if you want, rename the variables since the names from DAYMET are weird and long. (HINT: this is easiest if you rename() after pivoting to wide format)
#FINISH this code
phx_df <- phx %>%
pivot_wider(names_from = measurement, values_from = value) %>%
rename(dayl_s = dayl..s., prcp_mmday = prcp..mm.day., srad_wm2 = srad..W.m.2., swe_kgm2 = swe..kg.m.2., tmax_c = tmax..deg.c., tmin_c = tmin..deg.c., vp_pa = vp..Pa.) %>%
mutate(phx_mean=(tmax_c+tmin_c)/2) %>%
mutate(srad_mjm2=dayl_s/1000000) %>%
# Code for converting dates to year month and still make the date and yearmonth column plot-able, if you are using tidyverse/pipes style of coding
mutate(origin = as.Date(paste0(year, "-01-01"),tz = "UTC") -1 ,
date = as.Date(yday, origin = origin, tz="UTC")) %>%
mutate(yearmonth = as.Date(format(date, format = "%Y-%m-01")))
### Repeat for the Miami data, you can copy and paste code.
mia_df <- mia %>%
pivot_wider(names_from = measurement, values_from = value) %>%
rename(dayl_s = dayl..s., prcp_mmday = prcp..mm.day., srad_wm2 = srad..W.m.2., swe_kgm2 = swe..kg.m.2., tmax_c = tmax..deg.c., tmin_c = tmin..deg.c., vp_pa = vp..Pa.) %>%
mutate(mia_mean=(tmax_c+tmin_c)/2) %>%
mutate(srad_mjm2=dayl_s/1000000) %>%
# Code for converting dates to year month and still make the date and yearmonth column plot-able, if you are using tidyverse/pipes style of coding
mutate(origin = as.Date(paste0(year, "-01-01"),tz = "UTC") -1 ,
date = as.Date(yday, origin = origin, tz="UTC")) %>%
mutate(yearmonth = as.Date(format(date, format = "%Y-%m-01")))
### Our data is daily, but we only need monthly data to compute Potential Evapotranspiration (PET). For each city, compute monthly statistics by yearmonth. Make sure your final summary has columns for latitude (Hint: latitude is the same so think about how to "index" the first value in a column), monthly SUM of precipitation (mm), monthly MEAN max temperature, monthly MEAN min temperature, monthly MEAN mean temperature, monthly MEAN total daily solar radiation. FINISH this code.
phx_month <- phx_df %>%
group_by(yearmonth) %>%
summarise(latitude=latitude[1],month_sum=sum(phx_mean),monthmt_max=max(tmax_c),monthmt_mean=mean(phx_mean),srad_mean=mean(srad_wm2),monthmt_min=min(tmin_c),monthpre=sum(prcp_mmday)) %>%
ungroup()
mia_month <- mia_df %>%
group_by(yearmonth) %>%
summarise(latitude=latitude[1],month_sum=sum(mia_mean),monthmt_max=max(tmax_c),monthmt_mean=mean(mia_mean),srad_mean=mean(srad_wm2),monthmt_min=min(tmin_c),monthpre=sum(prcp_mmday)) %>%
ungroup()
###FINISH code: save the monthly data and the daily dataframes you made for each city so you do not have to download them again. Once you have succesfully saved the files, comment out this code
write_csv(mia_month, path ="mia.csv")## Warning: The `path` argument of `write_csv()` is deprecated as of readr 1.4.0.
## Please use the `file` argument instead.write_csv(phx_month, path ="phx.csv")
write_csv(phx_df, path ="phx.csv")
write_csv(mia_df, path ="mia.csv")Look at the data
Question 1: Why does Phoenix have so much more incoming solar radiation compared to Miami even though Miami is further south? What impact will this have on PET?
TYPE ANSWER HERE: Miami has more cloud cover. this could reduce PET in Miami compared to clear days in Miami
### FIGURE 1: plot a timeseries of daily mean temperature as a geom_line() with both cities on the same plot using a different color distinguish the cities. FINSIH this code
ggplot() +
geom_line(data=mia_df,aes(x=date, y=mia_mean),color="green")+
geom_line(data=phx_df,aes(x=date, y=phx_mean),color="blue") +
theme_bw()
### Figure 2: plot a timeseries of monthly solar radiation as a geom_line() with both cities on the same plot using a different color distinguish the cities. FINISH this code
ggplot()+
geom_line(data=mia_month,aes(x=yearmonth, y=srad_mean),color="green")+
geom_line(data=phx_month,aes(x=yearmonth, y=srad_mean),color="blue")+
theme_bw()
PET calculation
QUESTION 2: Since we are not using external (top of atmosphere) solar radiation to drive PET in Thornthwaite and Hargreaves calculations here, how do these methods figure out external radiation? (Hint: look at the full equations used)
TYPE ANSWER HERE: using latitude and geometry between sun and earth to get estimates of potential solar radiation
COMPARE PET
QUESTION 3: Estimating from the graphs, how similar are the three PET methods for Miami? And for Phoenix? IF the PET values are more similar for one of the cities, why do you think that is?
TYPE ANSWER HERE: Thornwaithe always has larger range, higher highs and lower lows. PET calculations are more similar for PHX. Less precip, its easier to estimate and the precip correction for Hargreaves method will not change PET that much
###FIGURE 3: For Miami, plot the monthly precipitation and PET from all three methods on the same plot using geom_line and different colors and or linetypes for each of the four time series. The units of precip and PET should already be in units of mm.
ggplot() +
geom_line(data=mia_month,aes(x=monthpre, y=mia_PET),color="green") +
geom_line(data=mia_month,aes(x=monthpre, y=mia_PET2),color="red") +
geom_line(data=mia_month,aes(x=monthpre, y=mia_PET3),color="blue") +
theme_bw() ## Don't know how to automatically pick scale for object of type ts. Defaulting to continuous.
### FIGURE 4: repeat for phoenix
ggplot() +
geom_line(data=phx_month,aes(x=monthpre, y=phx_PET),color="green") +
geom_line(data=phx_month,aes(x=monthpre, y=phx_PET2),color="red") +
geom_line(data=phx_month,aes(x=monthpre, y=phx_PET3),color="blue") +
theme_bw()## Don't know how to automatically pick scale for object of type ts. Defaulting to continuous.
Calculating drought index: SPEI
QUESTION 4: Looking at the 2 plots for Phoenix, describe the duration and intensity of drought (SPEI). Do you think there is a trend in drought frequency, duration, or intensity? Do the two PET methods provide similar values and interpretation of SPEI?
TYPE ANSWER HERE: long intense droughts. perhaps increasing frequency, but interpretation may differ depending on PET method used
QUESTION 5: Now looking at the plots for Miami, it is perhaps surprising that there are droughts here. What do droughts and wet periods look like in Miami (according to plot), and how can droughts happen in Miami FL where there is so much rain?
TYPE ANSWER HERE: very short droughts. they can happen cause ET is so high
### Using the spei() function calculate SPEI for both cities using PET data from 2 of the 3 methods:Thornthwaite and Hargreaves with precipitation correction method. Play around with different scale arguments (which are in units of months) and look at the output plots. Use a final scale of 3. SPEI requires inputs of water balance (Precipitation - PET) and scale.
# to see help on this function
?spei()## starting httpd help server ... done# Phoenix with 2 methods
phx_SPEI_tho <- spei(phx_PET,3)
phx_SPEI_har <- spei(phx_PET3,3)
# FIGURE 5-6: make plots of SPEI drought index over time for Phoenix
plot(phx_SPEI_tho)
plot(phx_SPEI_har)
# Miami with 2 methods
mia_SPEI_tho <- spei(mia_PET,3)
mia_SPEI_har <- spei(mia_PET3,3)
# FIGURE 7-8: make plots of SPEI drought index over time for Miami
plot(mia_SPEI_tho)
plot(mia_SPEI_har)
BIAS check
QUESTION 6: Inspect the plots of PET from Thornthwaite vs. Hargreaves methods. How much scatter is there in Miami compared to Phoenix? Would you recommend use of these PET calculation methods for both places?
TYPE ANSWER HERE: Phoenix has more scatter in Phoenix, but both fit along the 1 to 1 line reasonably well. Maybe one of the methods is less applicable in phoenix compared to miami
### The SPEI calculations are stored in a list with other outputs of the spei() function. We have to extract the "fitted" values from this list, which are the SPEI values calculated by fitting a probability distribution to the data.
# The SPEI values can be extracted from an object like this.
mia_SPEI_har$fitted## Jan Feb Mar Apr May Jun
## 1 NA NA -1.86943728 -1.89246196 -1.64306041 0.30835119
## 2 -0.44340961 0.73824729 1.13367846 1.13415764 -0.74797737 -0.84559282
## 3 -0.39636848 -1.28245803 -0.76775283 -0.38314026 1.36264357 0.74069797
## 4 -1.46222134 -1.89260359 -1.01962715 0.10859389 0.59734381 0.98912491
## 5 0.15656332 -0.75142603 -0.72035874 -1.33926676 -0.06023619 0.03824462
## 6 -1.43049346 0.01957990 -0.32954386 -0.03202449 0.46778973 0.13679998
## 7 1.69215994 2.92084722 1.79719755 1.25345339 -0.09485046 -0.91594306
## 8 1.28722866 0.83260527 -1.18851089 -0.88187654 -0.17447684 0.38150645
## 9 -0.07229090 -1.12105743 -1.71744416 -0.89208543 0.50439316 2.01570905
## 10 -0.96062057 -0.41516801 -0.55750444 0.81449140 1.75748486 0.45162381
## 11 -0.12941884 -0.14087479 -0.96672547 -0.85170569 -0.37582392 0.43870032
## 12 0.76024118 -0.10419052 0.58580378 1.13576720 0.47962097 -0.20918199
## 13 -0.61975399 0.71115481 0.78757970 -0.89183915 -1.51056623 -2.29853353
## 14 1.31841527 0.72873244 0.76840860 0.61774584 -0.58075679 -1.21859436
## 15 -0.86309788 -0.91374879 -1.57950297 -1.46636702 -0.83569977 -0.28754267
## 16 -0.08038028 0.12393876 -0.72776693 0.16455563 1.00760952 0.47750952
## 17 -0.21672848 -0.57189165 0.38501798 0.45835771 0.01578271 -0.36106373
## 18 -0.29817442 -0.89615192 -0.65574531 -0.11121582 -0.14713376 0.00298755
## 19 0.14862822 0.73098706 0.92688079 0.30554369 1.24021717 1.02063295
## 20 0.33820528 0.61633023 1.06179613 1.59846968 -0.39274793 0.65346349
## 21 1.62249534 0.91150133 0.50091255 -1.29702428 -1.25331597 -1.18303598
## 22 1.03768198 1.16288727 0.10918343 -1.24738314 -0.89621938 -0.85625329
## 23 -0.26226510 -0.20419160 -0.10004241 -0.55274557 -1.59585688 -1.16301770
## 24 -0.58443806 0.03285088 1.10912866 0.80323513 0.83903637 -0.66605099
## 25 -0.06926755 -0.38627104 0.50766169 -0.15022851 -0.59387621 -0.72148006
## 26 0.48839239 0.68287408 0.46173739 1.76958933 1.49640972 2.19659045
## 27 -2.41890051 -2.10796407 -0.59226338 -0.66163264 0.03391672 0.71806890
## 28 -0.98721279 -0.64363979 0.09433643 0.28273181 1.68980597 1.13786764
## 29 -0.23592854 -0.45855958 0.11093392 0.19652358 -1.53080542 -1.43461063
## 30 1.60640354 0.32732777 0.25427474 0.49043041 0.42593330 -0.05190476
## 31 1.12382883 1.09487457 2.03802600 1.55976364 0.60707172 0.37035115
## Jul Aug Sep Oct Nov Dec
## 1 1.11602237 1.54855345 1.48293657 2.05886328 1.54274884 0.85634550
## 2 -0.78798404 1.14810115 0.88623634 0.41214670 -0.22417593 -0.59058753
## 3 1.39871509 0.06717009 1.08229127 0.71600524 0.52706430 0.05585617
## 4 1.10888361 1.49755678 0.63500852 0.04550962 -0.33427286 0.05286876
## 5 0.92091960 -0.15780287 -1.18393439 -1.66541480 -1.55014340 -1.27761400
## 6 0.57467156 -0.10966823 0.74562584 -0.87797501 -0.43499176 -0.27365718
## 7 -0.57295831 0.87172667 1.76229410 1.28387960 0.51311941 0.14134986
## 8 0.34966295 0.07390800 -0.87768232 -0.46462960 -0.05869331 1.07550317
## 9 2.13030814 1.90859104 1.48169912 1.21699925 -0.04072011 -0.35110892
## 10 0.90904483 0.16435197 1.75361066 0.17929725 -0.72751043 -0.85777226
## 11 -0.37937898 -0.06552428 -0.83181070 -0.46816449 -0.02526785 0.27587233
## 12 -0.69129850 0.50619697 0.03026396 0.01492414 -1.17266657 -0.61307729
## 13 -2.05926813 -1.57703425 -0.14840328 0.28607676 0.71761929 0.71664491
## 14 -1.75741491 -1.91409278 -1.65934891 -1.14212435 -0.60439303 -0.35397955
## 15 -0.26632457 -0.89176807 -1.43174291 -1.30738197 -0.95055097 -0.14016626
## 16 -0.01986662 -0.99490915 0.06351627 1.29759810 1.23034733 1.02945738
## 17 -1.04870978 -1.00146697 -1.44068115 0.20548818 0.94987860 0.85810438
## 18 -0.35214862 0.09060463 0.44056013 -0.52567267 -1.02353036 -1.30393109
## 19 0.46693062 -0.81328259 -0.50111137 0.70833383 0.89351895 0.62723938
## 20 0.02276689 1.37878900 0.53641836 1.49571222 1.77939273 1.95678173
## 21 -0.42545067 0.11231188 -0.77985720 -1.11887537 -0.93425771 0.40151731
## 22 0.11001849 -0.20962917 -0.23777979 -0.91421944 -0.46012552 -0.96408799
## 23 -1.41986560 -1.20330950 -0.99900141 0.52503964 0.94694612 1.15693336
## 24 -0.51937443 -0.60518665 -0.03323821 -0.17023799 -0.22571826 -0.78412324
## 25 -0.49038971 -0.04914836 0.07199556 0.60825139 0.80435852 1.11539601
## 26 1.56780091 1.16233014 0.20275188 -0.71026530 -1.05755838 -1.87906692
## 27 0.72707185 -0.25563829 -0.19982017 0.01553791 0.24839183 -0.71904046
## 28 0.64971451 -1.64960062 -1.21280079 -0.61757369 -0.66766349 0.19966774
## 29 -1.33304302 1.16426777 0.66394386 1.11094404 1.44022467 1.78295224
## 30 0.09171917 -0.33578726 0.95557888 -0.19636815 1.01680333 -0.38491353
## 31 -0.05348971 0.28218639 -0.92451718 -2.13011207 -2.27319976 -2.01348745### Make a new dataframe for each city, which has 2 columns including the SPEI index calculated from Thornthwaite (eg. mia_SPEI_tho) and Hargreaves method estimated in the previous section. (HINT: use cbind(), extract the "fitted" values, and you may have to convert the SPEI values into a dataframe. You may also find it easier to rename the column names)
mia_spei <- cbind(PET_tho = as.data.frame(mia_SPEI_tho$fitted) %>% rename (tho=1),
PET_har = as.data.frame(mia_SPEI_har$fitted) %>% rename (har=1))
phx_spei <- cbind(PET_tho = as.data.frame(phx_SPEI_tho$fitted) %>% rename (tho=1),
PET_har = as.data.frame(phx_SPEI_har$fitted) %>% rename (har=1))
# FIGURE 9: plot the Thornthwaite vs. Hargreaves SPEI values for Miami also plotting the 1 to 1 line using geom_abline()
ggplot(mia_spei) +
geom_point(aes(x=tho, har)) +
geom_abline(slope=1, intercept=0)## Warning: Removed 2 rows containing missing values (geom_point).
# FIGURE 10: Repeat for Phoenix
ggplot(phx_spei) +
geom_point(aes(x=tho, har)) +
geom_abline(slope=1, intercept=0)## Warning: Removed 2 rows containing missing values (geom_point).
Note that the echo = FALSE parameter was added to the
code chunk to prevent printing of the R code that generated the
plot.