NSIDC_Data_R.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='top'></a>\n",
    "# Working with NSIDC DAAC/NOAA data in R\n",
    "---\n",
    "This notebook was built to provide basic use cases for accessing and working with NSIDC/NOAA data in R. The Jupyter Notebook uses an R kernel so that the R syntax and packages/libraries can be used. For the most part, all code can be copied and pasted directly into R Studio for running in a local environment.\n",
    "\n",
    "Please send any questions, comments, or use cases of interest to NSIDC User Services at nsidc@nsidc.org\n",
    "\n",
    "### Notebook Sections:\n",
    "_Click a link to jump to a section of interest_\n",
    "1. [Set up](#start) your JupyteR environment\n",
    "2. [Importing data](#import) from the web (HTTPS and FTP)\n",
    "3. Working with [NetCDF 4 (.nc)](#netcdf)\n",
    "4. Working with [HDF5 (.h5)](#hdf)\n",
    "5. Working with [HDF (.he5)](#he5)\n",
    "6. Working with [Binary (.bin)](#bin)\n",
    "7. Working with [multiple datasets](#mult)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "<a id='start'></a>\n",
    "## 1. **_Start Here_:** Steps for setting up your local JupyteR environment<h5>(Binder users can skip this step and move to the \"Load libraries\" cell)</h5>\n",
    "<b>Step 1.</b> Activate or create a conda environment with all of the packages required for this notebook (_note_: If you only need to run a specific section of the notebook, packages used are listed at the top of each section):<br />\n",
    "\n",
    "<i>conda create -n chooseyourenvname r-base r-essentials r-ncdf4 r-raster r-gdalUtils r-rgdal r-RCurl r-httr r-getPass r-sf r-rnaturalearth r-rnaturalearthdata</i>\n",
    "\n",
    "<b>Step 2.</b> Run <i>activate yourenvname</i>\n",
    "\n",
    "<b>Step 3.</b> Run <i>conda install -c conda-forge jupyterlab</i>\n",
    "\n",
    "<b>Step 4.</b> Run _jupyter lab_ to launch the notebook interface\n",
    "\n",
    "   <b><i>Notes:</i></b>\n",
    "   - Any R package you want to install with _conda_ should be written as \"r-<i>packagename</i>\"\n",
    "   - _r-essentials_ includes tidyverse, shiny, dplyr, and other useful packages so you don't need to install them individually.\n",
    "   - _rgdal_ can sometimes cause trouble for new users installing the package. If you do experience issues, learn more here: https://hdfeos.org/software/gdal.php"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "<b>Load libraries</b>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "libraries <- c(\"ncdf4\", \"rgdal\", \"gdalUtils\", \"raster\", \"RCurl\", \"R.utils\", \"httr\", \"getPass\", \"rnaturalearth\", \"rnaturalearthdata\", \"ggplot2\", \"sf\", \"sys\")\n",
    "invisible(lapply(libraries, library, character.only = TRUE))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<b>Note:</b> Each section starts with a cell setting a different working directory. You can skip these if you are working from a single directory."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "<a id='import'></a>\n",
    "## 2. **Importing NSIDC Datasets** from HTTPS and FTP Sources\n",
    "Package(s): **RCurl** & **httr** </b>Datasets: [**ATLAS/ICESat-2 L2A Global Geolocated Photon Data, Version 4**](https://nsidc.org/data/atl03), [**IMS Daily Northern Hemisphere Snow and Ice Analysis at 1 km, Version 1 (G02156)**](https://nsidc.org/data/G02156/versions/1)\n",
    "\n",
    "**LP DAAC** also has an R data access notebook: https://git.earthdata.nasa.gov/projects/LPDUR/repos/daac_data_download_r/browse/DAACDataDownload.R"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Set working directory and NASA Earthdata Authentication"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Check current working directory (usually the directory where Jupyter Lab was launched from)\n",
    "getwd()\n",
    "wd <- getwd()\n",
    "wd"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### [Create](https://nsidc.org/support/how/v0-programmatic-data-access-guide) a <b>.netrc</b> file that contains your NASA Earthdata login credientials manually and place it in your working directory.\n",
    "\n",
    "***Note:*** Creating a .netrc file from a Jupyter Notebook cell doesn't currently work as expected. You will need to manually create the file in the link above and run the code on your local computer to access data through the ****HTTPS Access method****."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set a path to .netrc file\n",
    "netrc_path <- file.path(wd,\".netrc\")\n",
    "path <- file.path(wd,\"data/\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### <u>HTTPS Access</u>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Single file:\n",
    "files <- \"https://n5eil01u.ecs.nsidc.org/ATLAS/ATL03.003/2018.11.23/ATL03_20181123000146_08450113_003_01.h5\"\n",
    "\n",
    "# A list of files:\n",
    "#files <- c(\"https://n5eil01u.ecs.nsidc.org/filepath/file.hdf\",\n",
    "#           \"https://n5eil01u.ecs.nsidc.org/filepath/file.hdf\")\n",
    "\n",
    "# Text file with list of links:\n",
    "#files <- readLines(\"C:/filepath/URL_file_list.txt\", warn = FALSE)\n",
    "\n",
    "# Loop through all files\n",
    "for (i in 1:length(files)) {\n",
    "  filename <-  tail(strsplit(files[i], '/')[[1]], n = 1) # Keep original filename\n",
    "\n",
    "  # Write file to disk (authenticating with netrc) using the current directory/filename\n",
    "  response <- GET(files[i], write_disk(file.path(path,filename), overwrite = TRUE), progress(),\n",
    "                  config(netrc = TRUE, netrc_file = netrc_path), set_cookies(\"LC\" = \"cookies\"))\n",
    "\n",
    "  # Check to see if file downloaded correctly\n",
    "  if (response$status_code == 200) {\n",
    "    print(sprintf(\"%s downloaded at %s\", filename, path))\n",
    "  } else {\n",
    "  print(sprintf(\"%s not downloaded. Verify that your username and password are correct in %s\", filename, netrc))\n",
    "  }\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### <u>NSIDC FTP Access</u>\n",
    "1. <b>Host Name:</b> sidads.colorado.edu <br />\n",
    "2. <b>User Name:</b> anonymous <br />\n",
    "3. <b>Password:</b> your e-mail address <br />\n",
    "4. <b>Directory:</b> /pub/DATASETS/XXXX <br />"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Getting Connected:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# This combination will be used throughout this section and shouldn't change\n",
    "usrpwd <- readline(prompt=\"Enter anonymous:youremailaddress \")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Your URL may change depending on what dataset you're interested in looking at. This will change in this section.\n",
    "url <- \"ftp://sidads.colorado.edu/pub/DATASETS/\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "filenames <- getURL(url, userpwd = usrpwd,\n",
    "             ftp.use.epsv = FALSE, dirlistonly = TRUE)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Uncomment the line below to get a list of sub-directories\n",
    "cat(filenames)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The data will be saved into your current working directory\n",
    "ims_url = \"ftp://sidads.colorado.edu/pub/DATASETS/NOAA/G02156/GIS/1km/2021/\"\n",
    "download.file(ims_url, destfile = \"ims2021001_1km_GIS_v1.3.tif.gz\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gunzip(\"ims2021001_1km_GIS_v1.3.tif.gz\", remove=FALSE)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now the unzipped file should be in your working directory and ready to use!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "<a id='netcdf'></a>\n",
    "## 3. Working with <b>NetCDF 4</b> files in R\n",
    "Package:<b> netcdf4 </b>Dataset:<b> [ATLAS/ICESat-2 L3B Daily and Monthly Gridded Sea Ice Freeboard, Version 2](https://nsidc.org/data/ATL20)</b>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set data path and filename\n",
    "path <- \"./data\"\n",
    "ncname <- \"/SNOWEX2020_EO_PLANE_2020Feb08_mosaicked_APLUW\"\n",
    "ncfname <- paste(path, ncname, \".nc\", sep=\"\")\n",
    "ncfname"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Open the netCDF file\n",
    "ncin <- nc_open(ncfname)\n",
    "\n",
    "# Only print if you want a very long output\n",
    "#print(ncin)\n",
    "\n",
    "# You can also use rgdal for netcdf files\n",
    "#ncin <- \"data/SNOWEX2020_EO_PLANE_2020Feb08_mosaicked_APLUW.nc\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get subdatasets so that you can grab the parameter slice you're interested in\n",
    "sds <- get_subdatasets(ncfname)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Convert your parameter slice to .tif and save to local directory\n",
    "gdal_translate(sds[2], dst_dataset = \"SNEX20_DEM_Save.tif\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Pull the new .tif raster back into R and examine the data\n",
    "dem <- raster(\"SNEX20_DEM_Save.tif\")\n",
    "coastlines <- ne_coastline(scale = 110, returnclass = c(\"sp\", \"sf\"))\n",
    "show(dem)\n",
    "plot(dem)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set up the projection for the data using WTK\n",
    "wkt_32613 <- 'PROJCS[\"WGS 84 / UTM zone 13N\",\n",
    "    GEOGCS[\"WGS 84\",\n",
    "        DATUM[\"WGS_1984\",\n",
    "            SPHEROID[\"WGS 84\",6378137,298.257223563,\n",
    "                AUTHORITY[\"EPSG\",\"7030\"]],\n",
    "            AUTHORITY[\"EPSG\",\"6326\"]],\n",
    "        PRIMEM[\"Greenwich\",0,\n",
    "            AUTHORITY[\"EPSG\",\"8901\"]],\n",
    "        UNIT[\"degree\",0.01745329251994328,\n",
    "            AUTHORITY[\"EPSG\",\"9122\"]],\n",
    "        AUTHORITY[\"EPSG\",\"4326\"]],\n",
    "    UNIT[\"metre\",1,\n",
    "        AUTHORITY[\"EPSG\",\"9001\"]],\n",
    "    PROJECTION[\"Transverse_Mercator\"],\n",
    "    PARAMETER[\"latitude_of_origin\",0],\n",
    "    PARAMETER[\"central_meridian\",-105],\n",
    "    PARAMETER[\"scale_factor\",0.9996],\n",
    "    PARAMETER[\"false_easting\",500000],\n",
    "    PARAMETER[\"false_northing\",0],\n",
    "    AUTHORITY[\"EPSG\",\"32613\"],\n",
    "    AXIS[\"Easting\",EAST],\n",
    "    AXIS[\"Northing\",NORTH]]'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "crs(dem) <- CRS(SRS_string=wkt_32613)\n",
    "extent(dem) <- c(-108.0000, 0.0000, -102.0000, 84.0000)\n",
    "show(dem)\n",
    "plot(dem)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Uncomment the line below to write the raster to your local directory\n",
    "#writeRaster(dem, \"SNEX20_DEM_EPSG32613_NK_QC_02102021.tif\", overwrite = TRUE)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "[**Back to top**](#top)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "<a id='hdf'></a>\n",
    "## 4. Working with <b> HDF5 </b> data in R\n",
    "Package:<b> rgdal </b>Dataset: **[ICESat-2 L3B Monthly Gridded Atmosphere, Version 3 (ATL16)](https://nsidc.org/data/ATL16)**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Convert a single file\n",
    "path <- \"./data\"\n",
    "is2_file <- 'ATL16_20190201001042_05280201_955a2_01.h5'\n",
    "filepath <- file.path(path, is2_file)\n",
    "print(filepath)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Be careful using this command as the list can be very long with unsubset HDF files.\n",
    "is2_sds <- get_subdatasets(filepath)\n",
    "list(is2_sds)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "od <- gdal_translate(is2_sds[9], dst_dataset = \"ATL16_global_od_WGS84_PanoplyFix.tif\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gaf <- raster(\"ATL16_global_od_WGS84_PanoplyFix.tif\")\n",
    "gaf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "wkt_wgs84 <- 'GEOGCS[\"WGS 84\",\n",
    "    DATUM[\"WGS_1984\",\n",
    "        SPHEROID[\"WGS 84\",6378137,298.257223563,\n",
    "            AUTHORITY[\"EPSG\",\"7030\"]],\n",
    "        AUTHORITY[\"EPSG\",\"6326\"]],\n",
    "    PRIMEM[\"Greenwich\",0,\n",
    "        AUTHORITY[\"EPSG\",\"8901\"]],\n",
    "    UNIT[\"degree\",0.01745329251994328,\n",
    "        AUTHORITY[\"EPSG\",\"9122\"]],\n",
    "    AUTHORITY[\"EPSG\",\"4326\"]]'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "coastlines <- ne_coastline(scale = 110, returnclass = c(\"sp\", \"sf\"))\n",
    "show(coastlines)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gaf.ext <- raster(\"ATL16_global_od_WGS84_PanoplyFix.tif\")\n",
    "crs(gaf.ext) <- CRS(\"+proj=longlat +datum=WGS84 +no_defs\")\n",
    "gaf.ext <- setExtent(gaf.ext, coastlines, keepres=FALSE)\n",
    "show(gaf.ext)\n",
    "#gaf.ext <- setExtent(gaf.ext, c(-180,180,-90,90), keepres = FALSE)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# As of 2/17/2021, ATL16 and ATL17 data needs to be flipped 180 degrees to properly plot the data\n",
    "gaf.ext.flip <- flip(gaf.ext, direction='y')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#jpeg(file=\"ATL16_Global_OD_Flip.jpeg\")\n",
    "\n",
    "plot(gaf.ext, \n",
    "      breaks = c(0.0,0.2, 0.4, 0.6, 0.8, 1.0, 3.0), \n",
    "      col = rainbow(5),\n",
    "      main=\"Global Column Optical Depth\",\n",
    "      legend=FALSE,\n",
    "      axes=TRUE)\n",
    "plot(coastlines, add=TRUE)\n",
    "\n",
    "#dev.off()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot(gaf.ext.flip, \n",
    "      breaks = c(0.0,0.2, 0.4, 0.6, 0.8, 1.0), \n",
    "      col = rainbow(5),\n",
    "      legend=TRUE,\n",
    "      axes=TRUE)\n",
    "plot(coastlines, add=TRUE)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Uncomment the following lines to save the plot\n",
    "# jpeg(file=\"ATL16_Global_OD_upsidedown.jpeg\")\n",
    "# plot(gaf.ext, \n",
    "#      breaks = c(0.2, 0.4, 0.6, 0.8, 1.0), \n",
    "#      col = rainbow(5),\n",
    "#      main=file,\n",
    "#      legend=FALSE,\n",
    "#      axes=TRUE)\n",
    "# plot(coastlines, add=TRUE)\n",
    "# dev.off()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Uncomment the following line if you want to save the final GeoTiff\n",
    "#writeRaster(gaf.ext.flip, 'ATL16_global_column_od_flip.tif', options=c('TFW=YES', 'XML=YES'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[**Back to top**](#top)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "<a id='he5'></a>\n",
    "## 5. Working with <b>.he5 </b>data in R\n",
    "Package:<b> rgdal </b>Dataset: **[AMSR-E/AMSR2 Unified L3 Daily 12.5 km Brightness Temperatures, Sea Ice Concentration, Motion & Snow Depth Polar Grids](https://nsidc.org/data/AU_SI12/versions/1)**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Convert a single file\n",
    "amsr_file <- \"./data/AMSR_U2_L3_SeaIce12km_B04_20210219.he5\"\n",
    "print(amsr_file)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# View subdatasets within the HDF file\n",
    "amsr_sds <- get_subdatasets(amsr_file)\n",
    "#amsr_sds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# convert your data slice to .tif and save new file to correct CRS and extent\n",
    "gdal_translate(amsr_sds[59], dst_dataset = \"SI_12km_SH_ICECON_DAY_20210219.tif\")\n",
    "amsr_gaf <- raster(\"SI_12km_SH_ICECON_DAY_20210219.tif\")\n",
    "show(amsr_gaf)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "wkt_south <- 'PROJCS[\"WGS 84 / NSIDC Sea Ice Polar Stereographic South\",\n",
    "    GEOGCS[\"WGS 84\",\n",
    "        DATUM[\"WGS_1984\",\n",
    "            SPHEROID[\"WGS 84\",6378137,298.257223563,\n",
    "                AUTHORITY[\"EPSG\",\"7030\"]],\n",
    "            AUTHORITY[\"EPSG\",\"6326\"]],\n",
    "        PRIMEM[\"Greenwich\",0,\n",
    "            AUTHORITY[\"EPSG\",\"8901\"]],\n",
    "        UNIT[\"degree\",0.0174532925199433,\n",
    "            AUTHORITY[\"EPSG\",\"9122\"]],\n",
    "        AUTHORITY[\"EPSG\",\"4326\"]],\n",
    "    PROJECTION[\"Polar_Stereographic\"],\n",
    "    PARAMETER[\"latitude_of_origin\",-70],\n",
    "    PARAMETER[\"central_meridian\",0],\n",
    "    PARAMETER[\"scale_factor\",1],\n",
    "    PARAMETER[\"false_easting\",0],\n",
    "    PARAMETER[\"false_northing\",0],\n",
    "    UNIT[\"metre\",1,\n",
    "        AUTHORITY[\"EPSG\",\"9001\"]],\n",
    "    AXIS[\"X\",EAST],\n",
    "    AXIS[\"Y\",NORTH],\n",
    "    AUTHORITY[\"EPSG\",\"3976\"]]'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "crs(amsr_gaf) <- CRS(SRS_string=wkt_south)\n",
    "extent(amsr_gaf) <- c(-3950000, 3950000, -3950000, 4350000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "show(amsr_gaf)\n",
    "plot(amsr_gaf)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "amsr_gaf_df <- as.data.frame(amsr_gaf, xy=TRUE)\n",
    "str(amsr_gaf_df)\n",
    "\n",
    "min <- minValue(amsr_gaf)\n",
    "max <- maxValue(amsr_gaf)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ggplot() +\n",
    "    geom_raster(data = amsr_gaf_df , aes(x = x, y = y, fill = SI_12km_SH_ICECON_DAY_20210219)) +\n",
    "    scale_fill_viridis_c() +\n",
    "    coord_quickmap()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#writeRaster(amsr_gaf, \"SI_12km_SH_ICECON_DAY_020821_extent.tif\", overwrite = TRUE)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[**Back to top**](#top)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "<a id='bin'></a>\n",
    "## 6. Working with <b>Binary (.bin) </b>data in R\n",
    "Package:<b> rgdal </b>Dataset: **[Sea Ice Concentrations from Nimbus-7 SMMR and DMSP SSM/I-SSMIS Passive Microwave Data, Version 1 (NSIDC-0051)](https://nsidc.org/data/NSIDC-0051/versions/1)**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Note:** Before you can convert Binary data to a raster, you must first create a _header file (.hdr)_ and place it in your working directory. NSIDC has a _How To_ article that includes the necessary header file information for the Northern and Southern Hemisphere versions of this example data product: https://nsidc.org/support/how/how-do-i-convert-nsidc-0051-sea-ice-concentration-data-binary-geotiff"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nt_file <- \"./data/nt_20201031_f17_v1.1_s.bin\"\n",
    "nt_file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gdal_translate(nt_file, \n",
    "               of = \"GTiff\",\n",
    "               dst_dataset = \"nt_test.tif\", \n",
    "               a_srs = '+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',\n",
    "               projwin_srs = '+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',\n",
    "               a_nodata = 255, \n",
    "               a_ullr = c(-3950000.0, 4350000.0, 3950000.0, -3950000.0),\n",
    "               scale = c(0, 255),\n",
    "               output_Raster = TRUE, verbose = FALSE);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nt_raster <- raster(\"t_20201031_f17_v1.1_s.tif\")\n",
    "plot(nt_raster)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[**Back to top**](#top)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "<a id='mult'></a>\n",
    "## 7. Working with **Multiple Datasets** in R\n",
    "Package(s): **rgdal, ggplot2** </b>Dataset: **[MEaSUREs inSAR-Based Antarctic Ice Velocity Map, Version 2 (NSIDC-0484)](https://nsidc.org/data/NSIDC-0484/versions/2), [GLIMS Glacier Database](http://glims.colorado.edu/glacierdata/)**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "path <- \"./data\"\n",
    "file_path <- file.path(path, \"glims_polygons.shp\")\n",
    "file_path\n",
    "pine_island <- st_read(file_path)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ggplot() + \n",
    "  geom_sf(data = pine_island, size = 1, color = \"black\", fill = \"cyan1\") + \n",
    "  ggtitle(\"Pine Island Glacier, Antarctica\") + \n",
    "  coord_sf()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "meas_file <- \"antarctica_ice_velocity_450m_v2.nc\"\n",
    "meas_path <- file.path(path, meas_file)\n",
    "meas_path"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# View subdatasets within the HDF file\n",
    "meas_sds <- get_subdatasets(meas_path)\n",
    "meas_sds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gdal_translate(meas_sds[3], dst_dataset = \"VX.tif\") #, a_srs = \"EPSG:3412\", s_srs = \"EPSG:3412\") #, projwin_srs = \"EPSG:3976\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "vx <- raster(\"VX.tif\")\n",
    "show(vx)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pine_island_prj <- st_transform(pine_island, crs = st_crs(vx))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot(vx)\n",
    "plot(pine_island_prj, add=TRUE)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "vx_crop <- crop(vx, pine_island_prj)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pine_island_prj_df <- fortify(pine_island_prj)\n",
    "vx_crop_df <- as.data.frame(vx_crop, xy = TRUE) \n",
    "str(vx_crop_df)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ggplot() +\n",
    "  geom_raster(data = vx_crop_df, \n",
    "              aes(x = x, y = y, \n",
    "                  fill = VX)) + \n",
    "  scale_fill_viridis_c() +  \n",
    "  scale_alpha(range = c(0.15, 0.65), guide = \"none\") +  \n",
    "  ggtitle(\"Pine Island Ice Velocity\") +\n",
    "  coord_quickmap()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "[**Back to the top**](#top)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "4.0.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}