OHI Science | Citation policy
This layer preparation script does the following for newly available SLR data.
This process is completed entirely within this script. The raw data is downloaded externally and held on a server at NCEAS. Although the raw data is not provided, this script can be used on the data downloaded from Aviso here. You will need to register with Aviso in order to get a username and password for data access.
One additional year of data, 2018, was added.
The source data are monthly mean sea level anomalies, in meters. These anomalies are calculated by subtracting the current absolute sea level for each month from the average sea level for that month calculated from 1993 - 2012.
Reference: The altimeter products were produced and distributed by Aviso (http://www.aviso.altimetry.fr/), as part of the Ssalto ground processing segment. AVISO MSLA heights, monthly means
Downloaded: July 30, 2019 (for 2018 data)
Description: Monthly mean sea level anomaly (meters above mean sea level)
Native data resolution: 0.25 degree grid cells
Time range: January 1993 - December 2018
Format: NetCDF
Citation information The altimeter products were produced and distributed by Aviso (http://www.aviso.altimetry.fr/), as part of the Ssalto ground processing segment. AVISO MSLA heights, monthly means
knitr::opts_chunk$set(fig.width = 6, fig.height = 4, fig.path = 'figs/',message = FALSE, warning = FALSE)
#setwd("~/github/ohiprep_v2019/globalprep/prs_slr/v2019") # change year to reflect current assessment year
## setting up provenance
# devtools::install_github('oharac/provRmd')
# library(provRmd)
# prov_setup()
## load and install packages
pkg <- c( "httr", "R.utils", "raster", "tidyverse", "sf", "RColorBrewer", "ggplot2", "googleVis",
"maps", "parallel", "foreach", "doParallel", "fasterize", "rasterVis")
new.pkg <- pkg[!(pkg %in% installed.packages())]
if (length(new.pkg)){install.packages(new.pkg)}
lapply(pkg, require, character.only = TRUE)
source('../../../workflow/R/common.R')
library(here)
## define paths and variables to use throughout data prep
scen_year <- 2019 # change to reflect assessment year!
dir_anx_aviso <- file.path(dir_M, "git-annex/globalprep/_raw_data/AVISO_slr") # raw data file path
dir_prs_slr <- sprintf("%s/git-annex/globalprep/prs_slr/v%s", dir_M, scen_year)
p4s_wgs84 <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0" # proj4string
cols <- rev(colorRampPalette(brewer.pal(9, "Spectral"))(255)) # rainbow color scheme for maps
mollCRS <- CRS("+proj=moll") # mollweide projection
## read in ocean raster with cells at 1km -- template for resampling (held on an NCEAS server)
ocean <- raster(file.path(dir_M, "git-annex/globalprep/spatial/v2017/ocean.tif"))## need AVISO username and password, define in console when prompted, don't save here!!
## if either password or username contain the @ symbol, replace with '%40' (percent-encoding for reserved characters)
userpasswd <- readline("Type AVISO username and password, separated by colon no space:")## the is a 'download_data.R' script in the Mazu raw data folder but that routine was not working for me...
year <- 2018 # data year to download
months <- str_pad(1:12, 2, pad = "0")
## download data from FTP and collect each month/year in raster stack
url <- "ftp://ftp-access.aviso.altimetry.fr/climatology/global/delayed-time/monthly_mean"
ub <- sprintf("ftp://%s@%s", userpasswd, substr(url, 7, 80))
## download is quite fast, takes less that 3 minutes for 12 files
for (mo in months){
#mo="11" # for testing
## retrieve the data (compressed, gzipped files)
u <- sprintf("%s/dt_global_allsat_msla_h_y2018_m%s.nc.gz", ub, mo)
u_filename <- file.path(dir_anx_aviso, "d2019", substr(u, 117, 155))
res <- httr::GET(u, write_disk(u_filename))
}
closeAllConnections()
zipfiles <- list.files(file.path(dir_anx_aviso, paste0("d", scen_year)),
full.names = TRUE, pattern = "*nc.gz")
for(zipfile in zipfiles){
message("Unzipping file: ", zipfile)
R.utils::gunzip(zipfile, remove = TRUE, skip = TRUE)
}All NetCDF files for each month are rasterized.
## d2016/msla_monthly_mean has data for 1993-2015
## then include list.files for d2017 through the data folder for current scenario year
nc_files <- c(list.files(file.path(dir_anx_aviso, "d2019"),
full.names = TRUE, pattern = ".nc"),
list.files(file.path(dir_anx_aviso, "d2018"),
full.names = TRUE, pattern = ".nc"),
list.files(file.path(dir_anx_aviso, "d2017"),
full.names = TRUE, pattern = ".nc"),
list.files(file.path(dir_anx_aviso, "d2016/msla_monthly_mean"),
full.names = TRUE, pattern = ".nc"))The raw monthly data looks like this:
plot(raster(nc_files[1]), col = cols, axes = F,
main = paste("Year", substr(nc_files[1], 90, 93), "Month", substr(nc_files[1], 96, 97)))The following code is used to:
The output is saved in the folder int/msla_monthly
registerDoParallel(10)
## parallel forloop function that rotates each monthly file, sets the long/lat projection, and keeps only coastal cells - saved to GitHub
foreach(file = nc_files) %dopar% {
m_yr <- substr(file, nchar(file)-10, nchar(file)-3)
## read in month raster
r <- raster(file) %>%
rotate()
## define projection of the raster before reprojecting
projection(r) <- "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"
## write raster to int folder in prs_slr
fp <- sprintf("%s/int/msla_monthly/msla_monthly_%s.tif", dir_prs_slr, m_yr)
writeRaster(r, filename = fp, overwrite = TRUE)
}Annual mean sea level anomaly rasters are calculated from the monthly data.
## will need to create 'msla_annual_mean' folder in 'int' (also 'msla_annual_mol' and 'msla_monthly' folders)
msla_files <- list.files(sprintf("%s/int/msla_monthly", dir_prs_slr),
full.names = TRUE)
maxyr <- substr(msla_files, 83, 86) %>% as.numeric() %>% max()
## stack all rasters for this year, and calc annual mean, then write as raster
registerDoParallel(6)
foreach(yr = c(1993:maxyr)) %dopar% {
files <- msla_files[str_detect(msla_files, as.character(yr))]
rast_annual_mean <- stack(files) %>%
calc(mean, na.rm = TRUE) %>%
writeRaster(filename = sprintf("%s/int/msla_annual_mean/msla_annual_%s.tif", dir_prs_slr, yr),
overwrite = TRUE)
}Since we are only interested in the increase in sea level near the coasts, we apply a mask to the raster layers that removes all cells farther than 3nm offshore. This mask was created previously for the OHI global 2016 assessment.
## 3nm offshore raster to select only nearshore cells
#ocean_mask_prev <- sprintf("%s/int/ocean_mask.tif", dir_prs_slr, scen_year-1)
ocean_mask_prev <- file.path(dir_M, "git-annex/globalprep/prs_slr/v2018/int/ocean_mask.tif")
if(file.exists(ocean_mask_prev)){
file.copy(ocean_mask_prev, file.path(dir_prs_slr, "int"))
} else {
poly_3nm <- read_sf(file.path(dir_M, "git-annex/Global/NCEAS-Regions_v2014/data"), "rgn_offshore3nm_mol")
poly_3nm[duplicated(poly_3nm$rgn_id), ] # check to make sure there are no duplicated regions
## create rasterize 3 nautical miles offshore rasters if cannot copy from previous assessment folder
s <- fasterize(poly_3nm, ocean, field = "rgn_id")
s <- calc(s, fun = function(x) {ifelse(x > 0, 1, NA)})
writeRaster(s, sprintf("%s/int/ocean_mask.tif", dir_prs_slr, scen_year))
}
s <- raster(sprintf("%s/int/ocean_mask.tif", dir_prs_slr, scen_year))
plot(s, col = "red")
## reproject to mollweide
annual_means <- list.files(file.path(dir_prs_slr, "int/msla_annual_mean"), full = TRUE)
foreach(file = annual_means) %dopar% { #file = annual_means[1]
yr <- str_sub(file, -8, -5)
msla_int <- file.path(dir_prs_slr, "int")
rast_data <- raster(file) %>%
projectRaster(crs = mollCRS, over = TRUE) %>%
raster::resample(ocean, method = "ngb",
filename = sprintf("%s/msla_annual_mol/mlsa_annual_mol_%s.tif",
msla_int, yr), overwrite = TRUE)
}
annual_mol <- list.files(file.path(dir_prs_slr, "int/msla_annual_mol"), full = TRUE)
foreach(file = annual_mol) %dopar% { # file = annual_mol[2]
yr <- str_sub(file,-8,-5)
rast <- raster(file)
mask(raster(file), s, filename = sprintf("%s/int/msla_annual_mol_coastal/msla_annual_mol_coastal_%s.tif",
dir_prs_slr, yr), overwrite = TRUE)
}
plot(raster(file.path(dir_prs_slr, "int/msla_annual_mol_coastal/msla_annual_mol_coastal_2010.tif")))The reference point is the 99.99th quantile of the entire data distribution from 1993 - 2015. (This value has been updated due to changes in the source data, previously was 0.246225 m, currently is 0.3359385 m).
coastal_rasts <- list.files(file.path(dir_prs_slr, "int/msla_annual_mol_coastal"), pattern = "tif", full.names = TRUE)
## get data across all years to 2015
## takes a really long times; added foreach dopar to speed...
## doesn't really need to be recalcuated each year unless theres reason to believe source updated past years data
registerDoParallel(8)
vals <- foreach(i = 1993:2015, .combine = c) %dopar% { # i = 1993
coastal_rasts[which(str_sub(coastal_rasts, -8, -5) == i)] %>%
raster() %>%
getValues() %>%
na.omit()
}
ref_point_slr <- quantile(vals, 0.9999)
dir_refpt <- "../../supplementary_information"
if(file.exists(sprintf("%s/v%s/reference_points_pressures.csv", dir_refpt, scen_year))){
## if already created and been partially updated for this assessment, don't want to overwrite with v2016 csv...
ref_tab <- read_csv(sprintf("%s/v%s/reference_points_pressures.csv", dir_refpt, scen_year))
} else {
## grab ref file from v2016 if doesn't exist yet in current assessment 'supplementary information' folder
ref_tab <- read_csv(file.path(dir_refpt, "v2016/reference_points_pressures.csv"))
}
ref_tab$ref_point[ref_tab$pressure == "Sea Level Rise"] <- ref_point_slr # set sea level rise reference to the 99.99 percentile
write.csv(ref_tab, sprintf("%s/v%s/reference_points_pressures.csv", dir_refpt, scen_year), row.names = FALSE)
## grab reference value from the supp_info csv
ref <- read_csv(sprintf("%s/v%s/reference_points_pressures.csv", dir_refpt, "2018")) %>%
filter(pressure == "Sea Level Rise") %>%
.$ref_point %>%
as.numeric()
### Ask about this!!!!
#### Stopped here on July 31 ####Each annual raster is recaled from 0 to 1 using the reference point. If a value is greater than the reference point, it is automatically given a value of 1.
# registerDoParallel(10) # reregister if restarted session between parallel steps above and here
foreach(file = coastal_rasts) %dopar% { # file = coastal_rasts[10]
yr <- str_sub(file, -8,-5)
raster::raster(file) %>%
calc(fun = function(x){ifelse(x < 0, 0, x)}) %>% # set all negative values to 0
calc(fun = function(x){ifelse(x > ref, 1, x/ref)}, # set equal to one if greater than ref, otherwise scale
filename = sprintf("%s/output/slr_%s.tif", dir_prs_slr, yr), overwrite = TRUE) }r <- raster(sprintf("%s/output/slr_%s.tif", dir_prs_slr, scen_year - 1))
plot(ocean, col = "cornsilk2", axes = FALSE, box = FALSE, main = "Sea Level Rise Pressure 2018", legend = FALSE)
plot(r, col = cols, axes = FALSE, box = FALSE, add = TRUE)
r_new <- raster(sprintf("%s/output/slr_%s.tif", dir_prs_slr, scen_year - 1))
r_old <- raster(sprintf("%s/git-annex/globalprep/prs_slr/v%s/output/slr_%s.tif", dir_M, scen_year - 1, scen_year - 2))
rasterVis::histogram(r_old, main = sprintf("Sea Level Pressure %s old data", scen_year - 2))
rasterVis::histogram(r_new, main = sprintf("Sea Level Pressure %s new data", scen_year - 2))## raster/zonal data, zones tifs created & spatial rgn updated in 2017
slr_loc <- file.path(dir_prs_slr, "output")
rasts <- list.files(slr_loc, full.names = TRUE) %>% str_subset(pattern = ".tif$")
stack_slr <- stack(rasts) # read in raster files
zones <- raster(file.path(dir_M, "git-annex/globalprep/spatial/v2017/regions_eez_with_fao_ant.tif"))
rgn_data <- read_sf(file.path(dir_M, "git-annex/globalprep/spatial/v2017"), "regions_2017_update") %>%
st_set_geometry(NULL) %>%
dplyr::filter(rgn_type == "eez") %>%
dplyr::select(rgn_id = rgn_ant_id, rgn_name)
## extract data for each region
## fyi takes awhile... about 2 hours for v2019....
regions_stats <- zonal(stack_slr, zones, fun = "mean", na.rm = TRUE, progress = "text") %>% data.frame()
rgn_stats_long <- foreach(i = zones, .packages = c("raster", "dplyr"), .combine = "rbind") %dopar% {
zonal(stack_slr, zones, fun = "mean", na.rm = TRUE, progress = "text") %>% data.frame()
} ## I don't think this line of code or any of the commented out are necessary... waiting for response from Mel on this.
# t0 <- Sys.time()
# zones_shp <- st_read(file.path(dir_M, "git-annex/globalprep/spatial/v2017"), "regions_2017_update") %>%
# filter(rgn_type == "eez" & rgn_id <= 250 & rgn_id != 213)
# rgn_data <- rgn_data %>% filter(rgn_id <= 250 & rgn_id != 213)
#
# rgn_stats_long <- foreach(i = 1:nlayers(stack_slr), .packages = c("raster", "dplyr"), .combine = "rbind") %dopar% {
# zonal_cmp_fun <- compiler::cmpfun(function(x){
# rgn_mean <- cellStats(rasterize(x, crop(stack_slr[[i]], as.matrix(extent(x))), mask = TRUE),
# stat = "mean", na.rm = TRUE)
# })
# foreach(j = rgn_data$rgn_id[1:20], .packages = c("dplyr"), .combine = "rbind") %dopar% {
# # rgn <- zones_shp %>% filter(rgn_id == j)
# rgn_stats_long <- c(as.integer(str_sub(rasts[i], -8, -5)), j,
# zonal_cmp_fun(zones_shp %>% filter(rgn_id == j)))
# }
# }
# regions_stats_test <- rgn_stats_long %>% data.frame() %>%
# rename(year = X1, rgn_id = X2, pressure_score = X3) # rownames not used ok
# Sys.time() - t0
setdiff(regions_stats$zone, rgn_data$rgn_id) # High Seas regions are in there, makes sense....no land
#[1] 260 261 262 263 264 266 267 269 270 272 273 274 275 276 277
setdiff(rgn_data$rgn_id, regions_stats$zone) #integer(0)
regions_stats <- regions_stats %>%
rename(rgn_id = zone) %>%
filter(rgn_id <= 250) %>%
gather("year", "pressure_score", -1) %>%
mutate(year = as.numeric(as.character(substring(year, 5, 8))))
write.csv(regions_stats, "output/slr.csv", row.names = FALSE)
regions_stats <- read_csv("output/slr.csv")
## visualize data
plotData <- regions_stats %>%
left_join(rgn_data, by = "rgn_id") %>%
dplyr::select(rgn_name, year, pressure_score) %>%
dplyr::arrange(rgn_name, year) %>%
data.frame()
Motion <- gvisMotionChart(plotData, idvar = "rgn_name", timevar = "year")
plot(Motion)
print(Motion, file = "slr.html")new_data <- read.csv("output/slr.csv") %>%
dplyr::select(rgn_id, year, new_pressure_score = pressure_score)
old <- read.csv(sprintf("../v%s/output/slr.csv", scen_year - 1)) %>%
left_join(new_data, by = c("year", "rgn_id"))
plot(old$pressure_score, old$new_pressure_score, ylab = "new score", xlab = "old score")
abline(0, 1, col = "red")There was no gapfilling for these data. Created gapfill files with values of 0.