map factor exploration notebook by @kgoebber

map_factor_calculations_with_MetPy.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e09f71e2",
   "metadata": {},
   "outputs": [],
   "source": [
    "from datetime import datetime\n",
    "\n",
    "import cartopy.crs as ccrs\n",
    "import cartopy.feature as cfeature\n",
    "import matplotlib.pyplot as plt\n",
    "import metpy.calc as mpcalc\n",
    "import metpy.constants as mpconst\n",
    "from metpy.units import units\n",
    "import numpy as np\n",
    "import pint\n",
    "from pyproj import Proj\n",
    "import xarray as xr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fa23793f",
   "metadata": {},
   "outputs": [],
   "source": [
    "def read_gdlist(name):\n",
    "    file = open(name, 'rb')\n",
    "\n",
    "    for _ in range(6):\n",
    "        file.readline()\n",
    "\n",
    "    grid_row = file.readline()\n",
    "    nlon = int(grid_row[-11:-5].strip())\n",
    "    nlat = int(grid_row[-4:].strip())\n",
    "    #print(nlon,nlat)\n",
    "\n",
    "    gempak_array = np.empty((nlat,nlon))\n",
    "\n",
    "    for _ in range(2):\n",
    "        file.readline()\n",
    "\n",
    "    scale_row = file.readline()\n",
    "    scale_factor = int(scale_row[19:].strip())\n",
    "    #print(scale_factor)\n",
    "\n",
    "    for _ in range(92):\n",
    "        file.readline()\n",
    "\n",
    "    col = 0\n",
    "    row = 0\n",
    "    for line in file.readlines():\n",
    "        if (line[:13].strip() != b'') & (col >= nlon):\n",
    "            col = 0\n",
    "            row+=1\n",
    "        for count in [0,1,2,3,4,5,6,7]:\n",
    "            gempak_array[row,col] = line[8+count*9:17+count*9]\n",
    "            col+=1\n",
    "\n",
    "    file.close()\n",
    "    \n",
    "    return (gempak_array, -scale_factor)\n",
    "\n",
    "def error_stats(metpy_data, gempak_data):\n",
    "    ignore = (gempak_data != -9999.0) * ~np.isnan(metpy_data)\n",
    "    \n",
    "    if isinstance(metpy_data.data, pint.Quantity):\n",
    "        print(f'Is Xarray with Quantity? True')\n",
    "        mdata = metpy_data.values[ignore]\n",
    "    elif isinstance(metpy_data, pint.Quantity):\n",
    "        print(f'Is Xarray with Quantity? False')\n",
    "        print(f'Is Metpy Unit Array? True')\n",
    "        mdata = metpy_data.m[ignore]\n",
    "    else:\n",
    "        print(f'Is Xarray with Quantity? False')\n",
    "        print(f'Is Metpy Unit Array? False')\n",
    "        print(f'I got back just a numpy array')\n",
    "        try:\n",
    "            mdata = metpy_data.values[ignore]\n",
    "        except:\n",
    "            mdata = metpy_data[ignore]\n",
    "    gempak_data = gempak_data[ignore]\n",
    "    \n",
    "    print()\n",
    "    print('Mean Comparison')\n",
    "    print(f'  Mean Values (MetPy): {mdata.mean()}')\n",
    "    print(f'  Mean Values (GEMPAK): {gempak_data.mean()}')\n",
    "    print()\n",
    "    print('Max Comparison')\n",
    "    print(f'  Max Values (MetPy): {mdata.max()}')\n",
    "    print(f'  Max Values (GEMPAK): {gempak_data.max()}')\n",
    "    print()\n",
    "    print('Min Comparison')\n",
    "    print(f'  Min Values (MetPy): {mdata.min()}')\n",
    "    print(f'  Min Values (GEMPAK): {gempak_data.min()}')\n",
    "\n",
    "    print()\n",
    "    print('Difference Array')\n",
    "    diff = mdata - gempak_data\n",
    "    print(diff)\n",
    "    print()\n",
    "    print('Various Statistical Analyses')\n",
    "    print(f'  Average Absolute Difference: {np.nanmean(np.abs(diff))}')\n",
    "    print(f'  RMS Error: {np.sqrt(np.nansum(diff**2))/len(diff.ravel())}')\n",
    "    print(f'  Standard Deviation of Difference: {np.nanstd(diff)}')\n",
    "    print(f'  Max Diff: {np.max(diff)}')\n",
    "    print(f'  Min Diff: {np.min(diff)}')\n",
    "    print(f'  Correlation: {np.corrcoef(mdata.ravel(), gempak_data.ravel())[0][1]}')\n",
    "    print(f'  Relative Magnitude Difference: {np.nanmean(np.abs(diff))/np.nanmax(mdata)}')\n",
    "    print()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "872eea3d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# To get data sources uncomment and run these lines\n",
    "#!wget http://bergeron.valpo.edu/testing_data/gempak_absolute_vorticity.out\n",
    "#!wget http://bergeron.valpo.edu/testing_data/gempak_vorticity.out\n",
    "#!wget http://bergeron.valpo.edu/testing_data/gfs_test_data.nc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3b6ed0d5",
   "metadata": {},
   "outputs": [],
   "source": [
    "gfs_data = xr.open_dataset('gfs_test_data.nc').metpy.parse_cf()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "99eb1c1e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Vorticity\n",
    "\n",
    "gempak_out_var, scale = read_gdlist('gempak_vorticity.out')\n",
    "\n",
    "out_var = gempak_out_var * 10**scale\n",
    "\n",
    "data_var_u = gfs_data['u-component_of_wind_isobaric'].sel(time=datetime(2018, 3, 8, 0), isobaric2=50000).metpy.quantify()\n",
    "data_var_v = gfs_data['v-component_of_wind_isobaric'].sel(time=datetime(2018, 3, 8, 0), isobaric2=50000).metpy.quantify()\n",
    "\n",
    "data_var_lat = gfs_data.lat.values\n",
    "\n",
    "vor = mpcalc.vorticity(data_var_u, data_var_v)\n",
    "\n",
    "print('Vorticity \\n')\n",
    "error_stats(vor[3:-3,3:-3]*1e5, out_var[3:-3,3:-3]*1e5)\n",
    "\n",
    "vor_corrected = vor + data_var_u / mpconst.earth_avg_radius * np.tan(np.deg2rad(data_var_lat[:, None]))\n",
    "\n",
    "print('\\nVorticity with spherical corrections \\n')\n",
    "error_stats(vor_corrected[3:-3,3:-3]*1e5, out_var[3:-3,3:-3]*1e5)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e696d61f",
   "metadata": {},
   "source": [
    "## Begin Testing Calcualtions with Map Factor"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a22afdaa",
   "metadata": {},
   "outputs": [],
   "source": [
    "xx, yy = np.meshgrid(gfs_data.lon.values, gfs_data.lat.values)\n",
    "factors = Proj(data_var_u.metpy.pyproj_crs).get_factors(xx, yy)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "061e99bc",
   "metadata": {},
   "outputs": [],
   "source": [
    "m_x = factors.meridional_scale\n",
    "m_y = factors.parallel_scale"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5eeeeb28",
   "metadata": {},
   "outputs": [],
   "source": [
    "dx, dy = mpcalc.lat_lon_grid_deltas(gfs_data.lon, gfs_data.lat)\n",
    "dudy = mpcalc.first_derivative(data_var_u, delta=dy, axis=-2)\n",
    "dvdx = mpcalc.first_derivative(data_var_v, delta=dx, axis=-1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e21c8aed",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculating Vorticity using Map Factors\n",
    "zeta = (m_x * dvdx - m_y * dudy\n",
    "        - data_var_v * (m_y/m_x) * mpcalc.first_derivative(m_x, delta=dx, axis=-1)\n",
    "        + data_var_u * (m_x/m_y) * mpcalc.first_derivative(m_y, delta=dy, axis=-2))\n",
    "avor = zeta + mpcalc.coriolis_parameter(gfs_data.lat * units.degree)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "16c1e80d",
   "metadata": {},
   "outputs": [],
   "source": [
    "error_stats(zeta[3:-3,3:-3]*1e5, out_var[3:-3,3:-3]*1e5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3b23baef",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.8.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}