08_Gridding_operations.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 8.1 Regular 1D grids"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "from netCDF4 import Dataset\n",
    "import h5py\n",
    "import scipy.interpolate\n",
    "from matplotlib import pyplot as plt\n",
    "from cartopy import crs as ccrs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fires = pd.read_csv(\"data/VIIRSNDE_global2018312.v1.0.txt\")\n",
    "fires.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Number of nx and ny points for the grid. 720 nx, 360 ny creates 1.0 degree grid.\n",
    "grid_size = 1.0\n",
    "coverage = [-180.0 , -90.0 , 180.0 , 90.0]\n",
    "\n",
    "num_points_x = int((coverage[2] - coverage[0])/grid_size)\n",
    "num_points_y = int((coverage[3] - coverage[1])/grid_size)\n",
    "\n",
    "print(num_points_x, num_points_y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Using a real step length will skip the end value\n",
    "print(np.mgrid[0:4:1.0])\n",
    "\n",
    "# Using a complex number of bins will include the end value\n",
    "print(np.mgrid[0:4:5j])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nx = complex(0, num_points_x)\n",
    "ny = complex(0, num_points_y)\n",
    "\n",
    "Xnew, Ynew = np.mgrid[coverage[0]:coverage[2]:nx, coverage[1]:coverage[3]:ny]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fire_count = np.zeros([num_points_x, num_points_y])\n",
    "\n",
    "for i, lon in enumerate(fires.Lon):\n",
    "    lat = fires.Lat[i]\n",
    "    \n",
    "    adjlat = ((lat + 90) / grid_size);\n",
    "    adjlon = ((lon + 180 ) / grid_size);\n",
    "    \n",
    "    latbin = int(adjlat)\n",
    "    lonbin = int(adjlon)\n",
    "    \n",
    "    fire_count[lonbin, latbin] = fire_count[lonbin, latbin]+1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fire_count[fire_count == 0] = np.nan"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, (ax1, ax2) = plt.subplots(nrows=2, figsize=[15, 15], \\\n",
    "           subplot_kw={'projection': ccrs.PlateCarree()})\n",
    "\n",
    "ax1.coastlines()\n",
    "ax1.set_global()\n",
    "ax1.scatter(fires['Lon'], fires['Lat'], s=1)\n",
    "\n",
    "ax2.coastlines()\n",
    "ax2.set_global()\n",
    "\n",
    "ax2.pcolormesh(Xnew, Ynew, fire_count, vmin=0, vmax=40)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 8.2 Regular 2D grids\n",
    "\n",
    "**Errata:** The book originally imports the operational file, using the path ```data/3B-HHR.MS.MRG.3IMERG.20160811-S233000-E235959.1410.V06B.HDF5```. This path should be changed to  ```data/3B-HHR.MS.MRG.3IMERG.20160811-S233000-E235959.1410.V06B_thinned.HDF5```; this is a smaller version of the operational file. The correct code is shown below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fname = 'data/3B-HHR.MS.MRG.3IMERG.20160811-S233000-E235959.1410.V06B_thinned.HDF5'\n",
    "imergv6 = h5py.File(fname, 'r')\n",
    "\n",
    "precip = imergv6['Grid/precipitationCal'][0,:,:]\n",
    "lat = imergv6['Grid/lat'][:]\n",
    "lon = imergv6['Grid/lon'][:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Xold, Yold = np.meshgrid(lon, lat, indexing='ij')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Number of nx and ny points for the grid. 720 nx, 360 ny creates 1.0 degree grid.\n",
    "def create_2d_grid(grid_size):\n",
    "    coverage = [-180.0 , -90.0 , 180.0 , 90.0]\n",
    "\n",
    "    num_points_x = int((coverage[2] - coverage[0])/grid_size)\n",
    "    num_points_y = int((coverage[3] - coverage[1])/grid_size)\n",
    "\n",
    "    nx = complex(0, num_points_x)\n",
    "    ny = complex(0, num_points_y)\n",
    "\n",
    "    Xnew, Ynew = np.mgrid[coverage[0]:coverage[2]:nx, coverage[1]:coverage[3]:ny]\n",
    "    return Xnew, Ynew"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Xnew, Ynew = create_2d_grid(0.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "precip.shape, Yold.shape, Ynew.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "values = precip.flatten()\n",
    "\n",
    "dims = (values.shape[0], 2)\n",
    "points =  np.zeros(dims)\n",
    "\n",
    "points[:, 0] = Xold.flatten()\n",
    "points[:, 1] = Yold.flatten()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gridOut = scipy.interpolate.griddata(points, values, (Xnew, Ynew), method='nearest')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(ncols=2, figsize=[15, 15], \\\n",
    "   subplot_kw={'projection': ccrs.PlateCarree()})\n",
    "\n",
    "for ax in axes:\n",
    "    ax.set_extent([-94.3, -88.8, 28.9, 33.1])\n",
    "    ax.coastlines('10m', color='orange')\n",
    "\n",
    "axes[0].set_title(\"Before regridding\")\n",
    "axes[0].pcolormesh(Xold, Yold, precip, vmin=0, vmax=20)\n",
    "\n",
    "axes[1].set_title(\"After regridding\")\n",
    "axes[1].pcolormesh(Xnew, Ynew, gridOut, vmin=0, vmax=20)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 8.3 Irregular 2D grids\n",
    "### 8.3.1 Resizing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fname = \"data/goes-meso/michael/OR_ABI-L1b-RadM1-M3C02_G16_s20182822019282_e20182822019339_c20182822019374.nc\"\n",
    "goesnc = Dataset(fname, 'r')\n",
    "C02 = goesnc.variables['Rad'][:,:]\n",
    "x = goesnc.variables['x'][:]\n",
    "y = goesnc.variables['y'][:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "proj_var = goesnc.variables['goes_imager_projection']\n",
    "\n",
    "sat_height = proj_var.perspective_point_height\n",
    "semi_major = proj_var.semi_major_axis\n",
    "semi_minor = proj_var.semi_minor_axis\n",
    "central_lon = proj_var.longitude_of_projection_origin\n",
    "globe = ccrs.Globe(semimajor_axis=semi_major, semiminor_axis=semi_minor)\n",
    "\n",
    "# Multiply the x, y coordinated by satellite height to get the pixel position\n",
    "X = x * sat_height\n",
    "Y = y * sat_height\n",
    "imgExtent = (X.min(), X.max(), Y.min(), Y.max())\n",
    "\n",
    "crs = ccrs.Geostationary(central_longitude=central_lon, satellite_height=sat_height, globe=globe)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "imgExtent"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "resized_C02 = C02[::20, ::20]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lims = [0, 500]\n",
    "\n",
    "fig, axes = plt.subplots(ncols=2,  figsize=[15, 15], \\\n",
    "     subplot_kw={'projection': crs})\n",
    "\n",
    "for ax in axes:\n",
    "    ax.set_extent([-90.0, -82.0, 22.0, 30.0])\n",
    "    ax.coastlines('10m', color='black')\n",
    "\n",
    "axes[0].set_title(\"Native Resolution\")\n",
    "im1 = axes[0].imshow(C02, cmap=plt.get_cmap(\"rainbow\"), origin='upper', \n",
    "            extent=imgExtent, vmin=lims[0], vmax=lims[1])\n",
    "\n",
    "axes[1].set_title(\"20x Lower Resolution\")\n",
    "im2 = axes[1].imshow(resized_C02, cmap=plt.get_cmap(\"rainbow\"), origin='upper', \n",
    "            extent=imgExtent, vmin=lims[0], vmax=lims[1])\n",
    "\n",
    "# Adding colorbars\n",
    "fig.colorbar(im1, pad=0.05, orientation='horizontal', ax=axes[0])\n",
    "fig.colorbar(im2, pad=0.05, orientation='horizontal', ax=axes[1])\n",
    "\n",
    "plt.show()\n",
    "# plt.savefig('C:\\\\Users\\\\rebekah\\\\Dropbox\\\\Python book\\\\AGU_Esmaili_9781119606888_RPP\\\\c08_rbin_flipped.png')\n",
    "# plt.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 8.3.2 Regridding"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.patches as patches"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fname = \"data/goes-meso/michael/OR_ABI-L1b-RadM1-M3C13_G16_s20182822019282_e20182822019350_c20182822019384.nc\"\n",
    "goesnc_sm = Dataset(fname, 'r')\n",
    "C13 = goesnc_sm.variables['Rad'][:]\n",
    "sm_x = goesnc_sm.variables['x'][:]\n",
    "sm_y = goesnc_sm.variables['y'][:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "vmin=0\n",
    "vmax=50\n",
    "\n",
    "delta = 0.005\n",
    "xbot = -0.025\n",
    "ybot =  0.075\n",
    "\n",
    "lower_left = (xbot, ybot)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(ncols=1, nrows=1, figsize=[10, 10])\n",
    "\n",
    "meso = axes.pcolormesh(sm_x, sm_y, C13, vmin=vmin, vmax=vmax, cmap=plt.get_cmap(\"tab20c\"))\n",
    "fig.colorbar(meso, label='Brightness Temperature (K)')\n",
    "zoom_box = patches.Rectangle(lower_left, delta, delta, \n",
    "             linewidth=4, edgecolor='black', facecolor='none')\n",
    "axes.add_patch(zoom_box)\n",
    "axes.set_aspect('equal')\n",
    "axes.axis('off')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Xold, Yold = np.meshgrid(sm_x, sm_y)\n",
    "Xnew, Ynew = np.meshgrid(x, y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "values = C13.flatten()\n",
    "\n",
    "dims = (values.shape[0], 2)\n",
    "points =  np.zeros(dims)\n",
    "\n",
    "points[:, 0] = Xold.flatten()\n",
    "points[:, 1] = Yold.flatten()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gridOut_nn = scipy.interpolate.griddata(points, values, (Xnew, Ynew), method='nearest')\n",
    "gridOut_lin = scipy.interpolate.griddata(points, values, (Xnew, Ynew), method='linear')\n",
    "gridOut_cube = scipy.interpolate.griddata(points, values, (Xnew, Ynew), method='cubic')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gridOut_rebin = np.repeat(np.repeat(C13, 4, axis=1), 4, axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gridOut = [[gridOut_rebin, gridOut_nn], [gridOut_lin, gridOut_cube]]\n",
    "labels = [[\"Re-binning\", \"Nearest Neighbor\"],[\"Linear\", \"Cubic\"]]\n",
    "plt.rcParams.update({'font.size': 18, 'figure.figsize': [20, 20]})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(ncols=2, nrows=2)\n",
    "\n",
    "cmap=plt.get_cmap(\"tab20c\")\n",
    "\n",
    "for i, axis in enumerate(axes):\n",
    "    for j, ax in enumerate(axis):\n",
    "        ax.set_xlim([xbot, xbot+delta])\n",
    "        ax.set_ylim([ybot, ybot+delta])\n",
    "        ax.set_aspect('equal')\n",
    "        ax.set_title(labels[i][j])\n",
    "        ax.pcolormesh(Xnew, Ynew, gridOut[i][j], vmin=vmin, vmax=vmax, cmap=cmap)\n",
    "        ax.axis('off')\n",
    "\n",
    "plt.show()\n",
    "# plt.savefig('C:\\\\Users\\\\rebekah\\\\Dropbox\\\\Python book\\\\AGU_Esmaili_9781119606888_RPP\\\\c08_nn_compare.png')\n",
    "# plt.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 8.3.3 Resampling"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pyresample import geometry\n",
    "from pyresample.kd_tree import resample_nearest"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fname='data/aod/JRR-AOD_v1r1_npp_s201808091955538_e201808091957180_c201808092049460_thinned.nc'\n",
    "file_id_NPP = Dataset(fname)\n",
    "\n",
    "aod = file_id_NPP.variables['AOD550'][:,:]\n",
    "lat = file_id_NPP.variables['Latitude'][:,:]\n",
    "lon = file_id_NPP.variables['Longitude'][:,:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Input list of swath points\n",
    "oldLonLat = geometry.SwathDefinition(lons=lon, lats=lat)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a new grid at 0.1 degree gridding\n",
    "x = np.arange(lon.min(), lon.max(), 0.1)\n",
    "y = np.arange(lat.min(), lat.max(), 0.1)\n",
    "newLon, newLat = np.meshgrid(x, y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# define the new grid using \n",
    "newLonLat = geometry.GridDefinition(lons=newLon, lats=newLat)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Resample the data\n",
    "newAOD = resample_nearest(oldLonLat, aod, newLonLat, radius_of_influence=5000, fill_value=None)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lims = [0,1]\n",
    "fig, axes = plt.subplots(nrows=2,  figsize=[15, 10], \\\n",
    "     subplot_kw={'projection': ccrs.PlateCarree()})\n",
    "\n",
    "for ax in axes:\n",
    "    ax.coastlines('10m', color='black')\n",
    "\n",
    "axes[0].set_title(\"Before regridding\")\n",
    "axes[0].pcolormesh(lon, lat, aod, vmin=lims[0], vmax=lims[1])\n",
    "\n",
    "axes[1].set_title(\"After regridding\")\n",
    "axes[1].pcolormesh(newLon, newLat, newAOD, vmin=lims[0], vmax=lims[1])\n",
    "\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "earthsciviz",
   "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.9.15"
  },
  "vscode": {
   "interpreter": {
    "hash": "48fc10cb14791ea88ba2696b1818b1bf7d4d99ed63d4104785d39779bf72a716"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}