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: August 10, 2021
Description: Monthly mean sea level anomaly (meters above mean sea level)
Native data resolution: 0.25 degree grid cells
Time range: January 1993 - December 2019 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)
## load packages
library(httr)
library(R.utils)
library(raster)
library(tidyverse)
library(sf)
library(RColorBrewer)
library(ggplot2)
library(googleVis)
library(maps)
library(parallel)
library(foreach)
library(doParallel)
library(fasterize)
library(rasterVis)
# 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 <- 2021 # 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"))These chunks can be used to download the data from AVISO. You will need an AVISO account to do this. If you want to view the files to see if a full previous year of data is available to download, go here: ftp://ftp-access.aviso.altimetry.fr/climatology/global/delayed-time/monthly_mean
or here: ftp://yourlogin@ftp-access.aviso.altimetry.fr/climatology
And you will be prompted to either open up a FTP application, like Cyberduck, or to just plug in your username and password.
If you don’t want to go through the trouble of running the download process in R, you can do it manually by just click-and-drag between two cyberduck-2 browsers (I found this method to be much quicker and easier).
# This layer lags by 2 years unfortunately, since we need a full year of data.
## 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:") ## there is a 'download_data.R' script in the Mazu raw data folder but that routine was not working for me...
year <- 2019 # data year to download; lags by 2 years, since 2020 only goes to may as of now... v2021
months <- str_pad(1:12, 2, pad = "0") #if they upload the rest just change the numbers here to reflect rest of months
## 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="01" # for testing
## retrieve the data (compressed, gzipped files)
u <- sprintf("%s/dt_global_allsat_msla_h_y2019_m%s.nc.gz", ub, mo)
u_filename <- file.path(dir_anx_aviso, "d2021", substr(u, 115, 155)) #updated 2020 to include full filename
res <- httr::GET(u, write_disk(u_filename))
}
closeAllConnections()Unzip the files you’ve just downloaded:
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, "d2021"),
full.names = TRUE, pattern = ".nc"),
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[3]), col = cols, axes = F,
main = paste("Year", substr(nc_files[3], 90, 93), "Month", substr(nc_files[3], 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/v2019/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)
## If not rerunning the above, use this (taken from v2019 reference point csv)
ref_point_slr <- 0.335938483476639
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, "2021")) %>%
filter(pressure == "Sea Level Rise") %>%
.$ref_point %>%
as.numeric()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[26]
yr <- str_sub(file, -8,-5)
if(file.exists(sprintf("%s/output/slr_%s.tif", dir_prs_slr, yr))){
message("skipping")
}else{
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 - 2))
plot(ocean, col = "cornsilk2", axes = FALSE, box = FALSE, main = "Sea Level Rise Pressure 2019", 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()
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 - 2)) %>%
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.
slr_gf <- read.csv("output/slr.csv")%>%
mutate(pressure_score = 0) %>%
rename(gapfilled = pressure_score)
write.csv(slr_gf, "output/slr_gf.csv", row.names=FALSE)