modis_prepro_final.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "from pyhdf.SD import SD, SDC, SDS\n",
    "from pyproj import CRS, Proj, Transformer\n",
    "import geopandas as gpd\n",
    "from shapely.geometry import Point\n",
    "import datetime"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def modis_preprocessor(filelist, location_data, storage_path, modis_data_path):\n",
    "    \"\"\"\n",
    "    filelist: list of modis files to be preprocessed. File list must be sorted by date of observation\n",
    "    location_data: numpy array containing latitude and longitude of location of interest\n",
    "\n",
    "    returns: \n",
    "        furnished: pandas dataframe with following columns - ['date', 'latitude', 'longitude', 'AOD47', 'AOD47_3', 'valid3', 'AOD47_5', 'valid5']\n",
    "    \n",
    "    \"\"\"\n",
    "\n",
    "    ls = []\n",
    "    count = 0\n",
    "\n",
    "    for i in filelist:\n",
    "        filename = str(i)\n",
    "\n",
    "        date = datetime.datetime.strptime(filename[11:16], '%y%j').date().strftime(\"%d-%m-%Y %H:%M\")\n",
    "        hdf_file = SD(os.path.join(modis_data_path, filename), SDC.READ)\n",
    "\n",
    "        #collecting required datasets\n",
    "        AOD_47 = hdf_file.select('Optical_Depth_047')\n",
    "        AOD_55 = hdf_file.select('Optical_Depth_055')\n",
    "        AOD_QA = hdf_file.select('AOD_QA')\n",
    "        CWV = hdf_file.select('Column_WV')\n",
    "\n",
    "        #seperation and parsing of meta data for grid-1 with 1km spatial resolution\n",
    "        struct_meta_str = hdf_file.attributes()['StructMetadata.0']\n",
    "        grid1 = struct_meta_str.split(\"END_GROUP=GRID_1\")[0]\n",
    "        grid1_meta = dict([x.split(\"=\") for x in grid1.split() if \"=\" in x])\n",
    "\n",
    "        ###############################################################\n",
    "\n",
    "        \"\"\" Creation of meshgrid\"\"\"\n",
    "\n",
    "        x0, y0 = eval(grid1_meta[\"UpperLeftPointMtrs\"])\n",
    "        x1, y1 = eval(grid1_meta[\"LowerRightMtrs\"])\n",
    "\n",
    "        # Interpolate points between corners, inclusive of bounds\n",
    "        shape = AOD_47.info()[2]\n",
    "        x = np.linspace(x0, x1, shape[2], endpoint=True)\n",
    "        y = np.linspace(y0, y1, shape[1], endpoint=True)\n",
    "\n",
    "        # Return two 2D arrays representing X & Y coordinates of all points\n",
    "        xv, yv = np.meshgrid(x, y)\n",
    "\n",
    "\n",
    "        ###############################################################\n",
    "\n",
    "        \"\"\" Reprojecting data into ESPG:4326 format \"\"\"\n",
    "\n",
    "        crs_info = {\n",
    "            \"crs\": grid1_meta[\"Projection\"],\n",
    "            \"crs_params\": eval(grid1_meta[\"ProjParams\"])\n",
    "        }\n",
    "\n",
    "        projfrom = \"+proj={} +lon_0=0 +x_0=0 +y_0=0 +a={} +units=m +no_defs\".format(\n",
    "            'sinu',crs_info['crs_params'][0])\n",
    "        projto = CRS.from_epsg(\"4326\")\n",
    "\n",
    "        transformer = Transformer.from_crs(\n",
    "                projfrom,\n",
    "                projto,\n",
    "                always_xy=True,\n",
    "            )\n",
    "        lon, lat = transformer.transform(xv, yv)\n",
    "\n",
    "        ###################################################################\n",
    "\n",
    "        \"\"\" Data calibration and collection\"\"\"\n",
    "\n",
    "        def data_calib(dataset):\n",
    "\n",
    "            shape = dataset.info()[2]\n",
    "            calib_info = dataset.attributes()\n",
    "            data = dataset.get()\n",
    "\n",
    "            invalid = (\n",
    "                    (data < calib_info[\"valid_range\"][0]) |\n",
    "                    (data > calib_info[\"valid_range\"][1]) |\n",
    "                    (data == calib_info[\"_FillValue\"])\n",
    "                )\n",
    "\n",
    "            masked_data = np.ma.masked_where(invalid, data)\n",
    "            \n",
    "\n",
    "            try:\n",
    "                masked_data = (\n",
    "                    (masked_data - calib_info[\"add_offset\"]) *\n",
    "                    calib_info[\"scale_factor\"]\n",
    "                )\n",
    "            except:\n",
    "                pass\n",
    "\n",
    "            if not calib_info['long_name'] == 'AOD_QA':\n",
    "                masked_data.fill_value = np.nan\n",
    "                masked_data = np.mean(masked_data, axis= 0)\n",
    "                #masked_data = np.max(masked_data, axis= 0)\n",
    "                \n",
    "                \n",
    "            return masked_data\n",
    "\n",
    "        ###################################################################\n",
    "\n",
    "        \"\"\" Searching for data near grid closest to given location \"\"\"\n",
    "\n",
    "        lat_min = lat.min()\n",
    "        lat_max = lat.max()\n",
    "        lon_min = lon.min()\n",
    "        lon_max = lon.max()\n",
    "\n",
    "\n",
    "        def search_loc(inp_lat, inp_lon):\n",
    "            R=6371007\n",
    "            lat1=np.radians(inp_lat)\n",
    "            lat2=np.radians(lat)\n",
    "            delta_lat=np.radians(lat-inp_lat)\n",
    "            delta_lon=np.radians(lon-inp_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",
    "            x,y=np.unravel_index(d.argmin(),d.shape)\n",
    "            return lat[x,y], lon[x,y], x,y\n",
    "\n",
    "\n",
    "        def grid_mean(calib_data,x,y,n):\n",
    "\n",
    "            \"\"\"calib data is preprocessed data\n",
    "                n is the grid size:\n",
    "                Example: n = 3 will given mean for grid of size 3x3\n",
    "                n should be odd\n",
    "\n",
    "            \"\"\"\n",
    "            num = int((n-1)/2)\n",
    "\n",
    "            while(x < num):\n",
    "                x+=1\n",
    "            while(x > data.shape[0]-num - 1) :\n",
    "                x+=1\n",
    "            while(y < num):\n",
    "                y+=1\n",
    "            while(y > data.shape[1]- num - 1):\n",
    "                y-=1\n",
    "            grid = calib_data[x-num:x+num + 1,y-num:y+num +1]\n",
    "            if grid.count() == 0:\n",
    "                value = 0\n",
    "                count = 0\n",
    "                var = 0\n",
    "            else:\n",
    "                value = grid.mean()\n",
    "                count = grid.count()\n",
    "                var = grid.var()\n",
    "            return value , count, var\n",
    "\n",
    "        data = data_calib(AOD_47)\n",
    "        cwv_data = data_calib(CWV)\n",
    "        data55 = data_calib(AOD_55)\n",
    "        for i in range(len(location_data)):\n",
    "            inp_lat = location_data[i][0]\n",
    "            inp_lon = location_data[i][1]\n",
    "            a,b,x,y = search_loc(inp_lat, inp_lon)\n",
    "            var = [date,\n",
    "                inp_lat,\n",
    "                inp_lon, \n",
    "                data[x,y],\n",
    "                grid_mean(data, x, y, 3)[0], \n",
    "                grid_mean(data, x, y, 3)[1], \n",
    "                grid_mean(data, x, y, 5)[0], \n",
    "                grid_mean(data, x, y, 5)[1],\n",
    "                cwv_data[x,y],\n",
    "                grid_mean(data, x, y, 3)[2],\n",
    "                grid_mean(data, x, y, 5)[2],\n",
    "                data55[x,y] ]\n",
    "            #print(data[x,y])\n",
    "            ls.append(var)\n",
    "        \n",
    "        count+=1\n",
    "        # if(count==2):\n",
    "        #     break\n",
    "\n",
    "    furnished = pd.DataFrame(ls, columns= ['date', 'latitude', 'longitude', 'AOD47', 'AOD47_3', 'valid3', 'AOD47_5', 'valid5', 'ColumnWV', 'AOD47_3var', 'AOD47_5var', 'AOD55'])\n",
    "    furnished.to_csv(storage_path)\n",
    "\n",
    "    return furnished\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "modis_data_path = 'Dataset/Modis-2020'\n",
    "filelist = os.listdir(modis_data_path)\n",
    "location_data = pd.read_csv('Dataset/CPCB/loc-lan-updated.csv')\n",
    "location_data = location_data.loc[:39, ['latitude', 'longitude']].astype(float).to_numpy()\n",
    "storage_location = 'Dataset/AOD-csv/aod-jan-Dec-2020.csv'\n",
    "doi = modis_preprocessor(filelist,location_data, storage_location, modis_data_path )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "modis_data_path = 'Dataset/Modis-2019'\n",
    "filelist = os.listdir(modis_data_path)\n",
    "location_data = pd.read_csv('Dataset/CPCB/loc-lan-updated.csv')\n",
    "location_data = location_data.loc[:39, ['latitude', 'longitude']].astype(float).to_numpy()\n",
    "storage_location = 'Dataset/AOD-csv/aod-jan-Dec-2019.csv'\n",
    "doi = modis_preprocessor(filelist,location_data, storage_location, modis_data_path )"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.10.4 ('dublin')",
   "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.10.4"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "475f2a5d657eb88784a6ef3f75fc3749514a0712ebeb5d3de7b4de4c6f3a72a9"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}