examples/NewGridMOM6.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "65cbc38c-6129-44d3-9e12-48302365de8f",
   "metadata": {},
   "source": [
    "# Read the Docs Tutorial\n",
    "\n",
    "This notebook is part of a tutorial called \"[Build and Edit a MOM6 Grid in Jupyter](https://mom6gridtools.readthedocs.io/en/latest/tutorials/jupyterMOM6.html)\"."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "48edee8b-ca4d-47d6-82c0-f066cd54a8c3",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load some basic python libraries\n",
    "import os, cartopy\n",
    "# This loads the GridUtils portion of the gridtools library\n",
    "from gridtools.gridutils import GridUtils"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1860a3c6-d564-4db8-8a5b-102431d9fbea",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a GridUtils object\n",
    "grd = GridUtils()\n",
    "\n",
    "# Create a grid generation application object\n",
    "grdGenApp = grd.app()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2013bf9b-2f1c-4632-9993-8f53b66fef9d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Launch the grid generation application within Jupyter\n",
    "display(grdGenApp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0196d2ad-ab34-4bc0-a663-b53f68dbea63",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Obtain roughness and topography grids\n",
    "runThisCodeBlock = False\n",
    "if runThisCodeBlock:\n",
    "    # Detach logger from application\n",
    "    grd.detachLoggingFromApplication()\n",
    "    \n",
    "    # Source of GEBCO 2020 topographic grid\n",
    "    highResTopographyFile = \"/import/AKWATERS/jrcermakiii/bathy/gebco/GEBCO_2020.nc\"\n",
    "\n",
    "    if os.path.isfile(highResTopographyFile):\n",
    "        topoGrids = grd.computeBathymetricRoughness(highResTopographyFile,\n",
    "            depthName='elevation',\n",
    "            maxMb=99, superGrid=False, useClipping=False,\n",
    "            auxVariables=['depth'])    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7f47ba60-e394-41dd-95d1-fcc8f402fd7c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Turn the diagnosed topography grid into an actual depth\n",
    "runThisCodeBlock = False\n",
    "if runThisCodeBlock:    \n",
    "    topoGrids['depth'] = -(topoGrids['depth'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "271f8346-a80e-42ae-9e10-5367aa6b8f18",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Write FMS coupler and mosaic files, roughness and topography\n",
    "runThisCodeBlock = False\n",
    "if runThisCodeBlock:\n",
    "    # Write current model grid files\n",
    "    wrkDir = \"/home/cermak/workdir/configs/zOutput\"\n",
    "    inputDir = os.path.join(wrkDir, \"INPUT\")\n",
    "    input2Dir = os.path.join(wrkDir, \"INPUT2\")\n",
    "\n",
    "    # Write FMS coupler and mosaic files\n",
    "    grd.makeSoloMosaic(\n",
    "        topographyGrid=topoGrids['depth'],\n",
    "        writeLandmask=True,\n",
    "        writeOceanmask=True,\n",
    "        inputDirectory=inputDir,\n",
    "        overwrite=True\n",
    "    )\n",
    "\n",
    "    # Write topographic variable\n",
    "    topoGrids.to_netcdf(os.path.join(inputDir, 'ocean_topog.nc'),\n",
    "            encoding=grd.removeFillValueAttributes(data=topoGrids))\n",
    "\n",
    "    # Write the model grid\n",
    "    grd.saveGrid(filename=os.path.join(inputDir, \"ocean_hgrid.nc\"))    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b98c1000-c4e5-4133-a51f-c25f4ba70180",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Write FMS coupler and mosaic files, roughness and topography\n",
    "runThisCodeBlock = False\n",
    "if runThisCodeBlock:\n",
    "    # Examine the topography grid \n",
    "    grd.setPlotParameters({\n",
    "        'figsize': (8,8),\n",
    "        'projection': {\n",
    "            'name': 'LambertConformalConic',\n",
    "            'lon_0': 230.0,\n",
    "            'lat_1': 25.0,\n",
    "            'lat_2': 55.0\n",
    "        },\n",
    "        'extent': [-160.0 ,-100.0, 20.0, 60.0],\n",
    "        'iLinewidth': 1.0,\n",
    "        'jLinewidth': 1.0,\n",
    "        'showGridCells': True,\n",
    "        'iColor': 'k',\n",
    "        'jColor': 'k',\n",
    "        'transform': cartopy.crs.PlateCarree(),\n",
    "        'satelliteHeight': 35785831.0\n",
    "    })\n",
    "    (figure, axes) = grd.plotGrid(showModelGrid = True,\n",
    "            plotVariables={\n",
    "            'depth': {\n",
    "                'values': topoGrids['depth'],\n",
    "                'title': 'Ocean topography (meters)',\n",
    "                'cbar_kwargs': {\n",
    "                    'orientation': 'horizontal',\n",
    "                }\n",
    "            }\n",
    "        })\n",
    "    display(figure)\n",
    "    \n",
    "    # Examine the ocean mask\n",
    "    oceanMask = grd.openDataset(os.path.join(inputDir, 'ocean_mask.nc'))\n",
    "\n",
    "    # Define our own color map (same used in mask editor)\n",
    "    import matplotlib.pyplot as plt\n",
    "    land_color = (0.6, 1.0, 0.6)\n",
    "    sea_color  = (0.6, 0.6, 1.0)\n",
    "    maskCM = plt.matplotlib.colors.ListedColormap(\n",
    "        [land_color, sea_color], name='land/sea')    \n",
    "    \n",
    "    # MOM6 places lon and lat in x and y\n",
    "    # x and y need to be lon and lat coordinates for the mask editor\n",
    "    oceanMask = oceanMask.rename({\n",
    "        'x': 'lon',\n",
    "        'y': 'lat'\n",
    "    })\n",
    "    oceanMask = oceanMask.set_coords(['lon', 'lat'])\n",
    "    \n",
    "    (figureMask, axesMask) = grd.plotGrid(showModelGrid = True,\n",
    "            plotVariables={\n",
    "            'mask': {\n",
    "                'values': oceanMask['mask'],\n",
    "                'title': 'Ocean mask (1 = ocean)',\n",
    "                'cmap': maskCM,\n",
    "                'cbar_kwargs': {\n",
    "                    'orientation': 'horizontal',\n",
    "                }\n",
    "            }\n",
    "        })\n",
    "    display(figureMask)\n",
    "    \n",
    "    # Zoom in to take a closer look\n",
    "    grd.setPlotParameters({\n",
    "        'extent': [-140.0 ,-120.0, 49.0, 59.0]\n",
    "    })\n",
    "    \n",
    "    (figureMaskZoom, axesMaskZoom) = grd.plotGrid(showModelGrid = True,\n",
    "            plotVariables={\n",
    "            'mask': {\n",
    "                'values': oceanMask['mask'],\n",
    "                'title': 'Ocean mask (1 = ocean): Zoom',\n",
    "                'cmap': maskCM,\n",
    "                'cbar_kwargs': {\n",
    "                    'orientation': 'horizontal',\n",
    "                }\n",
    "            }\n",
    "        })\n",
    "    display(figureMaskZoom)    \n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5cba18b8-a36f-490d-bb64-15b02de64cbd",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Use the editor to make some mask updates\n",
    "runThisCodeBlock = False\n",
    "if runThisCodeBlock:\n",
    "    # Load the mask editor application module from gridtools\n",
    "    from gridtools.app import maskEditor\n",
    "    \n",
    "    # Set a map projection for the mask editor to use\n",
    "    crs = cartopy.crs.Orthographic(-140, 45)\n",
    "\n",
    "    # Create the mask editor object\n",
    "    appObj = maskEditor(crs=crs, ds=oceanMask['mask'])\n",
    "    \n",
    "    # Create the mask editor application object\n",
    "    app = appObj.createMaskEditorApp()\n",
    "\n",
    "    # Launch the application\n",
    "    display(app)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "15a61fda-2e77-4a43-ab73-8c5024a34e56",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Save the new ocean mask\n",
    "runThisCodeBlock = False\n",
    "if runThisCodeBlock:\n",
    "    newMask = oceanMask['mask'].copy()\n",
    "    newMask = newMask.reset_coords(names = ['lat', 'lon'])\n",
    "    grd.saveDataset(os.path.join(inputDir, 'ocean_mask_new.nc'), newMask, \n",
    "                    overwrite=True, mapVariables = {'lon': 'x', 'lat': 'y'}, \n",
    "                    hashVariables = ['mask', 'x', 'y'])    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "37f244dc-5046-4ee0-a72f-5bf2abfe98f4",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Apply new ocean mask to ocean model grid\n",
    "runThisCodeBlock = False\n",
    "if runThisCodeBlock:\n",
    "    topoGrids['depth'] = grd.applyExistingOceanmask(topoGrids, 'depth',\n",
    "        os.path.join(inputDir, 'ocean_mask_new.nc'), 'mask',\n",
    "        MASKING_DEPTH=0.0, MINIMUM_DEPTH=0.0, MAXIMUM_DEPTH=-99999.0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eb020d0e-c555-4fdb-9c4c-993b4dfa62a4",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Rewrite FMS coupler and mosaic files\n",
    "runThisCodeBlock = False\n",
    "if runThisCodeBlock:\n",
    "    grd.makeSoloMosaic(\n",
    "        topographyGrid=topoGrids['depth'],\n",
    "        writeLandmask=True,\n",
    "        writeOceanmask=True,\n",
    "        inputDirectory=input2Dir,\n",
    "        overwrite=True,\n",
    "        MASKING_DEPTH=0.0, MINIMUM_DEPTH=0.0, MAXIMUM_DEPTH=-99999.0\n",
    "    )\n",
    "\n",
    "    # Be sure to save previously diagnosed `h2` grid\n",
    "    topoGrids.to_netcdf(os.path.join(input2Dir, 'ocean_topog.nc'),\n",
    "            encoding=grd.removeFillValueAttributes(data=topoGrids))    \n",
    "    \n",
    "    grd.saveGrid(filename=os.path.join(input2Dir, \"ocean_hgrid.nc\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7a933475-8fef-4075-96ef-fe2a27a8631d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Replot the ocean mask to check\n",
    "runThisCodeBlock = False\n",
    "if runThisCodeBlock:\n",
    "    (figureMaskZoom2, axesMaskZoom2) = grd.plotGrid(showModelGrid = True,\n",
    "            plotVariables={\n",
    "            'mask': {\n",
    "                'values': oceanMask['mask'],\n",
    "                'title': 'Ocean mask (1 = ocean): Zoom',\n",
    "                'cmap': maskCM,\n",
    "                'cbar_kwargs': {\n",
    "                    'orientation': 'horizontal',\n",
    "                }\n",
    "            }\n",
    "        })\n",
    "    display(figureMaskZoom2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "45c4c739-e84a-4832-8c29-b3340fe5cb90",
   "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
}