read_aod_and_calculate_pm25.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "read_aod_and_calculate_pm25.ipynb",
      "provenance": [],
      "collapsed_sections": []
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "bS43ungmHJdm"
      },
      "source": [
        "**Module**: read_aod_and_calculate_pm25.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 & Pawan Gupta, 2015\n",
        "\n",
        "**Modified for VIIRS**: Aavash Thapa, May 10 2019\n",
        "\n",
        "**Modified for Cartopy by**: Amanda Markert, University of Alabama in Huntsville, June 2019\n",
        "\n",
        "**Purpose**: To extract AOD data from a VIIRS NC file (or series of files), \n",
        "calculate PM 2.5 from the data, and create a map of the results. (***NOTE: the default slope used in this code is meant for demonstrative purposes only and is not meant to be used for quantitative analyses***)\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "5HsmJnE8QS6t"
      },
      "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": "OeZvUf35U1yF"
      },
      "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",
        "#from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER\n",
        "from matplotlib.colors import LinearSegmentedColormap\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": "wze-wuO8Ur-S"
      },
      "source": [
        "#!/usr/bin/python\n",
        "# =============================================================================\n",
        "# Inputs\n",
        "#This uses the file \"fileList.txt\", containing the list of files, in order to read the file\n",
        "\n",
        "try:\n",
        "    from google.colab import drive\n",
        "    drive.mount('/content/drive', force_remount=True)\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",
        "           #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",
        "        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",
        "        #get valid range for AOD SDS\n",
        "        range=data.getncattr(\"valid_range\")\n",
        "        min_range=min(range)\n",
        "        max_range=max(range)\n",
        "        \n",
        "        #get data within valid range\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",
        "        valid_data = data.ravel()\n",
        "        valid_data=[x for x in valid_data if x>=min_range]\n",
        "        valid_data=[x for x in valid_data if x<=max_range]\n",
        "        valid_data=np.asarray(valid_data)\n",
        "        #find the average\n",
        "        average=sum(valid_data)/len(valid_data)\n",
        "        #find the standard deviation\n",
        "        stdev=np.std(valid_data)\n",
        "        #print information\n",
        "        print('\\nThe valid range of values is: ',round(min_range,3), ' to ',round(max_range,3),'\\nThe average is: ',round(average,3),'\\nThe standard deviation is: ',round(stdev,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",
        "        \n",
        "\n",
        "        #asks user if they want to set PM2.5 calculation parameters\n",
        "        user_input=input('\\nWould you like to enter a slope and intercept for PM 2.5 calculation? (Y/N)')\n",
        "        if user_input == 'Y' or user_input == 'y': \n",
        "            slope=input('Please enter a slope: ')\n",
        "            intercept=input('Please enter an intercept: ')\n",
        "        else:\n",
        "            #if not, choose the following:\n",
        "            slope=29.4\n",
        "            intercept=8.8\n",
        "        valid_data=data\n",
        "        pm25=float(slope)*valid_data+float(intercept)\n",
        "        \n",
        "        \n",
        "        \n",
        "        #Asks user if they would like to see a map\n",
        "        is_map=input('\\nWould you like to create a map of this data? Please enter Y or N \\n')\n",
        "        #if user would like a map, view it\n",
        "        if is_map == 'Y' or is_map == 'y':\n",
        "            \n",
        "            #turn fillvalues to NaN\n",
        "            data=pm25.astype(float)\n",
        "            data[np.logical_and(data>=0,data <= 12)]=0\n",
        "            data[np.logical_and(data>12,data <= 35.4)]=1\n",
        "            data[np.logical_and(data>35.4,data <= 55.4)]=2\n",
        "            data[np.logical_and(data>55.4,data <= 150.4)]=3\n",
        "            data[np.logical_and(data>150.4,data <= 250.4)]=4\n",
        "            data[data>250.4]=5\n",
        "            data[data < 0] = np.nan\n",
        "            \n",
        "            #create the map\n",
        "            data = np.ma.masked_array(data, np.isnan(data))\n",
        "            \n",
        "            extent = (min_lon, max_lon, min_lat, max_lat)\n",
        "            m = plt.axes(projection=ccrs.PlateCarree())\n",
        "            m.set_extent(extent)\n",
        "            \n",
        "            my_cmap=LinearSegmentedColormap.from_list('mycmap', ['green','yellow','orange','red','purple','brown'],6)\n",
        "            plt.pcolormesh(lon, lat, data,cmap=my_cmap, transform=ccrs.PlateCarree())\n",
        "            plt.clim(0,6)\n",
        "             \n",
        "            #create colorbar\n",
        "            cb = plt.colorbar()\n",
        "            cb.set_label('AQI Category')\n",
        "            cb.set_ticks([.5, 1.5,2.5,3.5,4.5,5.5])  # force there to be only 7 ticks\n",
        "            cb.set_ticklabels(['Good', 'Moderate', 'Unhealthy for \\nSensitive Groups','Unhealthy','Very Unhealthy','Hazardous'])  # put text labels on them\n",
        "            \n",
        "            m.coastlines()\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",
        "            plt.autoscale()\n",
        "            \n",
        "            #title the plot\n",
        "            plotTitle=FILE_NAME[:-4]\n",
        "            plt.title('{0}\\n {1}'.format(plotTitle, 'PM 2.5'))\n",
        "            fig = plt.gcf()\n",
        "            # Show the plot window.\n",
        "            plt.show()\n",
        "            \n",
        "            #once you close the map it asks if you'd like to save it\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 like5\n",
        "                pngfile = '{0}.png'.format(plotTitle, 'PM 2.5')\n",
        "                fig.savefig('/content/drive/My Drive/Colab Notebooks/' + pngfile)\n",
        "            \n",
        "            \n"
      ],
      "execution_count": null,
      "outputs": []
    }
  ]
}