OHI Science | Citation policy
This data prep markdown calls a series of scripts to download and analyze sea ice data for the habitat subgoal and coastal protection goal. See data layers documentation for more information.
No changes to methods from previous assessment. Transitioned almost entirely to using the sf package for the spatial data; the one dependence on sp left in ObtainingData.R is OHIregion_points calculation since rasterToPoints() turns raster first into a spatialPointsDataFrame and then sf.
Reference:
Cavalieri, D.J., Parkinson, C.L., Gloersen, P. and Zwally, H. (1996). Sea ice concentrations from Nimbus-7 SMMR and DMSP SMM/I-SSMIS passive microwave data. 1979-2014. Boulder, Colorado, USA: NASA National Snow and Ice Data Center Distributed Active Archive Center (http://dx.doi.org/10.5067/8GQ8LZQVL0VL).
Downloaded: May 18, 2018
Description: Monthly sea ice extent data.
Data are in the following FTP directory: ftp://sidads.colorado.edu/pub/DATASETS/nsidc0051_gsfc_nasateam_seaice/final-gsfc/
Within the final-gsfc directory are north and south directories that contain data files, and a browse directory that contains browse image files. Daily and monthly data are further separated into directories named daily and monthly. For final daily data, there is also one directory for each year of available data. For example, all of the north daily data for 1990 are in a directory named /nsidc0051_gsfc_nasateam_seaice/final-gsfc/north/daily/1990/.
The directory structure is illustrated below; not all directory levels are shown. The structure for each south directory matches that of the corresponding north directory. Each browse directory is divided into a structure that reflects that of the data. In this illustration, the year directories underneath final-gsfc are representative placeholders; there are actually many such directories, each named for the year of data it contains, such as 1987, 2000, etc.
/nsidc0051_gsfc_nasateam_seaice
. . /final-gsfc
. . . . /browse
. . . . . . /north
. . . . . . . . /daily
. . . . . . . . . . /year
. . . . . . . . /monthly
. . . . . . /south
. . . . /north
. . . . . . /daily
. . . . . . . . /year
. . . . . . /monthly
. . . . /south
For complete documentation and more information about data access, please see:
http://nsidc.org/data/nsidc-0051.html
If you wish to be notified of updates or corrections to these data, please register with NSIDC User Services by sending e-mail to: nsidc@nsidc.org
Identify yourself as a user of “Sea Ice Concentrations from Nimbus-7 SMMR and DMSP SSM/I-SSMIS Passive Microwave Data (NSIDC-0051).” Include your name, e-mail address, postal address, and telephone number.
If you have questions, please contact NSIDC User Services.
CONTACT INFORMATION: User Services National Snow and Ice Data Center CIRES, 449 UCB University of Colorado Boulder, CO USA 80309-0449 Phone: +1 303-492-6199 Fax: +1 303-492-2468 E-mail: nsidc@nsidc.org
Time range: 1970-2017
Load all relevant libraries, establish/define parameters and commonly used pathnames. Manually change scenario and data years in file pathnames code chunk to reflect the most recent data (d) and current assessment year (v) in setup code chunk.
## install R packages where necessary
pkg <- c("raster", "fasterize", "sf", "sp", "rgdal", "fields") # "fields" for colors in Status_Trend.R
new.pkg <- pkg[!(pkg %in% installed.packages())]
if (length(new.pkg)){install.packages(new.pkg)}
if (!("ohicore" %in% installed.packages())){devtools::install_github("ohi-science/ohicore")}
## load libraries, set directories
lapply(c(pkg, "ohicore"), require, character.only = TRUE)
## UPDATE THESE!
assessYear <- "v2018" # change to reflect assessment year
previous_yr <- "v2017" # previous assessment year
data_yr_gci <- "d2018" # change to reflect year of most recently downloaded data
final.year <- 2017 # final year of data (all months); double check nsidc0051 monthly dataset was updated at ftp://sidads.colorado.edu/pub/DATASETS/
## comment out when knitting
#setwd(paste0("globalprep/hab_seaice/", assessYear))source("../../../src/R/common.R") # directory locations
source("../../../src/R/spatial_common.R")## Reading layer `regions_2017_update' from data source `/home/shares/ohi/git-annex/globalprep/spatial/v2017' using driver `ESRI Shapefile'
## Simple feature collection with 526 features and 7 fields
## geometry type: MULTIPOLYGON
## dimension: XY
## bbox: xmin: -18040060 ymin: -9020048 xmax: 18040080 ymax: 9020048
## epsg (SRID): NA
## proj4string: +proj=moll +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defsThese maps of the OHI regions were made by the PreparingSpatialFiles.R script. If there are changes to the OHI regions, the PreparingSpatialFiles.R script will need to be run. Additionally, it is critical to walk through the ObtainingData.R script if any of the spatial files have been modified (this saves the files as spatial points).
The original polygon files are used from: git-annex/globalprep/_raw_data/NSIDC_SeaIce/v2015. However, we use the fasterize package to rasterize the polygons and save them to: git-annex/globalprep/_raw_data/NSIDC_SeaIce/v2017.
maps <- file.path(dir_M, "git-annex/globalprep/_raw_data/NSIDC_SeaIce/v2018") # change year to current assessmentEstablish: CRS, website to collect data, data selection parameters. Filename format for final monthly data is nt_YYYYMM_SSS_vVV_R.bin. Parameters will be used to scrape the data from the web in ObtainingData script.
pixel = 25000 # pixel dimension in meters for both x and y
prj.n = "+proj=stere +lat_0=90 +lat_ts=70 +lon_0=-45 +k=1 +x_0=0 +y_0=0 +a=6378273 +b=6356889.449 +units=m +no_defs"
prj.s = "+proj=stere +lat_0=-90 +lat_ts=-70 +lon_0=0 +k=1 +x_0=0 +y_0=0 +a=6378273 +b=6356889.449 +units=m +no_defs"
prj.mol = "+proj=moll +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"
# URL base (ub), filename format for final monthly data is nt_YYYYMM_SSS_vVV_R.bin
ub.n = "ftp://sidads.colorado.edu/pub/DATASETS/nsidc0051_gsfc_nasateam_seaice/final-gsfc/north/monthly"
ub.s = "ftp://sidads.colorado.edu/pub/DATASETS/nsidc0051_gsfc_nasateam_seaice/final-gsfc/south/monthly"
poles = c("n","s")
years = c(1979:final.year) # full range of data
months = 1:12
n.pym = length(poles)*length(years)*length(months)
i.pym = 0Collects the data for each month/year from the website and add to raster stack that is saved in tmp folder as: n_rasters_points.rdata or s_rasters_points.rdata. And, if it doesn’t already exist, it converts the region shapefile into a raster points file. See ObtainingData.R script for more details.
source("ObtainingData.R")Using the data from the .rdata files created from Function 1 with the ObtainingData.R script, this function calculates status and trend for shoreline ice and ice edge habitat. Data is saved in intermediate (int) folder:
n_IceEdgeHabitat.csv, s_IceEdgeHabitat.csvn_IceShoreProtection.csv, s_IceShoreProtection.csvref.years <- 1979:2000
source("Status_Trend.R") # calculates status and trend for shoreline ice and ice edge habitatRead in ice edge habitat and ice shore protection csv-format datasets, remove anamolous eez regions with minimal ice cover, remove disputed regions. Bind these datasets and convert to units of km^2. Save seaice health, extent, trend, and extent data.
n_edge <- read.csv("int/n_IceEdgeHabitat_ref1979to2000.csv")
s_edge <- read.csv("int/s_IceEdgeHabitat_ref1979to2000.csv")
edge <- rbind(n_edge, s_edge)
edge <- edge %>%
dplyr::filter(Reference_avg1979to2000monthlypixels != 0) %>%
dplyr::filter(!(rgn_id %in% c(59, 141, 219, 4, 172, 94))) %>% # anomolous eez regions with very little ice cover
dplyr::filter(!(rgn_id %in% c("248300","258510","258520","258600","258700"))) %>% # cut due to minimal ice (<200km2/yr - avg of months)
dplyr::filter(rgn_nam != "DISPUTED") %>%
dplyr::mutate(habitat="seaice_edge")
n_shore <- read.csv("int/n_IceShoreProtection_ref1979to2000.csv")
s_shore <- read.csv("int/s_IceShoreProtection_ref1979to2000.csv")
shore <- rbind(n_shore, s_shore)
shore <- shore %>%
dplyr::filter(Reference_avg1979to2000monthlypixels != 0) %>%
dplyr::filter(!(rgn_id %in% c(59, 89, 177, 178))) %>% # anomolous eez regions with very little ice cover
dplyr::filter(rgn_nam != "DISPUTED") %>%
dplyr::mutate(habitat = "seaice_shoreline")
data <- rbind(edge, shore)
data <- data %>%
dplyr::mutate(km2 = Reference_avg1979to2000monthlypixels/12 * (pixel/1000)^2)
write.csv(data, "int/sea_ice.csv", row.names = FALSE)
## health data
health <- data %>%
dplyr::filter(rgn_typ == "eez") %>%
dplyr::select(rgn_id, habitat, dplyr::starts_with("pctdevR")) %>%
tidyr::gather("year", "health", -(1:2)) %>%
dplyr::mutate(year = substring(year, 9, 12)) %>%
dplyr::mutate(year = as.numeric(year)) %>%
dplyr::mutate(health = ifelse(health > 1, 1, health))
write.csv(health, "output/hab_ice_health_eez.csv", row.names = FALSE)
## trend data
trend <- data %>%
dplyr::filter(rgn_typ == "eez") %>%
dplyr::select(rgn_id, habitat, dplyr::starts_with("Trend")) %>%
tidyr::gather("year", "trend", -(1:2)) %>%
dplyr::mutate(year = substring(year, 13, 16)) %>%
dplyr::mutate(year = as.numeric(year)) %>%
dplyr::mutate(trend = trend * 5) %>%
dplyr::mutate(trend = ifelse(trend > 1, 1, trend)) %>%
dplyr::mutate(trend = ifelse(trend < (-1), -1, trend))
write.csv(trend, "output/hab_ice_trend_eez.csv", row.names = FALSE)
## seaice extent data
extent <- data %>%
dplyr::filter(rgn_typ == "eez") %>%
dplyr::mutate(year = 2016) %>% # extent not updated each year (historic extent of the sea ice habitat); updated last 2016 bc of source methods
dplyr::select(rgn_id, habitat, year, km2)
write.csv(extent, "output/hab_ice_extent_eez.csv", row.names = FALSE) # save seaice extent data## health
health <- read.csv("output/hab_ice_health_eez.csv")
ice_health_eez <- read.csv(sprintf("../%s/output/hab_ice_health_eez.csv", previous_yr)) %>% # compare to last year's data
dplyr::rename(health_prev_assess = health) %>%
dplyr::left_join(health, by = c("rgn_id", "habitat", "year")) %>%
na.omit("health")
plot(ice_health_eez$health_prev_assess, ice_health_eez$health)
abline(0, 1, col="red")
## trend
trend <- read.csv("output/hab_ice_trend_eez.csv")
ice_trend_eez <- read.csv(sprintf("../%s/output/hab_ice_trend_eez.csv", previous_yr)) %>%
dplyr::rename(trend_prev_assess = trend) %>%
dplyr::left_join(trend, by = c("rgn_id", "habitat", "year")) %>%
na.omit("trend")
plot(ice_trend_eez$trend_prev_assess, ice_trend_eez$trend)
abline(0, 1, col="red")
## extent
extent <- read.csv("output/hab_ice_extent_eez.csv", stringsAsFactors = FALSE) # compare to last year's data
ice_extent_eez <- read.csv(sprintf("../%s/output/hab_ice_extent_eez.csv", previous_yr), stringsAsFactors = FALSE) %>%
dplyr::rename(km2_prev_assess = km2) # not sure why left_join won't work- results in NAs in km2 col
plot(ice_extent_eez$km2_prev_assess, extent$km2)
abline(0, 1, col="red")
## pair with region names
regions <- rgn_master %>%
dplyr::select(rgn_id = rgn_id_2013, rgn_name = rgn_nam_2013) %>%
unique()
## make sure sea ice health across regions makes intuitive sense...
data <- read.csv("output/hab_ice_health_eez.csv") %>%
dplyr::left_join(regions, by = "rgn_id") %>%
dplyr::arrange(habitat, health)There was no gapfilling for these data. Created gapfill files with values of 0. Note: all layers need a gf file, eventhough if there was no gapfilling. In this case the gapfill value is 0 for every region.
hab_ice_health_gf <- read.csv("output/hab_ice_health_eez.csv")%>%
mutate(gapfilled = 0) %>%
select(rgn_id, year, gapfilled)
write.csv(hab_ice_health_gf, "output/hab_ice_health_eez_gf.csv", row.names=FALSE)
hab_ice_extent_gf <- read.csv("output/hab_ice_extent_eez.csv")%>%
mutate(gapfilled = 0) %>%
select(rgn_id, year, gapfilled)
write.csv(hab_ice_health_gf, "output/hab_ice_extent_eez_gf.csv", row.names=FALSE)
hab_ice_trend_gf <- read.csv("output/hab_ice_trend_eez.csv")%>%
mutate(gapfilled = 0) %>%
select(rgn_id, year, gapfilled)
write.csv(hab_ice_health_gf, "output/hab_ice_trend_eez_gf.csv", row.names=FALSE)