read_and_map_viirs.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "read_and_map_viirs.ipynb",
      "provenance": [],
      "collapsed_sections": []
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "im8zRbMslQLL"
      },
      "source": [
        "**Module:** read_and_map_viirs.ipynb\n",
        "\n",
        "**Disclaimer**: The code is for demonstration purposes only. Users are responsible to check for accuracy and revise to fit their objective.\n",
        "\n",
        "**Organization**: NASA ARSET\n",
        "\n",
        "**Author**: Justin Roberts-Pierel and Pawan Gupta, 2015.\n",
        "\n",
        "**Modified to work with netCDF** : Vikalp Mishra, 2019 \n",
        "\n",
        "**Modified to work with VIIRS data**: Aavash Thapa, 2020\n",
        "\n",
        "**Purpose**: To extract variables from a VIIRS netCDF4 file and create and save a map of the data"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Y7U02gkECz6c"
      },
      "source": [
        "#Mount drive to save files there\n",
        "#clone the repository to access files from there\n",
        "#pull the latest\n",
        "from google.colab import drive\n",
        "drive.mount('/content/drive', force_remount=True)\n",
        "! git clone https://github.com/NASAARSET/VIIRS_NASA.git\n",
        "! git -C VIIRS_NASA/ pull"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "4HOszY1RByrT"
      },
      "source": [
        "! pip install netCDF4\n",
        "import numpy as np\n",
        "import pandas as pd\n",
        "import sys\n",
        "from netCDF4 import Dataset\n",
        "import matplotlib.pyplot as plt\n",
        "#Colab requires specific installation of cartopy\n",
        "!apt-get -qq install python-cartopy python3-cartopy;\n",
        "!pip uninstall -y shapely;    # cartopy and shapely aren't friends (early 2020)\n",
        "!pip install shapely --no-binary shapely;\n",
        "import cartopy.crs as ccrs\n",
        "import cartopy.feature as cfeature\n",
        "from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter\n",
        "from matplotlib.axes import Axes\n",
        "from cartopy.mpl.geoaxes import GeoAxes\n",
        "GeoAxes._pcolormesh_patched = Axes.pcolormesh\n",
        "from textwrap import wrap"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "HNl5bq70Coha"
      },
      "source": [
        "#!/usr/bin/python\n",
        "try:\n",
        "    fileList = open('VIIRS_NASA/fileList.txt', 'r')\n",
        "except:\n",
        "    print('Did not find a text file containing file names (perhaps name does not match)')\n",
        "    sys.exit()\n",
        "\n",
        "#loops through all files listed in the text file\n",
        "for FILE_NAME in fileList:\n",
        "    FILE_NAME=FILE_NAME.strip()\n",
        "    #change 'raw_input' to 'input' if an error is shown about the input\n",
        "    user_input=input('\\nWould you like to process\\n' + FILE_NAME + '\\n\\n(Y/N)')\n",
        "    if(user_input == 'N' or user_input == 'n'):\n",
        "        print('Skipping...')\n",
        "        continue\n",
        "    else:\n",
        "        file = Dataset('VIIRS_NASA/'+ FILE_NAME, 'r')\n",
        "# read the data\n",
        "        ds=file\n",
        "        #grp='PRODUCT'        \n",
        "        lat= ds.variables['Latitude'][:][:]\n",
        "        lon= ds.variables['Longitude'][:][:]\n",
        "        if 'AERDB' in FILE_NAME:\n",
        "            \n",
        "            #The user has a choice of 5 sds variable and has to input a number to choose.\n",
        "            #The loop keeps repeating until the user inputs a value between 1-5 inclusive.\n",
        "            while  True:\n",
        "              choice = input(\"\"\"Pick the number with the corresponding sds variable of your choice: \n",
        "              1) Aerosol_Optical_Thickness_550_Land\n",
        "              2) Aerosol_Optical_Thickness_550_Land_Ocean_Best_Estimate\n",
        "              3) Aerosol_Optical_Thickness_QA_Flag_Land\n",
        "              4) Aerosol_Type_Land_Ocean\n",
        "              5) Angstrom_Exponent_Land_Ocean_Best_Estimate \"\"\")\n",
        "              \n",
        "              if choice in ['1', '2', '3', '4', '5']:\n",
        "                break\n",
        "              else:\n",
        "                print(\"Please input a valid response!\")\n",
        "\n",
        "\n",
        "            if choice == '1':\n",
        "              sds_name='Aerosol_Optical_Thickness_550_Land'\n",
        "            elif choice =='2':\n",
        "              sds_name='Aerosol_Optical_Thickness_550_Land_Ocean_Best_Estimate'\n",
        "            elif choice =='3':\n",
        "              sds_name='Aerosol_Optical_Thickness_QA_Flag_Land'\n",
        "            elif choice =='4':\n",
        "              sds_name='Aerosol_Type_Land_Ocean'\n",
        "            elif choice =='5':\n",
        "              sds_name='Angstrom_Exponent_Land_Ocean_Best_Estimate'\n",
        "\n",
        "        data= ds.variables[sds_name]\n",
        "        map_label = sds_name\n",
        "        map_label = map_label.replace('_', ' ')\n",
        "        map_label = '\\n'.join(wrap(map_label, 40)) \n",
        "        #get necessary attributes \n",
        "        fv=data._FillValue\n",
        "\n",
        "        #get lat and lon information \n",
        "        min_lat=np.min(lat)\n",
        "        max_lat=np.max(lat)\n",
        "        min_lon=np.min(lon)\n",
        "        max_lon=np.max(lon)\n",
        "        \n",
        "        # set map labels\n",
        "        #map_label = data.units\n",
        "        map_title = data.long_name\n",
        "        #print(data.units)\n",
        "    \n",
        "        #get the data as an array and mask fill/missing values\n",
        "        dataArray=np.array(data[:][:])\n",
        "        dataArray = np.multiply(dataArray, 1.0)\n",
        "        fv = fv*1.0\n",
        "        dataArray[dataArray==fv]=np.nan\n",
        "        data=dataArray\n",
        "        \n",
        "        \n",
        "        #get statistics about data\n",
        "        average=np.nanmean(dataArray)\n",
        "        stdev=np.nanstd(dataArray)\n",
        "        median=np.nanmedian(dataArray)\n",
        "        vmax = np.nanmax(dataArray)\n",
        "        \n",
        "        #print statistics \n",
        "        print('The average of this data is: ',round(average,3),'\\nThe standard deviation is: ',round(stdev,3),'\\nThe median is: ',round(median,3))\n",
        "        print('The range of latitude in this file is: ',min_lat,' to ',max_lat, 'degrees \\nThe range of longitude in this file is: ',min_lon, ' to ',max_lon,' degrees')\n",
        "        #change 'raw_input' to 'input' if an error is shown about the input\n",
        "        is_map=input('\\nWould you like to create a map of this data? Please enter Y or N \\n')\n",
        "        \n",
        "        #if user would like a map, view it\n",
        "        if is_map == 'Y' or is_map == 'y':\n",
        "          ax = plt.axes(projection=ccrs.LambertConformal())\n",
        "          extent = [min_lon, max_lon, min_lat, max_lat]\n",
        "          ax.coastlines() #outlines the continents          \n",
        "          ax.set_extent(extent)\n",
        "          \n",
        "          my_cmap = plt.cm.get_cmap('jet')\n",
        "          my_cmap.set_under('w')\n",
        "          plt.pcolormesh(lon, lat, data, vmin=0.0, vmax=vmax, cmap=my_cmap, transform=ccrs.PlateCarree())\n",
        "          plt.autoscale()\n",
        "          #title the plot\n",
        "          plt.title('{0}\\n {1}'.format(FILE_NAME, map_title))\n",
        "          fig = plt.gcf()\n",
        "          cb = plt.colorbar(shrink = 0.7)\n",
        "          cb.set_label(map_label, fontsize =9, wrap=True)\n",
        "          \"\"\"\n",
        "            grd = m.gridlines(crs=ccrs.PlateCarree(), draw_labels=True, linewidth=2, color='gray', alpha=0.5, linestyle='--')\n",
        "            grd.xlabels_top = None\n",
        "            grd.ylabels_right = None\n",
        "            grd.xformatter = LONGITUDE_FORMATTER\n",
        "            grd.yformatter = LATITUDE_FORMATTER\n",
        "            \"\"\"\n",
        "          # Show the plot window.\n",
        "          plt.show()\n",
        "          #once you close the map it asks if you'd like to save it\n",
        "          #change 'raw_input' to 'input' if an error is shown about the input    \n",
        "          is_save=str(input('\\nWould you like to save this map? Please enter Y or N \\n'))\n",
        "          if is_save == 'Y' or is_save == 'y':\n",
        "            #saves as a png if the user would like\n",
        "            pngfile = '{0}.png'.format(FILE_NAME[:-3])\n",
        "            fig.savefig('/content/drive/My Drive/Colab Notebooks/' + pngfile, dpi = 300, bbox_inches='tight')\n",
        "        #close the hdf5 file \n",
        "        file.close()"
      ],
      "execution_count": null,
      "outputs": []
    }
  ]
}