forked from jagephart/FishPrint
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathwrangle_evap.R
127 lines (96 loc) · 5.34 KB
/
wrangle_evap.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
# Author: Kelvin Gorospe
# Download NOAA global evaporation dataset for evaporative water loss model
rm(list=ls())
library(RCurl) # for getURL
library(tidyverse)
library(stars) # newer package for working with raster data; better integration with sf package
library(sf)
library(raster) # needed for rotate to reset climate data from 0 to 360 to -180 to 180
library(rnaturalearth)
# GET FILES from NOAA FTP:
# https://ftp.cpc.ncep.noaa.gov/wd51yf/global_monthly
# Download Climatology Data:
climatology_evap <- "/Volumes/jgephart/BFA Environment 2/Data/Evaporation/Monthly Climatology 1981-2010"
url <- "ftp://ftp.cpc.ncep.noaa.gov/wd51yf/global_monthly/GeoTIFF/clim/"
# Website: https://ftp.cpc.ncep.noaa.gov/wd51yf/global_monthly/GeoTIFF/monthly/
filenames <- getURL(url, ftp.use.epsv = FALSE, dirlistonly = TRUE)
filenames <- strsplit(filenames, "\n")
filenames <- unlist(filenames)
for (filename in filenames) {
download.file(url = paste(url, filename, sep = ""), destfile = file.path(climatology_evap, filename))
}
# Process NOAA global evaporation dataset
####################################################### READ IN monthly climatology files
outdir <- "/Volumes/jgephart/BFA Environment 2/Outputs"
outdir_evap <- "/Volumes/jgephart/BFA Environment 2/Outputs/Evaporation Outputs"
# Use rotate() to transform coordinate system to have standard x coordinates: xmin = -180 to xmax = 180
# Convert longitude from 0 - 360 to -180 to 180 (Standard transformation for climate data)
# https://stackoverflow.com/questions/25730625/how-to-convert-longitude-from-0-360-to-180-180
clim_files <- list.files(file.path(climatology_evap))
clim_raster <- lapply(clim_files, function(i){raster(file.path(climatology_evap, i))}) # read in with raster() function so that rotate() function works
clim_rotate <- lapply(clim_raster, function(i){rotate(i)})
clim_stars_list <- lapply(clim_rotate, function(i){st_as_stars(i)})
# Generate expression for c() based on length(clim_stars) and evaluate with eval()
stars_objects <- paste("clim_stars_list[[", 1:length(clim_stars_list), "]]", sep = "", collapse = ", ")
clim_stars <- eval(parse(text = paste("c(", stars_objects, ", along = list(month = clim_files))", sep = "")))
# Plot with ggplot:
ggplot() +
geom_stars(data = clim_stars) +
facet_wrap(~month)
ggsave(file.path(outdir_evap, "clim_panels.png"), width = 6, height = 4, unit = "in")
# Plot each month separately
clim_dat_list <- read_stars(file.path(climatology_evap, clim_files))
for (i in 1:length(clim_dat_list)){
ggplot() +
geom_stars(data = clim_dat_list[i])
ggsave(file.path(outdir_evap, paste("clim_list_", i, ".png", sep = "")), width = 6, height = 4, unit = "in")
}
####################################################### Process raster data
# <calculate mean across climatology months>
# Compute mean for each pixel (x,y) across all months
# use sum(is.na(clim_dat_list[[i]])) to check for NA's; no NAs found but keep na.rm = TRUE anyway
clim_mean_stars <- st_apply(clim_stars, c("x", "y"), mean, na.rm = TRUE)
# Plot with ggplot:
ggplot() +
geom_stars(data = clim_mean_stars)
ggsave(file.path(outdir_evap, "clim_mean.png"), width = 6, height = 4, unit = "in")
####################################################### Add country polygons
# Get countries from naturalearth package and set to same crs as clim_mean_stars
# datum = WGS84 i.e., standard lat/long and GPS coordinate system
# projection (instructions on how to transform datum onto a flat surface)
world_sf<- ne_countries(returnclass='sf') %>%
dplyr::select(admin) %>%
#st_transform("+proj=moll")
st_transform(st_crs(clim_mean_stars))
st_crs(world_sf) == st_crs(clim_mean_stars)
ggplot() +
geom_stars(data = clim_mean_stars) +
#geom_stars(data = clim_mean_stars, mapping = aes(x = x, y = y, fill = mean)) +
geom_sf(data = world_sf, alpha = 0.1, color = "black")
ggsave(file.path(outdir_evap, "clim_mean_world_map.png"), width = 6, height = 4, unit = "in")
# CROP raster and replot
clim_mean_stars <- st_crop(clim_mean_stars, world_sf)
ggplot() +
geom_stars(data = clim_mean_stars) +
geom_sf(data = world_sf, alpha = 0.1, color = "black")
ggsave(file.path(outdir_evap, "clim_mean_world_map_2.png"), width = 6, height = 4, unit = "in")
## although coordinates are longitude/latitude, st_intersects assumes that they are planar
# See: https://www.r-spatial.org/r/2020/06/17/s2.html (The Earth is no longer flat in r-spatial)
####################################################### Summarize data at the country-level
clim_mean_sf <- st_as_sf(clim_mean_stars)
clim_world_sf <- st_intersection(clim_mean_sf, world_sf)
clim_world_sf <- set_units(clim_world_sf$mean, mm)
# Calculate ANNUAL mean within each country (multiply by 12)
clim_by_country <- clim_world_sf %>%
group_by(admin) %>%
summarise(mean_evap_mm = mean(mean, na.rm = TRUE) * 12) %>% # units are in "mm"/ year
ungroup()
ggplot(clim_by_country) +
geom_sf(mapping = aes(fill = mean_evap_mm)) +
labs(fill = "Country-level evaporation (mm / year)")
ggsave(file.path(outdir_evap, "clim_summarise_by_country.png"), width = 6, height = 4, unit = "in")
write.csv(clim_by_country %>% st_set_geometry(NULL), file = file.path(outdir, "clim_summarise_by_country.csv"), quote = FALSE)
# BAR GRAPH:
#ggplot(clim_by_country) +
# geom_col(mapping = aes(x = reorder(admin, desc(country_level)), y = country_level))
#country_level_clim <- aggregate(clim_world_sf, by = admin, FUN = mean)