examples/mkGridsExample04a.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9a3eb609",
   "metadata": {},
   "outputs": [],
   "source": [
    "# conda: gridTools"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5e36ee05",
   "metadata": {},
   "outputs": [],
   "source": [
    "# This is the same as Example#4 except that it uses dask for\n",
    "# array manipulation to keep the memory footprint down.\n",
    "\n",
    "# NOTE: THIS EXAMPLE IS NOT COMPLETE!\n",
    "\n",
    "# In this example, we recreate the polar sterographic IBCAO grid.\n",
    "# Part A: outside the gridtools library\n",
    "# Part B: using the gridtools library\n",
    "# IBCAO grid details\n",
    "# Cartesian grid\n",
    "#   dx: 2500 meters\n",
    "#   dy: 2500 meters\n",
    "# Grid extent:\n",
    "#   X: -2902500 to 2902500\n",
    "#   Y: -2902500 to 2902500\n",
    "#\n",
    "# Projection\n",
    "#   Stereographic\n",
    "#   lat_0=90.0 (latitude of projection center)\n",
    "#   lat_ts=75.0 (latitude where scale is not distorted)\n",
    "#   lon_0 (not specified, at lat_0=90 it can be any value, proj defaults to 0.0)\n",
    "\n",
    "# For proper orientation for plotting to match the IBCAO technical manual\n",
    "# please use lon_0=0.0 for a plotting parameter.\n",
    "\n",
    "# Example 5 is a python script that is similar to this notebook.\n",
    "# Example 5a is the same example but using a slightly different earth radius."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "25873e5f",
   "metadata": {},
   "source": [
    "# Part A: Compute the grid ourselves"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "46603a08",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import xarray as xr\n",
    "import dask.array as da"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4809995c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# IBCAO\n",
    "# Working in cartesian coordinates, all values are in meters\n",
    "# NOTE: To create a true MOM6 supergrid, the cell spacing\n",
    "# is half the length of a full grid cell.\n",
    "dx = 2500. / 2.\n",
    "dy = 2500. / 2.\n",
    "x = da.arange(-2902500., 2902500. + dx, dx, dtype=np.float32, chunks = (100))\n",
    "y = da.arange(-2902500., 2902500. + dy, dy, dtype=np.float32, chunks = (100))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "45cba2b0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# This is the grid size\n",
    "x.shape, y.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dfeb2609",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create cartesian grid (y, x)\n",
    "# Confirm grid points and x, y coordinates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e8814ac7",
   "metadata": {},
   "outputs": [],
   "source": [
    "xx, yy = da.meshgrid(x, y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4f56b3dd",
   "metadata": {},
   "outputs": [],
   "source": [
    "# From the command line we can convert single points\n",
    "# gdaltransform -s_srs \"+ellps=WGS84 +proj=stere +lat_0=90 +lat_ts=75\" -t_srs EPSG:4326\n",
    "# The pyproj library provides for 2D array transformation of coordinates.\n",
    "\n",
    "from pyproj import CRS, Transformer\n",
    "\n",
    "PROJSTRING = \"+proj=stere +lat_0=90 +lat_ts=75 +ellps=WGS84\"\n",
    "\n",
    "# create the coordinate reference system\n",
    "crs = CRS.from_proj4(PROJSTRING)\n",
    "# create the projection from lon/lat to x/y\n",
    "proj = Transformer.from_crs(crs.geodetic_crs, crs)\n",
    "\n",
    "# compute the lon/lat\n",
    "# ORIGINAL\n",
    "#lon, lat = proj.transform(yy, xx, direction='INVERSE')\n",
    "\n",
    "# DASK\n",
    "def gufoo(xx, yy):\n",
    "    return proj.transform(xx, yy, direction='INVERSE')\n",
    "lon, lat = da.apply_gufunc(gufoo, \"(j),(i)->(j),(i)\", xx, yy, vectorize=True, allow_rechunk=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8bd04005",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Confirm we have the correct grid points and lat lon values\n",
    "# Expected answers:\n",
    "# -2902500.0 -2902500.0 53.8170746379705 -45.0\n",
    "# 2902500.0 2902500.0 53.8170746379705 135.0\n",
    "# REF: https://www.ngdc.noaa.gov/mgg/bathymetry/arctic/IBCAO_TechnicalReference.PDF\n",
    "print(float(yy[0,0]), float(xx[0,0]), float(lat[0,0]), float(lon[0,0]))\n",
    "print(float(yy[y.shape[0]-1, x.shape[0]-1]), float(xx[y.shape[0]-1, x.shape[0]-1]), float(lat[y.shape[0]-1, x.shape[0]-1]), float(lon[y.shape[0]-1, x.shape[0]-1]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e132d4ca",
   "metadata": {},
   "outputs": [],
   "source": [
    "from gridtools.gridutils import GridUtils\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "62b8e59b",
   "metadata": {},
   "outputs": [],
   "source": [
    "grd = GridUtils()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f8926a23",
   "metadata": {},
   "outputs": [],
   "source": [
    "grd.clearGrid()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3a047ef9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define IBCAO grid for gridtools library\n",
    "grd.setGridParameters({\n",
    "    'projection': {\n",
    "        'name': \"Stereographic\",\n",
    "        'ellps': 'WGS84',\n",
    "        'lon_0': 0.0,\n",
    "        'lat_0': 90.0,\n",
    "        'lat_ts': 75.0,\n",
    "        'dx': 5000.0,\n",
    "        'dy': 5000.0,\n",
    "        'dxUnits': 'meters',\n",
    "        'dxUnits': 'meters'\n",
    "    }\n",
    "})\n",
    "grd.grid['x'] = (('nyp','nxp'), lon)\n",
    "grd.grid['y'] = (('nyp','nxp'), lat)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4e440953",
   "metadata": {},
   "outputs": [],
   "source": [
    "grd.grid"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f59e57a3",
   "metadata": {},
   "outputs": [],
   "source": [
    "grd.setPlotParameters(\n",
    "    {\n",
    "        'figsize': (8,8),\n",
    "        'projection': {\n",
    "            'name': 'Stereographic',    \n",
    "            'lon_0': 0.0,\n",
    "            'lat_0': 90.0\n",
    "        },\n",
    "        'extent': [-180, 180, 60, 90],\n",
    "        'iLinewidth': 1,\n",
    "        'jLinewidth': 1,\n",
    "        'showGrid': True,\n",
    "        'title': 'Stereographic: IBCAO',\n",
    "        'iColor': 'y',\n",
    "        'jColor': 'k'\n",
    "    }\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f64ba2a9",
   "metadata": {},
   "outputs": [],
   "source": [
    "(figure, axes) = grd.plotGrid()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "790fc9aa",
   "metadata": {},
   "outputs": [],
   "source": [
    "figure"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4467d94f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute grid metrics\n",
    "grd.computeGridMetricsSpherical()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b2dda06c",
   "metadata": {},
   "outputs": [],
   "source": [
    "grd.grid"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3ddfa3b8",
   "metadata": {},
   "source": [
    "# Part B: Use the gridtools library"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "03a5c42d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Clear the grid and plot parameters\n",
    "grd.clearGrid()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bcedc406",
   "metadata": {},
   "outputs": [],
   "source": [
    "# From the technical documentation, we know:\n",
    "#  * the center point,\n",
    "#  * the grid extent\n",
    "#  * the grid resolution\n",
    "#  * the projection: standard ellipse and polar projection\n",
    "#  * the latitude of true scale\n",
    "#\n",
    "# The full size of the grid is 2902500 + 2902500 meters (5805000 meters)\n",
    "# North Pole: lat_0 = 90.0\n",
    "# Latitude of true scale: lat_ts = 75.0\n",
    "grd.setGridParameters({\n",
    "    'projection': {\n",
    "        'name': \"Stereographic\",\n",
    "        'ellps': 'WGS84',\n",
    "        'lon_0': 0.0,\n",
    "        'lat_0': 90.0,\n",
    "        'lat_ts': 75.0,\n",
    "    },\n",
    "    'centerX': 0.0,\n",
    "    'centerY': 90.0,\n",
    "    'cneterUnits': 'degrees',\n",
    "    'dx': 5805000.0,\n",
    "    'dy': 5805000.0,\n",
    "    'dxUnits': 'meters',\n",
    "    'dyUnits': 'meters',\n",
    "    'gridResolution': 2500.0,\n",
    "    'gridResolutionUnits': 'meters',\n",
    "    'tilt': 0.0,\n",
    "    'gridMode': 2,\n",
    "    'gridType': 'MOM6',\n",
    "    'ensureEvenI': True,\n",
    "    'ensureEvenJ': True\n",
    "})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "542449a3",
   "metadata": {},
   "outputs": [],
   "source": [
    "grd.makeGrid()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9bc93692",
   "metadata": {},
   "outputs": [],
   "source": [
    "grd.grid"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "175d123c",
   "metadata": {},
   "outputs": [],
   "source": [
    "grd.setPlotParameters(\n",
    "    {\n",
    "        'figsize': (8,8),\n",
    "        'projection': {\n",
    "            'name': 'Stereographic',    \n",
    "            'lon_0': 0.0,\n",
    "            'lat_0': 90.0\n",
    "        },\n",
    "        'extent': [-180, 180, 60, 90],\n",
    "        'iLinewidth': 1,\n",
    "        'jLinewidth': 1,\n",
    "        'showGrid': True,\n",
    "        'title': 'Stereographic: IBCAO',\n",
    "        'iColor': 'y',\n",
    "        'jColor': 'k'\n",
    "    }\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a7f60fde",
   "metadata": {},
   "outputs": [],
   "source": [
    "(figure, axes) = grd.plotGrid()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "384e1c92",
   "metadata": {},
   "outputs": [],
   "source": [
    "figure"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d7a8a96b-5a4e-493c-b1d8-6296ff72910c",
   "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.7.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}