read_viirs_at_a_location.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "read_viirs_at_a_location.ipynb",
      "provenance": [],
      "collapsed_sections": []
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "6TjFLdFV7o3R"
      },
      "source": [
        "\n",
        "\n",
        "**Module**: read_viirs_at_a_location.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, 2015 \n",
        "\n",
        "**Modified**: Vikalp Mishra & Pawan Gupta, May 10 2019 to read netcdf data\n",
        "\n",
        "**Modified**: Aavash Thapa, May 10 2019 to read VIIRS data\n",
        "\n",
        "**Purpose**: Identify points that are closest to the user-input latitude and longitude, and output the value of the closest collocation, and the means of a 3x3 and 5x5 pixel square centered around the input.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "cXtxOJN-ob8S"
      },
      "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": "gG9b8yQK3s3O"
      },
      "source": [
        "! pip install netCDF4\n",
        "#import necessary modules\n",
        "import numpy as np\n",
        "import sys\n",
        "from numpy import unravel_index\n",
        "from netCDF4 import Dataset"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "QcnG3Lo-3SWy"
      },
      "source": [
        "#!/usr/bin/python\n",
        "\n",
        "#This finds the user's current path so that all hdf4 files can be found\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",
        "    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",
        "        if 'AERDB' in FILE_NAME:\n",
        "            print('This is a VIIRS Deep Blue file.')\n",
        "            #this is how you access the data tree in an hdf5 file\n",
        "            #SDS_NAME='Aerosol_Optical_Thickness_550_Land'\n",
        "            SDS_NAME='Aerosol_Optical_Thickness_550_Land_Best_Estimate'\n",
        "        ds=file\n",
        "        #grp='PRODUCT'        \n",
        "        lat= ds.variables['Latitude'][:][:]\n",
        "        lon= ds.variables['Longitude'][:][:]\n",
        "        data= ds.variables[SDS_NAME] \n",
        "        \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 the data as an array and mask fill/missing values\n",
        "        dataArray=np.array(data[:][:])\n",
        "        dataArray[dataArray==fv]=np.nan\n",
        "        data=dataArray\n",
        "        \n",
        "        #get statistics about data\n",
        "        average=np.nanmean(dataArray)\n",
        "        stdev=np.nanstd(dataArray)\n",
        "        median=np.nanmedian(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",
        "        \n",
        "        user_lat=float(input('\\nPlease enter the latitude you would like to analyze (Deg. N): '))\n",
        "        user_lon=float(input('Please enter the longitude you would like to analyze (Deg. E): '))\n",
        "        #Continues to ask for lat and lon until the user enters valid values\n",
        "        while user_lat < min_lat or user_lat > max_lat:\n",
        "            user_lat=float(input('The latitude you entered is out of range. Please enter a valid latitude: '))\n",
        "        while user_lon < min_lon or user_lon > max_lon:\n",
        "            user_lon=float(input('The longitude you entered is out of range. Please enter a valid longitude: '))\n",
        "            \n",
        "        #calculation to find nearest point in data to entered location (haversine formula)\n",
        "        R=6371000#radius of the earth in meters\n",
        "        lat1=np.radians(user_lat)\n",
        "        lat2=np.radians(lat)\n",
        "        delta_lat=np.radians(lat-user_lat)\n",
        "        delta_lon=np.radians(lon-user_lon)\n",
        "        a=(np.sin(delta_lat/2))*(np.sin(delta_lat/2))+(np.cos(lat1))*(np.cos(lat2))*(np.sin(delta_lon/2))*(np.sin(delta_lon/2))\n",
        "        c=2*np.arctan2(np.sqrt(a),np.sqrt(1-a))\n",
        "        d=R*c\n",
        "        #gets (and then prints) the x,y location of the nearest point in data to entered location, accounting for no data values\n",
        "        x,y=np.unravel_index(d.argmin(),d.shape)\n",
        "        print('\\nThe nearest pixel to your entered location is at: \\nLatitude:',lat[x,y],' Longitude:',lon[x,y])\n",
        "\n",
        "        if np.isnan(dataArray[x,y]):\n",
        "            print('The value of ',SDS_NAME,'at this pixel is', fv ,'\\n')\n",
        "        elif dataArray[x,y] != fv:\n",
        "            print('The value of ', SDS_NAME, 'at this pixel is ', dataArray[x,y])\n",
        "            \n",
        "        #calculates mean, median, stdev in a 3x3 grid around nearest point to entered location\n",
        "        if x < 1:\n",
        "            x+=1\n",
        "        if x > dataArray.shape[0]-2:\n",
        "            x-=2\n",
        "        if y < 1:\n",
        "            y+=1\n",
        "        if y > dataArray.shape[1]-2:\n",
        "            y-=2\n",
        "        three_by_three=dataArray[x-1:x+2,y-1:y+2]\n",
        "        three_by_three=three_by_three.astype(float)\n",
        "        three_by_three[three_by_three==float(fv)]=np.nan\n",
        "        nnan=np.count_nonzero(~np.isnan(three_by_three))\n",
        "        if nnan == 0:\n",
        "            print ('There are no valid pixels in a 3x3 grid centered at your entered location.')\n",
        "        else:\n",
        "            three_by_three_average=np.nanmean(three_by_three)\n",
        "            three_by_three_std=np.nanstd(three_by_three)\n",
        "            three_by_three_median=np.nanmedian(three_by_three)\n",
        "            if nnan == 1:\n",
        "                npixels='is'\n",
        "                mpixels='pixel'\n",
        "            else:\n",
        "                npixels='are'\n",
        "                mpixels='pixels'\n",
        "            print('There',npixels,nnan,'valid',mpixels,'in a 3x3 grid centered at your entered location.')\n",
        "            print('The average value in this grid is: ',round(three_by_three_average,3),' \\nThe median value in this grid is: ',round(three_by_three_median, 3),'\\nThe standard deviation in this grid is: ', round(three_by_three_std, 3))\n",
        "        \n",
        "        #calculates mean, median, stdev in a 5x5 grid around nearest point to entered location\n",
        "        if x < 2:\n",
        "            x+=1\n",
        "        if x > dataArray.shape[0]-3:\n",
        "            x-=1\n",
        "        if y < 2:\n",
        "            y+=1\n",
        "        if y > dataArray.shape[1]-3:\n",
        "            y-=1\n",
        "        five_by_five=dataArray[x-2:x+3,y-2:y+3]\n",
        "        five_by_five=five_by_five.astype(float)\n",
        "        five_by_five[five_by_five==float(fv)]=np.nan\n",
        "        nnan=np.count_nonzero(~np.isnan(five_by_five))\n",
        "        if nnan == 0:\n",
        "            print ('There are no valid pixels in a 5x5 grid centered at your entered location. \\n')\n",
        "        else:\n",
        "            five_by_five_average=np.nanmean(five_by_five)\n",
        "            five_by_five_std=np.nanstd(five_by_five)\n",
        "            five_by_five_median=np.nanmedian(five_by_five)\n",
        "            if nnan == 1:\n",
        "                npixels='is'\n",
        "                mpixels='pixel'\n",
        "            else:\n",
        "                npixels='are'\n",
        "                mpixels='pixels'\n",
        "            print('\\nThere',npixels,nnan,' valid',mpixels,' in a 5x5 grid centered at your entered location. \\n')\n",
        "            print('The average value in this grid is: ',round(five_by_five_average, 3),' \\nThe median value in this grid is: ',round(five_by_five_median, 3),'\\nThe standard deviation in this grid is: ',round(five_by_five_std, 3))"
      ],
      "execution_count": null,
      "outputs": []
    }
  ]
}