ohi logo
OHI Science | Citation policy

1 Summary

This layer preparation script does the following for newly available SLR data.

  • Clips all monthly rasters to the coast using a 3 nautical mile offshore buffer
  • Calculates annual mean sea level anomaly rasters from monthly data
  • Rescales values from 0 to 1 using the reference point
  • Sets to zero all negative values, indicating decreases in mean sea level
  • Resamples raster to ~ 1km2 and reproject to Molleweide

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.

2 Updates from previous assessment

One additional year of data, 2018, was added.


3 Data

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


4 Methods

4.3 Data Prep

4.3.1 Clip data to coastal cells

All NetCDF files for each month are rasterized.

The raw monthly data looks like this:

The following code is used to:

  1. Rasterize each monthly NetCDF file
  2. Rotate each raster so that the Atlantic Ocean is centered in the raster

The output is saved in the folder int/msla_monthly

4.5 Changing the projection and masking

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")))

4.6 Reference Point

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 ####

5 Results

## 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")