07_Making_maps.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 7.3 Cylindrical Maps\n",
    "### 7.2.1 Global Plots"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "from cartopy import crs as ccrs\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Options to print figures into notebook/increase size\n",
    "plt.rcParams.update({'font.size': 16, 'figure.figsize': [12, 6]})"
   ]
  },
  {
   "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": [
    "fig = plt.figure(figsize=[20,20])\n",
    "\n",
    "ax = plt.subplot(projection=ccrs.PlateCarree())\n",
    "ax.coastlines()\n",
    "ax.set_global()\n",
    "\n",
    "plt.scatter(fires['Lon'], fires['Lat'])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 7.2.2 Changing Projections"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=[20,20])\n",
    "\n",
    "ax = plt.subplot(projection=ccrs.LambertAzimuthalEqualArea())\n",
    "\n",
    "ax.coastlines()\n",
    "ax.gridlines()\n",
    "ax.set_global()\n",
    "\n",
    "plt.scatter(fires['Lon'], fires['Lat'], transform=ccrs.PlateCarree())\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 7.2.3 Regional Plots"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "extent = [-125, -120, 38, 44]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "to_proj = ccrs.PlateCarree()\n",
    "from_proj = ccrs.PlateCarree()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lonLabels = np.arange(-180, 180, 1.5)\n",
    "latLabels = np.arange(-90, 90, 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER\n",
    "from matplotlib import ticker"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure()\n",
    "ax = plt.subplot(projection=to_proj)\n",
    "\n",
    "ax.coastlines('50m')\n",
    "ax.set_extent(extent)\n",
    "\n",
    "frp = ax.scatter(fires['Lon'], fires['Lat'], c=fires['frp(MW)'],\n",
    "    s=fires['frp(MW)'], transform=from_proj)\n",
    "cbar = plt.colorbar(frp)\n",
    "cbar.set_label(\"Fire Radiative Power (MW)\", rotation='vertical')\n",
    "\n",
    "# 1) Maps the gridlines to the variable gl\n",
    "gl = ax.gridlines(crs=to_proj, draw_labels=True)\n",
    "\n",
    "# 2) Adds two attributes to gl, which are xlocator and ylocator\n",
    "gl.xlocator = ticker.FixedLocator(lonLabels)\n",
    "gl.ylocator = ticker.FixedLocator(latLabels)\n",
    "\n",
    "# 3) Changes labels to show degrees North/South and East/West\n",
    "gl.xformatter = LONGITUDE_FORMATTER\n",
    "gl.yformatter = LATITUDE_FORMATTER\n",
    "\n",
    "# 4) Removed labels from top and right side\n",
    "# comment out if you want to include\n",
    "gl.xlabels_top = False\n",
    "gl.ylabels_right = False\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 7.2.4 Swath Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from netCDF4 import Dataset\n",
    "\n",
    "fname='data/aod/JRR-AOD_v1r1_npp_s201808091955538_e201808091957180_c201808092049460_thinned.nc'\n",
    "file_id = Dataset(fname)\n",
    "file_id.variables.keys()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(file_id.variables['AOD550'].shape, \n",
    "      file_id.variables['Latitude'].shape,\n",
    "      file_id.variables['Longitude'].shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "aod = file_id.variables['AOD550'][:,:]\n",
    "lat = file_id.variables['Latitude'][:,:]\n",
    "lon = file_id.variables['Longitude'][:,:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(file_id.variables['AOD550'].valid_range)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "levs = np.arange(0, 1.8, 0.1)\n",
    "levs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=[12, 12])\n",
    "ax = plt.subplot(projection=ccrs.PlateCarree())\n",
    "\n",
    "ax.coastlines('50m')\n",
    "\n",
    "extent = [-125, -120, 38, 44]\n",
    "ax.set_extent(extent)\n",
    "\n",
    "x1 = ax.contourf(lon, lat, aod, levs, extend='both')\n",
    "\n",
    "fig.colorbar(x1, extend='both',\n",
    "    orientation=\"horizontal\", fraction=0.05)\n",
    "\n",
    "# Adds the active fire scatter plot on top\n",
    "ax.scatter(fires['Lon'], fires['Lat'], color='red', s=50)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 7.2.5 Quality Flag Filtering"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import quality flag\n",
    "quality_flag = file_id.variables['QCAll'][:,:]\n",
    "\n",
    "# Keep all but the \"best\" quality using masked arrays\n",
    "maskHQ = (quality_flag != 0)\n",
    "aodHQ = np.ma.masked_where(maskHQ, aod)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "file_id.variables"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(aod.count()-aodHQ.count())/aod.count()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Top plot\n",
    "fig = plt.figure()\n",
    "\n",
    "upper_axis = plt.subplot(2,1,1, projection=ccrs.PlateCarree())\n",
    "upper_axis.set_title(\"All Quality\")\n",
    "\n",
    "upper_axis.coastlines('50m')\n",
    "\n",
    "upper_fig = upper_axis.contourf(lon, lat, aod, levs, extend='both')\n",
    "fig.colorbar(upper_fig, ax=upper_axis, extend='both')\n",
    "\n",
    "# Bottom plot\n",
    "lower_axis = plt.subplot(2,1,2, projection=ccrs.PlateCarree())\n",
    "lower_axis.set_title(\"High Quality\")\n",
    "\n",
    "lower_axis.coastlines('50m')\n",
    "\n",
    "lower_fig = lower_axis.contourf(lon, lat, aodHQ, levs, extend='both')\n",
    "fig.colorbar(lower_fig, ax=lower_axis, extend='both')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 7.3 Polar Stereographic Maps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Cryo-Sat (ESA satellite/data) source: https://n5eil01u.ecs.nsidc.org/ICEBRIDGE/RDEFT4.001/2020.01.01/\n",
    "import numpy.ma as ma\n",
    "import xarray as xr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fname = \"data/RDEFT4_20200131.nc\"\n",
    "ice = xr.open_dataset(fname)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(ice['sea_ice_thickness'].min())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ice_masked = ice.where(ice['sea_ice_thickness'] != -9999.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "to_proj = ccrs.NorthPolarStereo()\n",
    "from_proj = ccrs.PlateCarree()\n",
    "\n",
    "plt.figure(figsize=[10,10])\n",
    "ax = plt.subplot(projection=to_proj)\n",
    "ax.coastlines('50m')\n",
    "ax.gridlines()\n",
    "\n",
    "ice_plot = ax.pcolormesh(ice['lon'], ice['lat'], ice_masked['sea_ice_thickness'], transform=from_proj)\n",
    "plt.colorbar(ice_plot)\n",
    "\n",
    "ax.set_extent([-180, 180, 60, 90], crs=ccrs.PlateCarree())\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 7.4 Geostationary Maps"
   ]
  },
  {
   "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",
    "file_id = Dataset(fname)\n",
    "list(file_id.variables)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "file_id.variables['x'][0:10], file_id.variables['y'][0:10]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "proj_var = file_id.variables['goes_imager_projection']\n",
    "print(proj_var)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define the satellite height and central longitude for plots\n",
    "# Can vary depending on the geo satellite\n",
    "sat_height = proj_var.perspective_point_height\n",
    "\n",
    "# Get the radiance data\n",
    "rad = file_id.variables['Rad'][:,:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define the globe projection\n",
    "semi_major = proj_var.semi_major_axis\n",
    "semi_minor = proj_var.semi_minor_axis\n",
    "\n",
    "globe = ccrs.Globe(semimajor_axis=semi_major, semiminor_axis=semi_minor)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "central_lon = proj_var.longitude_of_projection_origin"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "crs = ccrs.Geostationary(central_longitude=central_lon, satellite_height=sat_height, globe=globe)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = file_id.variables['x'][:] * sat_height\n",
    "Y = file_id.variables['y'][:] * sat_height\n",
    "\n",
    "imgExtent = (X.min(), X.max(), Y.min(), Y.max())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(10,10))\n",
    "\n",
    "ax = plt.subplot(projection=crs)\n",
    "ax.coastlines('10m', color='orange', linewidth=2)\n",
    "\n",
    "ax.imshow(rad, origin='upper', cmap='gray_r', extent=imgExtent)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "glmfname = 'data/GLM-L2-LCFA_2018_282_20_OR_GLM-L2-LCFA_G16_s20182822000200_e20182822000400_c20182822000427.nc'\n",
    "file_id_glm = Dataset(glmfname)\n",
    "file_id_glm.variables.keys()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "glmLon = file_id_glm.variables['event_lon'][:]\n",
    "glmLat = file_id_glm.variables['event_lat'][:]\n",
    "area = file_id_glm.variables['event_energy'][:]\n",
    "\n",
    "glmDF = pd.DataFrame({'lat': glmLat, 'lon': glmLon, 'area': area })"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "crs = ccrs.Geostationary(central_longitude=central_lon, satellite_height=sat_height, globe=globe)\n",
    "\n",
    "from_proj =  ccrs.PlateCarree()\n",
    "plt.figure(figsize=[10,10])\n",
    "\n",
    "ax = plt.subplot(projection=crs)\n",
    "ax.coastlines('10m', color='orange', linewidth=2)\n",
    "ax.set_extent([-90.0, -82.0, 22.0, 30.0])\n",
    "\n",
    "# Plot ABI\n",
    "ax.imshow(rad, origin='upper',cmap='gray_r', extent=imgExtent)\n",
    "\n",
    "# Add GLM data\n",
    "plt.scatter(glmDF.lon, glmDF.lat, c=glmDF.area, transform=from_proj,\n",
    "    s = 300, marker='x')\n",
    "plt.colorbar(extend='both')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 7.5 Plotting with OpenDAP"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import cartopy.feature as cfeature\n",
    "import xarray as xr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "baseURL = 'http://www.esrl.noaa.gov'\n",
    "catalogURL = '/psd/thredds/dodsC/Datasets/noaa.ersst.v5/sst.mnmean.nc'\n",
    "sstID = xr.open_dataset(baseURL+catalogURL)\n",
    "print(sstID)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sst = sstID.sst\n",
    "sst"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mostRecent = len(sst.time.values)-1\n",
    "recentSST = sst.isel(time=mostRecent)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sstmin = 0\n",
    "sstmax = 30"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=[10,5])\n",
    "ax = plt.subplot(projection=ccrs.Orthographic(-90, 0))\n",
    "recentSST.plot(cmap=plt.get_cmap('plasma'),\n",
    "    vmin=sstmin, vmax=sstmax, transform=ccrs.PlateCarree())\n",
    "ax.coastlines('50m')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "levels = np.arange(sstmin, sstmax, 2)\n",
    "\n",
    "fig = plt.figure(figsize=[10,5])\n",
    "ax = plt.subplot(projection=ccrs.Orthographic(-90, 0))\n",
    "recentSST.plot.contourf(levels=levels, cmap=plt.get_cmap('plasma'), transform=ccrs.PlateCarree())\n",
    "ax.coastlines('50m')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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
}