Georef_AMSR2_SWE_PolarGrid_Data_Recipe.ipynb

{"cells":[{"cell_type":"markdown","metadata":{},"source":["  <br>\n","   This code takes GHRC LANCE AMSR2 Snow Water Equivalent Data <br>\n","   formatted in 25km HDF-EOS5 files, georeferences it to <br>\n","   NSIDC Polar Grids, and creates four GeoTIFF files.<br>\n","  <br>\n","  Inputs: GHRC NRT AMSR2 Daily L3 Global Snow Water Equivalent EASE-Grids<br>\n","  Outputs: Four GeoTIFF files of North and South Pole Snow Water Equivalent<br>\n","          (SWE) and flags.<br>\n","   <br>\n","   Authors: Kel Markert (1), Amanda Weigel(2)<br>\n","<br>\n","   (1) Earth System Science Center (ESSC)    <br>\n","   (2)Information and Technology Systems Center (ITSC)<br>\n","   University of Alabama in Huntsville<br>\n","<br>\n","  Data Edited: Date: April 07, 2017<br>\n","<br>"]},{"cell_type":"code","execution_count":2,"metadata":{},"outputs":[{"name":"stdout","output_type":"stream","text":["Requirement already satisfied: numpy in /home/davis_dev/anaconda3/lib/python3.9/site-packages (1.21.5)\n","Requirement already satisfied: tables in /home/davis_dev/anaconda3/lib/python3.9/site-packages (3.6.1)\n","Requirement already satisfied: numpy>=1.9.3 in /home/davis_dev/anaconda3/lib/python3.9/site-packages (from tables) (1.21.5)\n","Requirement already satisfied: numexpr>=2.6.2 in /home/davis_dev/anaconda3/lib/python3.9/site-packages (from tables) (2.8.1)\n","Requirement already satisfied: packaging in /home/davis_dev/anaconda3/lib/python3.9/site-packages (from numexpr>=2.6.2->tables) (21.3)\n","Requirement already satisfied: pyparsing!=3.0.5,>=2.0.2 in /home/davis_dev/anaconda3/lib/python3.9/site-packages (from packaging->numexpr>=2.6.2->tables) (3.0.4)\n","Requirement already satisfied: pyproj in /home/davis_dev/anaconda3/lib/python3.9/site-packages (3.3.1)\n","Requirement already satisfied: certifi in /home/davis_dev/anaconda3/lib/python3.9/site-packages (from pyproj) (2021.10.8)\n"]},{"name":"stderr","output_type":"stream","text":["/home/davis_dev/anaconda3/envs/fresh/lib/python3.8/site-packages/tables/attributeset.py:291: DataTypeWarning: Unsupported type for attribute 'DIMENSION_LIST' in node 'Flags_NorthernDaily'. Offending HDF5 class: 9\n","  value = self._g_getattr(self._v_node, name)\n","/home/davis_dev/anaconda3/envs/fresh/lib/python3.8/site-packages/tables/attributeset.py:291: DataTypeWarning: Unsupported type for attribute 'DIMENSION_LIST' in node 'SWE_NorthernDaily'. Offending HDF5 class: 9\n","  value = self._g_getattr(self._v_node, name)\n","/home/davis_dev/anaconda3/envs/fresh/lib/python3.8/site-packages/tables/attributeset.py:291: DataTypeWarning: Unsupported type for attribute 'DIMENSION_LIST' in node 'Flags_SouthernDaily'. Offending HDF5 class: 9\n","  value = self._g_getattr(self._v_node, name)\n","/home/davis_dev/anaconda3/envs/fresh/lib/python3.8/site-packages/tables/attributeset.py:291: DataTypeWarning: Unsupported type for attribute 'DIMENSION_LIST' in node 'SWE_SouthernDaily'. Offending HDF5 class: 9\n","  value = self._g_getattr(self._v_node, name)\n","/home/davis_dev/anaconda3/envs/fresh/lib/python3.8/site-packages/pyproj/crs/crs.py:130: FutureWarning: '+init=<authority>:<code>' syntax is deprecated. '<authority>:<code>' is the preferred initialization method. When making the change, be mindful of axis order changes: https://pyproj4.github.io/pyproj/stable/gotchas.html#axis-order-changes-in-proj-6\n","  in_crs_string = _prepare_from_proj_string(in_crs_string)\n","/home/davis_dev/anaconda3/envs/fresh/lib/python3.8/site-packages/pyproj/crs/crs.py:130: FutureWarning: '+init=<authority>:<code>' syntax is deprecated. '<authority>:<code>' is the preferred initialization method. When making the change, be mindful of axis order changes: https://pyproj4.github.io/pyproj/stable/gotchas.html#axis-order-changes-in-proj-6\n","  in_crs_string = _prepare_from_proj_string(in_crs_string)\n","/tmp/ipykernel_2718/1342170419.py:122: DeprecationWarning: This function is deprecated. See: https://pyproj4.github.io/pyproj/stable/gotchas.html#upgrading-to-pyproj-2-from-pyproj-1\n","  ul_x,ul_y = transform(inProj,outProj,ullon,ullat)\n","ERROR 1: unhandled axis direction: UNKNOWN\n","/home/davis_dev/anaconda3/envs/fresh/lib/python3.8/site-packages/pyproj/crs/crs.py:130: FutureWarning: '+init=<authority>:<code>' syntax is deprecated. '<authority>:<code>' is the preferred initialization method. When making the change, be mindful of axis order changes: https://pyproj4.github.io/pyproj/stable/gotchas.html#axis-order-changes-in-proj-6\n","  in_crs_string = _prepare_from_proj_string(in_crs_string)\n","/tmp/ipykernel_2718/1342170419.py:122: DeprecationWarning: This function is deprecated. See: https://pyproj4.github.io/pyproj/stable/gotchas.html#upgrading-to-pyproj-2-from-pyproj-1\n","  ul_x,ul_y = transform(inProj,outProj,ullon,ullat)\n","ERROR 1: unhandled axis direction: UNKNOWN\n","/home/davis_dev/anaconda3/envs/fresh/lib/python3.8/site-packages/pyproj/crs/crs.py:130: FutureWarning: '+init=<authority>:<code>' syntax is deprecated. '<authority>:<code>' is the preferred initialization method. When making the change, be mindful of axis order changes: https://pyproj4.github.io/pyproj/stable/gotchas.html#axis-order-changes-in-proj-6\n","  in_crs_string = _prepare_from_proj_string(in_crs_string)\n","/tmp/ipykernel_2718/1342170419.py:122: DeprecationWarning: This function is deprecated. See: https://pyproj4.github.io/pyproj/stable/gotchas.html#upgrading-to-pyproj-2-from-pyproj-1\n","  ul_x,ul_y = transform(inProj,outProj,ullon,ullat)\n","ERROR 1: unhandled axis direction: UNKNOWN\n","/home/davis_dev/anaconda3/envs/fresh/lib/python3.8/site-packages/pyproj/crs/crs.py:130: FutureWarning: '+init=<authority>:<code>' syntax is deprecated. '<authority>:<code>' is the preferred initialization method. When making the change, be mindful of axis order changes: https://pyproj4.github.io/pyproj/stable/gotchas.html#axis-order-changes-in-proj-6\n","  in_crs_string = _prepare_from_proj_string(in_crs_string)\n","/tmp/ipykernel_2718/1342170419.py:122: DeprecationWarning: This function is deprecated. See: https://pyproj4.github.io/pyproj/stable/gotchas.html#upgrading-to-pyproj-2-from-pyproj-1\n","  ul_x,ul_y = transform(inProj,outProj,ullon,ullat)\n","ERROR 1: unhandled axis direction: UNKNOWN\n"]}],"source":["! pip install numpy\n","! pip install tables\n","! pip install pyproj\n","\n","import numpy\n","from osgeo import gdal, osr\n","from osgeo.gdalnumeric import *  \n","from osgeo.gdalconst import *\n","import tables\n","from pyproj import Proj, transform\n","\n","# Path to input file\n","\n","infile = 'test_files/georef/AMSR_U2_L3_DailySnow_P02_20220823.he5'\n","\n","# Path to output file directory\n","\n","outfolder = 'test_files/out/' # must have '/' at the end\n","\n","# Open HDF5 file\n","\n","h5file = tables.open_file(infile)\n","\n","# Define metadata path to data layer of interest within file\n","\n","north_grid_flags = h5file.get_node('/HDFEOS/GRIDS/Northern Hemisphere/Data Fields/Flags_NorthernDaily').read()\n","north_grid_SWE = h5file.get_node('/HDFEOS/GRIDS/Northern Hemisphere/Data Fields/SWE_NorthernDaily').read()\n","\n","south_grid_flags = h5file.get_node('/HDFEOS/GRIDS/Southern Hemisphere/Data Fields/Flags_SouthernDaily').read()\n","south_grid_SWE = h5file.get_node('/HDFEOS/GRIDS/Southern Hemisphere/Data Fields/SWE_SouthernDaily').read()\n","\n","# Close HDF5 file\n","\n","h5file.close()\n","\n","# Make sure data layers are listed in the same order as they appear in file\n","\n","file_dict = {'NQAF':north_grid_flags,'NSWE':north_grid_SWE,'SQAF':south_grid_flags,'SSWE':south_grid_SWE}\n","             \n","products = ['NQAF','NSWE','SQAF','SSWE',]\n","\n","# Set output projections\n","\n","northProj = '''PROJCS[\"NSIDC EASE-Grid North\",\n","    GEOGCS[\"Unspecified datum based upon the International 1924 Authalic Sphere\",\n","        DATUM[\"Not_specified_based_on_International_1924_Authalic_Sphere\",\n","            SPHEROID[\"International 1924 Authalic Sphere\",6371228,0,\n","                AUTHORITY[\"EPSG\",\"7057\"]],\n","            AUTHORITY[\"EPSG\",\"6053\"]],\n","        PRIMEM[\"Greenwich\",0,\n","            AUTHORITY[\"EPSG\",\"8901\"]],\n","        UNIT[\"degree\",0.01745329251994328,\n","            AUTHORITY[\"EPSG\",\"9122\"]],\n","        AUTHORITY[\"EPSG\",\"4053\"]],\n","    UNIT[\"metre\",1,\n","        AUTHORITY[\"EPSG\",\"9001\"]],\n","    PROJECTION[\"Lambert_Azimuthal_Equal_Area\"],\n","    PARAMETER[\"latitude_of_center\",90],\n","    PARAMETER[\"longitude_of_center\",0],\n","    PARAMETER[\"false_easting\",0],\n","    PARAMETER[\"false_northing\",0],\n","    AUTHORITY[\"EPSG\",\"3408\"],\n","    AXIS[\"X\",UNKNOWN],\n","    AXIS[\"Y\",UNKNOWN]]'''    #NSIDC North Pole Ease-Grid\n","    \n","    \n","southProj = '''PROJCS[\"NSIDC EASE-Grid South\",\n","    GEOGCS[\"Unspecified datum based upon the International 1924 Authalic Sphere\",\n","        DATUM[\"Not_specified_based_on_International_1924_Authalic_Sphere\",\n","            SPHEROID[\"International 1924 Authalic Sphere\",6371228,0,\n","                AUTHORITY[\"EPSG\",\"7057\"]],\n","            AUTHORITY[\"EPSG\",\"6053\"]],\n","        PRIMEM[\"Greenwich\",0,\n","            AUTHORITY[\"EPSG\",\"8901\"]],\n","        UNIT[\"degree\",0.01745329251994328,\n","            AUTHORITY[\"EPSG\",\"9122\"]],\n","        AUTHORITY[\"EPSG\",\"4053\"]],\n","    UNIT[\"metre\",1,\n","        AUTHORITY[\"EPSG\",\"9001\"]],\n","    PROJECTION[\"Lambert_Azimuthal_Equal_Area\"],\n","    PARAMETER[\"latitude_of_center\",-90],\n","    PARAMETER[\"longitude_of_center\",0],\n","    PARAMETER[\"false_easting\",0],\n","    PARAMETER[\"false_northing\",0],\n","    AUTHORITY[\"EPSG\",\"3409\"],\n","    AXIS[\"X\",UNKNOWN],\n","    AXIS[\"Y\",UNKNOWN]]'''    #NSIDC South Pole EASE-Grid\n","    \n","spat_res  = 25000 #pixel spatial resolution in meters\n","inProj = Proj(init='epsg:4326')\n","\n","# Load GDAL GeoTIFF dataset driver for dataset format\n","\n","drv = gdal.GetDriverByName(\"GTiff\")\n","\n","# Looping over 4 hdf data layers within each data file\n","\n","for i in range(len(products)):\n","    #grab data from file_dict (list of data layer arrays)\n","    indata = file_dict[products[i]]\n","    #output file path using names in file_dict  \n","    outfile = outfolder+products[i]+'.tif'\n","    \n","    #Checks if North or South NSIDC Ease Grid pPojection applies\n","    #set upper left lon and lat extent dimensions for polar data\n","    \n","    #Northern Hemisphere\n","    if products[i][0] == 'N':\n","        ullon = -135 #this is an approximate value\n","        ullat = -86.5 #this is an approximate value\n","        dest_wkt = northProj\n","        outProj = Proj(init='epsg:3408')\n","    \n","    #Southern Hemisphere\n","    else:\n","        ullon = -45 #this is an approximate value\n","        ullat = 86.5 #this is an approximate value\n","        dest_wkt = southProj\n","        outProj = Proj(init='epsg:3409')\n","    \n","    #Converts coordinates from WGS 1984 (entered above) to polar grid coordinates   \n","    ul_x,ul_y = transform(inProj,outProj,ullon,ullat)\n","    \n","    # Create blank ouput dataset.  Set (outfile, x-dimension,,y-dimension, number of bands, data type).\n","    dsOut = drv.Create(outfile, indata.shape[1],  indata.shape[0], 1, gdal.GDT_Byte)\n","    \n","    # Specify and set output file's geotransform. \n","    # geotransform = [upper left longitude, x-dimension resolution, x skew, upper left latitude, y skew, y-dimension resolution (usually negative)\n","    gt = [ul_x, spat_res, 0, ul_y, 0, -spat_res]\n","    dsOut.SetGeoTransform(gt)\n","    \n","    # Set output dataset projection\n","    dsOut.SetProjection(dest_wkt)\n","    \n","    # Write gridded swath to the output file\n","    bandOut=dsOut.GetRasterBand(1)\n","    BandWriteArray(bandOut, indata)\n","    \n","    # Release output dataset from memory\n","    dsOut = None\n","    bandOut = None \n","    \n","#once the file GeoTIFFs have been created, pull into ArcMap and define the projection \n","#as NSIDC EASE Grid North and NSIDC EASE Grid South\n"]}],"metadata":{"kernelspec":{"display_name":"Python 3.8.13 ('fresh')","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.8.13"},"vscode":{"interpreter":{"hash":"41b8fb484b64eab2e0ad5b1c1d559b9e124087993d39def2a4dca1e0bb9d0d7c"}}},"nbformat":4,"nbformat_minor":2}