UsefulAssets.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "640189da",
   "metadata": {},
   "source": [
    "# Melt-Pond Classification Functions\n",
    "\n",
    "In this script, we provide a few important functions used for analyzing melt pond fractions from the sea-ice classification algorithm using Sentinel-2 imagery. The main functions are as follows:\n",
    "\n",
    "1. Assigning projection information from classified HDF files.\n",
    "2. Mosaicing .tif and .h5 layers.\n",
    "3. Calculating Melt Pond Fraction and Sea-Ice concentration."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ac4f8a96",
   "metadata": {},
   "outputs": [],
   "source": [
    "from osgeo import gdal\n",
    "import os\n",
    "\n",
    "import rasterio\n",
    "from rasterio import mask\n",
    "from pylab import *\n",
    "import numpy as np\n",
    "import sys\n",
    "import pandas as pd\n",
    "import h5py\n",
    "from PIL import Image\n",
    "import csv"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "71a5d687",
   "metadata": {},
   "source": [
    "## 1. Assigning projection information for classified HDF files\n",
    "\n",
    "The classified HDF files lack georeferencing in their native format. This script uses GDAL to georeference the files to the tiff format. This will allow for easier comparisons with the satellite imagery. To properly run the tool, users will have to specify the spatial extent of the imagery.\n",
    "\n",
    "1. Upper Left Longitude\n",
    "2. Upper Left Latitude\n",
    "3. Lower Right Longitude\n",
    "4. Lower Right Latitude\n",
    "\n",
    "To find extent coordinates, users may view properties of the image in their preferred GIS. \n",
    "\n",
    "In this example, we convert the HDF to a georefferenced .tif file. Users may change the output format based on their preference. \n",
    "\n",
    "**Note: Users should make sure to assign the correct EPSG code for proper projection**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "745801ae",
   "metadata": {},
   "outputs": [],
   "source": [
    "def GDAL_georef(directory, HDF_File, extent, output):\n",
    "    os.chdir(directory)\n",
    "    rasterFiles = os.listdir(os.getcwd())\n",
    "\n",
    "    # We will convert the HDF to a georefferenced tif\n",
    "    fileExtension = \".tif\"\n",
    "\n",
    "    ### Open the HDF File ###\n",
    "    hdflayer = gdal.Open(rasterFiles[0], gdal.GA_ReadOnly)\n",
    "    rlayer = gdal.Open(HDF_File, gdal.GA_ReadOnly)\n",
    "\n",
    "    # naming output objects\n",
    "    outputRaster = directory + output + fileExtension\n",
    "\n",
    "    # Georeferencing with EPSG:32620 (WGS 84 UTM zone 20N)\n",
    "    ### NOTE: Make sure you assign an appropriate EPSG code for your projection region !!! ###\n",
    "    EPSG = \"-a_srs EPSG:32620\"\n",
    "\n",
    "    translateOptionText = EPSG + \" -a_ullr \" + extent\n",
    "    translateoptions = gdal.TranslateOptions(gdal.ParseCommandLine(translateOptionText))\n",
    "    gdal.Translate(outputRaster, rlayer, options = translateoptions)\n",
    "    hdflayer = None\n",
    "    rlayer = None\n",
    "    return \"Classification georefferenced with EPSG 32620\"\n",
    "\n",
    "# These are what we need to pass into the function\n",
    "# Please make sure to change parameters to you local machine.\n",
    "directory = 'INSERT FULL DIRECTORY PATH HERE'\n",
    "HDF_File = 'T20XNS_20220804T192921_classification.h5' #Example HDF output\n",
    "\n",
    "# Extent is defined by GDAL bounding coordinates i.e. ullon ullat lrlon lrlat\n",
    "extent = '499980 9200040 609780 9090240' #Provided as an example\n",
    "output = 'GeorefferencedClassification'\n",
    "\n",
    "#How to call the function\n",
    "GDAL_options = GDAL_georef(directory, HDF_File, extent, output)\n",
    "print(GDAL_options)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "99128eea",
   "metadata": {},
   "source": [
    "## 2. Mosaic function for different classification/imagery types\n",
    "\n",
    "Classification HDF outputs are processed tile-by-tile based on individual imagery files. Users may want to view results as one continuous layer. With this function, users may create mosaics that patch classification results or satellite imagery as one layer. The function searches files with specific file extensions within a specified directory and creates a mosaic of all matching files. In this notebook we provide an example for mosaicking '.jp2' files. However, users may mosaic HDF files, NetCDF, .tiff, .jpeg, etc. \n",
    "\n",
    "**Note: Items in the mosaic function require native georefferencing. The Sentinel-2 imagery contains this from the base download. However, classified HDF's require georefferencing. Please refer to the previous function for details.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d7caea4f",
   "metadata": {},
   "outputs": [],
   "source": [
    "def Mosaic(base_direc, output_tif):\n",
    "    res = [] #Empty list to store files\n",
    "    # Iterate directory (based on base_direc)\n",
    "    for file in os.listdir(base_direc):\n",
    "       #Check only .jp2 files (Could also change to .h5, .nc, or even combine .tif formats)\n",
    "       if file.endswith('.jp2'):\n",
    "            res.append(base_direc + file)\n",
    "\n",
    "    output_src = base_direc + output_tif\n",
    "    g = gdal.Warp(output_src, res, format = \"GTiff\", options = [\"COMPRESS=LZW\", \"TILED=YES\"])\n",
    "    g = None #Close file and flush to disk\n",
    "\n",
    "    return \"Bands Merged\"\n",
    "\n",
    "base_direc = 'INSERT FULL PATH DIRECTORY HERE'\n",
    "output_tif = 'MethodTest.tif' #Provided as an example\n",
    "\n",
    "# How to call the function\n",
    "print(Mosaic(base_direc, output_tif))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c90436dc",
   "metadata": {},
   "source": [
    "## 3. Calculating Melt Pond Fraction (MPF) and Sea-Ice Concentration (SIC)\n",
    "\n",
    "This function is based off the \"S2_Results_Analysis\" script but might be more efficient for users who have processed their results in subsetted temporal ranges i.e. weekly, monthly, etc. Users can pass in text files created from the classification script and make the required calculations for Melt Pond Fraction and Sea-Ice Concentration with exports to .csv format. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bfeec92e",
   "metadata": {},
   "outputs": [],
   "source": [
    "def MPF(output_loc):\n",
    "    data = pd.read_csv(output_loc, delimiter='\\t',names=['Tile','IMG_Pixel','Border','Ice','OW','MP','Other'])\n",
    "\n",
    "    # Calculate MPF and SIC:\n",
    "    data['MPF'] = data.MP/(data.MP + data.Ice)\n",
    "    data['SIC'] = (data.Ice + data.MP)/(data.Ice + data.MP + data.OW)\n",
    "    data['Surface'] = data.IMG_Pixel - data.Border\n",
    "\n",
    "    # We don't want to calculate melt pond fraction for tiles with a sea-ice concentration less than 15%.\n",
    "    # For those values, we convert to nan.\n",
    "\n",
    "    for b in data[data.SIC<.15].index.values:\n",
    "        data.at[b, 'MPF'] = np.nan\n",
    "\n",
    "    # Day of Year (For Datetime objects)\n",
    "    # We parse out the datetime objects based on the file naming conventions\n",
    "    # The \"tile\" dataframe includes the tile and the \"date\" the image was obtained.\n",
    "    # By parsing out the position of the date values, we can extract datetime objects.\n",
    "    data['DOY']=[datetime.datetime.strptime(a,'%Y%m%d').timetuple().tm_yday for a in [a[7:15] for a in data.Tile]]\n",
    "\n",
    "    return data\n",
    "\n",
    "# Can export MPF objects to csv for instance:\n",
    "classification_loc = 'ENTER FULL TEXT FILE PATH' # i.e. '~/directory/output.txt'\n",
    "MPF(classification_loc).to_csv('ENTER FILE NAME YOU WISH TO SAVE CSV WITH FILE EXTENSION')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}