OHI Science | Citation policy
The Ultraviolet Radiation pressure layer is generated from daily data on Local Noon Erythemal UV Irradiance (mW/m2) derived from satellite observations.
One additional year of data was added to mazu.
https://disc.gsfc.nasa.gov/datasets/OMUVBd_V003/summary https://acdisc.gesdisc.eosdis.nasa.gov/data/Aura_OMI_Level3/OMUVBd.003/ Reference: The Ultraviolet Radiation pressures layer uses the Aura OMI Global Surface UVB Data Product.
Native Data Resolution: 1 degree
Values: Level-3 OMI Surface UV Irradiance and Erythemal Dose- OMUVBd
Time Range: Daily data from 2005 - 2020 (10/1/2004 through 08/10/2021, but only full years of data are used) Format: NetCDF HDF5 (.he5.nc) Downloaded: August 10, 2021
knitr::opts_chunk$set(message = FALSE, warning = FALSE)
## rhdf5 package for working with HDF5 files, from bioconductor: http://bioconductor.org/packages/release/bioc/html/rhdf5.html
#if (!requireNamespace("BiocManager", quietly = TRUE))
# install.packages("BiocManager")
#BiocManager::install()
#BiocManager::install("rhdf5")
library(raster)
source('http://ohi-science.org/ohiprep_v2021/workflow/R/common.R')
library(ncdf4)
library(rgdal)
library(sf) # use simple features rather than sp
library(rhdf5)
library(ggplot2)
library(RColorBrewer)
library(foreach)
library(doParallel)
library(dplyr)
library(readr)
library(stringr)
library(httr)
library(lubridate)
library(googleVis)
library(animation)
library(plotly)
## update these 3 to reflect current assessment year, or whichever year data is being used!!!
data_yr <- "d2021"
raw_data_dir <- file.path(dir_M, "git-annex/globalprep/_raw_data/NASA_OMI_AURA_UV", data_yr)
int_sp_data <- file.path(dir_M, "git-annex/globalprep/prs_uv/v2021/int") # intermediate spatial data location
out_dir <- file.path(dir_M, sprintf("git-annex/globalprep/prs_uv/v2021/output"))
## years of data we are using for this data layer
yrs <- c(2005:2020)
mths <- str_pad(1:12, 2, "left", pad = "0")
days_full <- seq(1, 358, 7)
## global ocean raster at 1km for resampling/projecting purposes
ocean <- raster(file.path(dir_M, "model/GL-NCEAS-Halpern2008/tmp/ocean.tif"))
ocean_shp <- st_read(file.path(dir_M, "git-annex/globalprep/spatial/d2014/data"), layer = "regions_gcs")
land <- ocean_shp %>% filter(rgn_typ %in% c("land", "land-disputed", "land-noeez")) %>% st_geometry()
## define mollweide projection CRS
mollCRS <- crs("+proj=moll +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs")
## define colors in spectral format for plotting rasters -- rainbow color scheme
cols <- rev(colorRampPalette(brewer.pal(9, "Spectral"))(255))Data files can be found in the GES DISC EARTHDATA Data Archive. An EarthData account will be required, and will need to be linked to the GES DISC data archive, see instructions here. Download a links list for the variable “ErythemalDailyDose”, in NetCDF format. To get this list, go here, and click on “Subset/Get Data”. For the download method, select “Get file subsets using OPeNDAP”; you don’t need to change the date range (that will happen in the code below) and you don’t need to change the region. From the variables drop down select just the “Erythemal Daily Dose” option and the file format should be netCDF. Note this list is valid for only 2 days so a new one must be generated if data is to be downloaded after that time frame. Once the links list .txt file has been downloaded, the code below will use the list and earthdata login info (username and password) to download the new data files into the raw data directory. The naming convention of the downloaded files: ‘OMI-Aura’ refers to the instrument, ‘L3’ means it is a level 3 data product, ‘OMUVBd’ is the measurement, the first date is when the data was recorded, the second date and time corresponds to modification/upload of the data.
## need username and password, define in console when prompted (or read from secure file), don't save here!!
usrname <- readline("Type earthdata username:")
pass <- readline("Type earthdata password:")## This took 103.5 minutes in 2019...
## read in file list .text, downloaded from earthdata & saved in destination folder (same as raw_data_dir)
file_list_raw <- read_delim(file.path("/home/shares/ohi/git-annex/globalprep/_raw_data/NASA_OMI_AURA_UV/d2021", "file_list.txt"), delim = "\\", col_names = FALSE)
file_list <- file_list_raw %>%
mutate(url_str = as.character(X1)) %>%
mutate(check_netcdf = str_match(url_str, pattern = "http.*OMI-Aura_L3-OMUVBd.*nc")) %>%
filter(!is.na(check_netcdf)) %>%
select(url_str)
## set up timekeeping for data download
t0 = Sys.time()
n.pym = length(file_list$url_str)
i.pym = 0
## download the data
for(i in 1:nrow(file_list)){
#i = 1
url = as.character(file_list[i,])
name_raw_file = substr(url, 88, 144)
if(file.exists(file.path(raw_data_dir, "data", name_raw_file)) != TRUE){
x = httr::GET(url, authenticate(usrname, pass, type = "basic"), verbose(info = TRUE, ssl = TRUE))
bin = content(x, "raw")
writeBin(bin, file.path(raw_data_dir, "data", name_raw_file)) # gnutls_handshake() failed: Handshake failed
i.pym <- i.pym + 1
min.done <- as.numeric(difftime(Sys.time(), t0, units="mins"))
min.togo <- (n.pym - i.pym) * min.done/i.pym
print(sprintf("Retrieving %s of %s. Minutes done=%0.1f, to go=%0.1f",
i.pym, n.pym, min.done, min.togo)) # approx time remaining for data download
}else{
print(sprintf("Skipping %s of %s. Already done.",
i.pym, n.pym))
}
}
## tip: after downloading, check all .he5.nc files are about the same size i.e. they all downloaded properly/fully
# v2021: all look good ## list and check missing dates in NetCDF data files
## list all files from raw data folder
files <- list.files(file.path(raw_data_dir, "data"), pattern = "OMI-Aura_L3-OMUVBd.*.he5.nc$", full.names = TRUE) # netcdf not hdf
files <- files[substr(files, 96, 99) %in% yrs] # select only files for yrs we want
## check all days are there; should go from Oct 2004 through Dec 31 of last full year of data
files_df <- files %>%
data.frame() %>%
rename(fullname = ".") %>% # View(files_df)
mutate(post_modify_date = substr(fullname, 111, 119),
year = substr(fullname, 96, 99),
mo = substr(fullname, 101, 102),
dy = substr(fullname, 103, 104)) %>%
mutate(date = paste(year, mo, dy, sep = "-"),
wk_of_yr = lubridate::week(as_date(date))) %>%
group_by(year, wk_of_yr) %>%
mutate(nday_in_wk = n()) %>%
ungroup()
check_ndays <- files_df %>%
group_by(year, mo) %>% # group_by(year) %>%
summarize(ndays = n()) #%>% # View(check_ndays)
#filter(ndays < 28)
## For some reason June 2016 is missing 14 days (June 1-14, 2016).. I checked the website and those files don't exist. It is just a gap in the data. The rest seem ok. - v2021Calculate weekly means and standard deviations across all years:
## for every week in each year of the time series, calculate weekly mean and standard deviation
registerDoParallel(3)
## note: 2016 wk 22 has only 4 layers (4th layer all NAs) and length(days)=52; missing all week 23
## note: so far some 2005 files appear to have been downloaded incorrectly... replace them manually or just retry download??
foreach (yr = yrs) %dopar% { # j = 22; yr = 2016
#yr = 2005; #j = 1
l <- files[substr(files, 96, 99) == yr]
days_df <- files_df %>%
filter(year == yr) %>%
select(wk_of_yr, nday_in_wk) %>%
distinct() %>% # select just distinct weeks with number of data days they contain
tidyr::complete(wk_of_yr = seq(1:53)) %>% # possible max 53 weeks
mutate(nday_in_wk = replace(nday_in_wk, is.na(nday_in_wk), 0)) %>% # zeros if no data
mutate(lag_nday = lag(nday_in_wk),
lag_nday = replace(lag_nday, is.na(lag_nday), 1),
doy = cumsum(lag_nday)) # day-of-year for start of each week of data
days <- days_df$doy
weeks <- days_df$wk_of_yr
for (j in weeks[-length(weeks)]) { # print(days[j]:(days[j+1]-1)) # checking without foreach+dopar
## gapfill for weeks with 1 or fewer days using prev + subseq. weeks
if(days_df$nday_in_wk[j] < 2){
wk_files <- l[days[j-1]:(days[j+2]-1)] # gapfilling
} else {
wk_files <- l[days[j]:(days[j+1]-1)]
}
rasters <- raster(wk_files[1], varname = "ErythemalDailyDose")
for(i in wk_files[-1]){
r <- raster(i, varname = "ErythemalDailyDose")
rasters <- stack(rasters, r)
}
uv_week <- rasters
week = str_pad(weeks[j], 2, "left", pad = "0")
week_mean <- calc(uv_week, fun = function(x) {mean(x, na.rm = TRUE)},
filename = sprintf("%s/weekly_means/weekly_means_%s_%s.tif",
int_sp_data, yr, week), overwrite = TRUE)
week_sd <- calc(uv_week, fun = function(x) {sd(x, na.rm = TRUE)},
filename = sprintf("%s/weekly_sd/weekly_sd_%s_%s.tif",
int_sp_data, yr, week), overwrite = TRUE)
week_mean_sd <- overlay(week_mean, week_sd, fun = function(x, y) {x + y},
filename = sprintf("%s/weekly_mean_sd/weekly_mean_sd_%s_%s.tif",
int_sp_data, yr, week), overwrite = TRUE)
}
}
## get weekly climatologies across all years in the time series
names_weekly <- list.files(file.path(int_sp_data, "weekly_means"), full.names = TRUE)
match_weeks <- substr(names_weekly, 87, 92) %>% unique()
## check all weeks expected to be there are there
names_weekly_df <- names_weekly %>%
data.frame() %>%
rename(fullname = ".") %>%
mutate(yr = substr(fullname, 82, 85),
wk = substr(fullname, 87, 88)) # View(names_weekly_df)
tmp <- names_weekly_df %>% # View(tmp)
select(yr, wk) %>%
group_by(yr) %>%
summarize(maxwk = max(wk))
foreach(i = match_weeks) %dopar% {
w = names_weekly[(substr(names_weekly, 87, 92) == i)] %>% stack()
m <- calc(w, fun = function(x){mean(x, na.rm = TRUE)},
filename = sprintf("%s/weekly_climatologies/mean_week_%s", int_sp_data, i),
overwrite = TRUE)
sd <- calc(w, fun = function(x){sd(x, na.rm = TRUE)},
filename = sprintf("%s/weekly_climatologies/sd_week_%s", int_sp_data, i),
overwrite = TRUE)
m_sd <- overlay(m, sd, fun = function(x, y){x + y},
filename = sprintf("%s/weekly_climatologies/mean_sd_week_%s", int_sp_data, i),
overwrite = TRUE) ## climatologies based on mean & sd of all years, additional year each assessment...
}Compare each week in each year to the climatology for that week. The climatology is equal to the mean plus one standard deviation.
## loop to calculate annual positive anomalies
foreach (i = yrs) %dopar% {
match_weeks <- names_weekly_df %>% filter(yr == i)
s = stack()
for(j in match_weeks$wk) {
w_mean = raster(sprintf("%s/weekly_means/weekly_means_%s_%s.tif", int_sp_data, i, j)) # mean UV for week j, year i
w_anom = raster(sprintf("%s/weekly_climatologies/mean_sd_week_%s.tif", int_sp_data, j)) # week j climatology
count = overlay(w_mean, w_anom, fun = function(x, y){ifelse(x > y, 1, 0)}) # compare to average anomaly for that week
s = stack(s, count)
}
year = calc(s, fun = function(x){sum(x, na.rm = TRUE)},
filename = sprintf("%s/annual_anomalies_diff/annual_pos_anomalies_%s.tif", int_sp_data, i),
overwrite = TRUE) ## each assessment get new year of data, another layer in this calculation...
}Calculate the difference in total number of anomalies over a 5 year period compared to the first 5 years (2005-2009)
l <- list.files(file.path(int_sp_data, "annual_anomalies_diff"),
pattern = "anomalies", full.names = TRUE)
## reference period is 2005-2009
ref <- l[1:5] %>% stack() %>%
calc(., fun = function(x){sum(x, na.rm = TRUE)})
plot(ref, col = cols, axes = FALSE, main = "Total Number of Anomalies 2005-2009")
plot(land, add = TRUE)
registerDoParallel(4)
foreach(i = 2005:(max(yrs) - 4)) %dopar% {
## calculate difference between total number of anomalies in recent and historical (2005-2009) time periods
y <- i:(i + 4)
s <- stack(l[str_sub(l, -8, -5) %in% y]) %>% sum(., na.rm = TRUE)
diff <- s - ref
projection(diff) <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
writeRaster(diff,
filename = sprintf("%s/annual_anomalies_diff/annual_diff_%s_%s.tif",
int_sp_data, y[1], y[5]),
overwrite = TRUE)
}ref_files <- list.files(file.path(int_sp_data, "annual_anomalies_diff"),
full.names = TRUE, pattern = "diff")
## this says reference point is 41
read.csv("../../supplementary_information/v2016/reference_points_pressures.csv", stringsAsFactors = FALSE) %>%
filter(pressure == "Ultraviolet Radiation Anomalies") %>%
dplyr::select(ref_point) %>%
as.numeric(.$ref_point)
## get the reference point (resc_num = rescale number, excludes baseline 2005_2009 tif)
vals <- c()
for(i in 2006:(max(yrs) - 4)){
max_yr <- i + 4
m <- raster(ref_files[str_sub(ref_files, -13, -10) == i]) %>%
mask(as(land, "Spatial"), inverse = TRUE) %>% getValues()
vals <- c(vals, m)
}
resc_num <- stats::quantile(vals, prob = 0.9999, na.rm = TRUE) # 42 for v2018; 46 for v2019, 49 for v2020, 47 for v2021
## rescale using reference point
for(file in ref_files){
print(file)
the_name <- gsub(".tif", "", basename(file))
m <- raster(file) %>%
calc(fun = function(x){ifelse(x > 0, ifelse(x > resc_num, 1, x/resc_num), 0)},
filename = file.path(int_sp_data, sprintf("rescaled/%s_rescaled.tif", the_name)),
overwrite = TRUE)
}
resc_files <- list.files(file.path(int_sp_data, "rescaled"),
full.names = TRUE, pattern = "rescaled.tif")
## resample to ocean raster scale (~1km) and then mask
registerDoParallel(4)
foreach(i = 2005:(max(yrs) - 4)) %dopar% {
max_yr <- i + 4
if(file.exists(file.path("/home/shares/ohi/git-annex/globalprep/prs_uv/v2021/output",
sprintf("uv_%s_%s-2005_2009_mol_1km.tif", i, max_yr))) != TRUE){
mol1km_masked <- resc_files[str_sub(resc_files, -22, -19) == i] %>%
raster() %>%
projectRaster(crs = mollCRS, over = TRUE, method = "ngb") %>%
resample(ocean, method = "ngb") %>%
mask(ocean,
filename = file.path("/home/shares/ohi/git-annex/globalprep/prs_uv/v2021/output",
sprintf("uv_%s_%s-2005_2009_mol_1km.tif", i, max_yr)),
overwrite = TRUE)
}else{
print(sprintf("Skipping %s to %s. Already done.",
i, max_yr))
}
}mol1km_masked <- list.files(out_dir, pattern = "uv_.*_mol_1km.tif", full.names = TRUE)
out <- raster(mol1km_masked[length(mol1km_masked)])
plot(out, box = FALSE, col = cols, axes = FALSE, main = "UV Pressure Layer \n OHI 2021")## load raster/zonal data
zones <- raster(file.path(dir_M, "git-annex/globalprep/spatial/v2017/regions_eez_with_fao_ant.tif"))
rgn_data <- read.csv("~/github/ohi-global/eez/spatial/regions_list.csv") %>%
filter(rgn_id <= 250)
## get raster data
rasts <- list.files(out_dir, full = TRUE, pattern = "uv_.*_mol_1km.tif")
pressure_stack <- stack(rasts)saveGIF({
for(i in 1:nlayers(pressure_stack)){
n = sprintf("UV Pressure %s",
substr(names(pressure_stack[[i]]), 4, 12))
plot(pressure_stack[[i]],
zlim = c(0, 1), # don't forget to fix the zlimits
axes = FALSE, box = FALSE, col = cols,
main = n)}},
ani.width = 750,
ani.height = 400,
movie.name = "uv.gif")## extract data for each region
regions_stats <- zonal(pressure_stack, zones, fun = "mean", na.rm = TRUE,
progress = "text") %>% data.frame()
write.csv(regions_stats, "int/uv_mean_rgn.csv", row.names = FALSE)
## check regions are all present or missing as expected
setdiff(regions_stats$zone, rgn_data$rgn_id) # high seas and Antarctica
# 260 261 262 263 264 266 267 269 270 272 273 274 275 276 277 248100 248200 248300 248400
# 248500 248600 258410 258420 258431 258432 258441 258442 258510 258520 258600 258700 288100 288200 288300
setdiff(rgn_data$rgn_id, regions_stats$zone) # Antarctica is 213
data <- merge(rgn_data, regions_stats, all.y = TRUE, by.x = "rgn_id", by.y = "zone") %>%
tidyr::gather("year", "pressure_score", starts_with("uv")) %>%
filter(rgn_id <= 250) # filter out non OHI global regions
uv_data <- data %>%
mutate(year = substring(year, 9, 12)) %>%
mutate(year = as.numeric(year)) %>%
dplyr::select(rgn_id, rgn_name, year, pressure_score)## visualize data using googleVis plot
plotData <- uv_data %>%
dplyr::select(rgn_name, year, pressure_score)
motion = gvisMotionChart(plotData,
idvar = "rgn_name",
timevar = "year")
plot(motion)
print(motion, file = "uv.html")
## save data layer
uv_data <- uv_data %>%
dplyr::select(rgn_id, year, pressure_score)
write.csv(uv_data, "output/uv.csv", row.names = FALSE)uv_data <- read.csv("output/uv.csv")
old <- uv_data %>%
filter(year == 2020) %>% # data years lag assessment yrs by 1
select(-year, new_pressure = pressure_score) %>%
left_join(read.csv("../v2020/output/uv.csv") %>% filter(year == 2019),
by = c("rgn_id")) %>%
rename(old_pressure = pressure_score)
plot_diff <- ggplot(old, aes(new_pressure, old_pressure, label = rgn_id)) +
geom_point() +
geom_abline()
ggplotly(plot_diff)
# ggplot(old %>% mutate(difference = new_pressure - old_pressure), aes(difference)) + geom_histogram(binwidth = 0.002)Niilo Kalakoski, Panu Lahtinen, Jari Hovila (May 2016) OMI/Aura Surface UVB Irradiance and Erythemal Dose Daily L3 Global Gridded 1.0 degree x 1.0 degree V3, NASA Goddard Space Flight Center.