DO ecology extraction.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# DO Ecology Extraction\n",
    "\n",
    "This notebook takes Ecology's bounding scenario model output files in NetCDF format (available at https://fortress.wa.gov/ecy/ezshare/EAP/SalishSea/SalishSeaModelBoundingScenarios.html) and extracts the bottom DO values at just the nodes of interest (based on the domain nodes shapefile created using [ProcessGrid](ProcessGrid.ipynb)). A new NetCDF file is created.\n",
    "\n",
    "The Ecology NetCDF files are over 370 GB each. The CDF file output by this notebook is much smaller, about 450 MB depending on the size of the domain, which makes later processing steps easier to repeat."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "exist_cdf = \"/home/benr/wqmodels/ssm/ecology_bounding_outputs/2008_SSM4_WQ_Exist1_nodes.nc\"\n",
    "reference_cdf = \"/home/benr/wqmodels/ssm/ecology_bounding_outputs/2008_SSM4_WQ_Ref1_nodes.nc\"\n",
    "domain_nodes_shp = \"gis/ssm domain nodes.shp\"\n",
    "masked_nodes_txt = \"gis/masked nodes.txt\"\n",
    "\n",
    "do_output_cdf = \"model_results/bottom do 2008.nc\"\n",
    "\n",
    "from netCDF4 import Dataset\n",
    "import geopandas as gpd\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "\n",
    "from fvcom.grid import FvcomGrid\n",
    "from fvcom.control_volume import ControlVolume"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load the shapefile containing the domain nodes as a GeoDataFrame. This gives us the node IDs to extract data for"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "ERROR 1: PROJ: proj_create_from_database: Open of /home/benr/mambaforge/envs/ssm-analysis/share/proj failed\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>depth</th>\n",
       "      <th>geometry</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>node_id</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>4369</th>\n",
       "      <td>45.184</td>\n",
       "      <td>POINT (515771.000 5333564.000)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4370</th>\n",
       "      <td>51.814</td>\n",
       "      <td>POINT (515865.000 5334886.000)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4371</th>\n",
       "      <td>51.814</td>\n",
       "      <td>POINT (516496.000 5336059.000)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4372</th>\n",
       "      <td>55.545</td>\n",
       "      <td>POINT (517099.000 5337226.000)</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4373</th>\n",
       "      <td>60.431</td>\n",
       "      <td>POINT (518039.000 5338339.000)</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "          depth                        geometry\n",
       "node_id                                        \n",
       "4369     45.184  POINT (515771.000 5333564.000)\n",
       "4370     51.814  POINT (515865.000 5334886.000)\n",
       "4371     51.814  POINT (516496.000 5336059.000)\n",
       "4372     55.545  POINT (517099.000 5337226.000)\n",
       "4373     60.431  POINT (518039.000 5338339.000)"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "domain_nodes = gpd.read_file(domain_nodes_shp).set_index('node_id')\n",
    "domain_nodes.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "grid = FvcomGrid.from_mesh('SSM_Grid/ssm_grid.2dm')\n",
    "masked_nodes = np.loadtxt(masked_nodes_txt).astype(int)\n",
    "cv = ControlVolume(grid=grid, nodes=set(domain_nodes.index)) - set(masked_nodes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<class 'netCDF4._netCDF4.Dataset'>\n",
       "root group (NETCDF4 data model, file format HDF5):\n",
       "    dimensions(sizes): IJK(160120), Time(8760)\n",
       "    variables(dimensions): float32 Var_1(Time, IJK), float32 Var_2(Time, IJK), float32 Var_3(Time, IJK), float32 Var_4(Time, IJK), float32 Var_5(Time, IJK), float32 Var_6(Time, IJK), float32 Var_7(Time, IJK), float32 Var_8(Time, IJK), float32 Var_9(Time, IJK), float32 Var_10(Time, IJK), float32 Var_11(Time, IJK), float32 Var_12(Time, IJK), float32 Var_13(Time, IJK), float32 Var_14(Time, IJK), float32 Var_15(Time, IJK), float32 Var_16(Time, IJK), float32 Var_17(Time, IJK), float32 Var_18(Time, IJK), float32 Var_19(Time, IJK), float32 Var_20(Time, IJK), float32 Var_21(Time, IJK), float32 Var_22(Time, IJK), float32 Var_23(Time, IJK), float32 Var_24(Time, IJK), float32 Var_25(Time, IJK), float32 Var_26(Time, IJK), float32 Var_27(Time, IJK), float32 Var_28(Time, IJK), float32 Var_29(Time, IJK), float32 Var_30(Time, IJK), float32 Var_31(Time, IJK), float32 Var_32(Time, IJK), float32 Var_33(Time, IJK), float32 Var_34(Time, IJK), float32 Var_35(Time, IJK), float32 Var_36(Time, IJK), float32 Var_37(Time, IJK), float32 Var_38(Time, IJK), float32 Var_39(Time, IJK), float32 Var_40(Time, IJK), float32 Var_41(Time, IJK), float32 Var_42(Time, IJK), float32 Var_43(Time, IJK), float32 Var_44(Time, IJK), float32 Var_45(Time, IJK), float32 Var_46(Time, IJK), float32 Var_47(Time, IJK), float32 Var_48(Time, IJK), float32 Var_49(Time, IJK), float32 Var_50(Time, IJK), float32 Var_51(Time, IJK), float32 Var_52(Time, IJK)\n",
       "    groups: "
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "exist = Dataset(exist_cdf)\n",
    "exist"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define the values of the dimension IJK we want.\n",
    "\n",
    "IJK is a representation of the 10 depth points per node, zero-indexed, so to get the bottom points we need to multiply the node number minus 1 by 10 and add 9. For instance, if we wanted the bottom point of node 1 we'd get IJK index 9, and for node 2 we'd get IJK index 19. This simplifies to the expression in the cell below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 43689,  43699,  43709, ..., 159939, 159999, 160029])"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3164\n"
     ]
    }
   ],
   "source": [
    "node_ids = np.array(cv.nodes_list)\n",
    "ijk_index = node_ids * 10 - 1\n",
    "display(ijk_index)\n",
    "print(len(ijk_index))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Extract the existing-condition bottom DO"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "masked_array(\n",
       "  data=[[ 5.88764,  6.36058,  6.74573, ..., 11.1078 , 10.9568 , 10.7358 ],\n",
       "        [ 5.92326,  6.37053,  6.84681, ..., 10.6618 , 11.3917 , 11.1215 ],\n",
       "        [ 5.22892,  5.70077,  6.14503, ..., 10.6665 , 11.1288 , 11.2379 ],\n",
       "        ...,\n",
       "        [ 6.49599,  6.451  ,  6.48955, ..., 10.9775 , 11.2026 , 11.215  ],\n",
       "        [ 6.59279,  6.74018,  7.01533, ..., 11.068  , 11.2184 , 11.2079 ],\n",
       "        [ 6.63089,  6.77024,  7.02655, ..., 11.1304 , 11.2046 , 11.2178 ]],\n",
       "  mask=False,\n",
       "  fill_value=1e+20,\n",
       "  dtype=float32)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "exist_bottom_do = exist['Var_10'][:,ijk_index]\n",
    "exist_bottom_do"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(8760, 3164)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "exist_bottom_do.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now repeat this for the reference condition"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "masked_array(\n",
       "  data=[[ 5.89682,  6.3745 ,  6.76485, ..., 11.1365 , 10.9926 , 10.7745 ],\n",
       "        [ 5.93117,  6.3823 ,  6.86365, ..., 10.704  , 11.4083 , 11.1629 ],\n",
       "        [ 5.23647,  5.712  ,  6.15984, ..., 10.7015 , 11.1511 , 11.2795 ],\n",
       "        ...,\n",
       "        [ 6.50992,  6.46475,  6.50368, ..., 11.0227 , 11.2494 , 11.2664 ],\n",
       "        [ 6.60842,  6.75745,  7.03503, ..., 11.1131 , 11.2674 , 11.261  ],\n",
       "        [ 6.64774,  6.78907,  7.04787, ..., 11.1765 , 11.2557 , 11.2727 ]],\n",
       "  mask=False,\n",
       "  fill_value=1e+20,\n",
       "  dtype=float32)"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ref = Dataset(reference_cdf)\n",
    "ref_bottom_do = ref['Var_10'][:,ijk_index]\n",
    "ref_bottom_do"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<class 'netCDF4._netCDF4.Dataset'>\n",
       "root group (NETCDF4 data model, file format HDF5):\n",
       "    dimensions(sizes): time(8760), node(3164)\n",
       "    variables(dimensions): int32 node(node), float32 time(time), float32 existing(time, node), float32 reference(time, node)\n",
       "    groups: "
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "do_output = Dataset(do_output_cdf, \"w\")\n",
    "timeDim = do_output.createDimension('time', exist_bottom_do.shape[0])\n",
    "nodeDim = do_output.createDimension('node', exist_bottom_do.shape[1])\n",
    "nodeVar = do_output.createVariable('node', \"i4\", ('node'))\n",
    "do_output['node'][:] = node_ids\n",
    "# Time values are not given in the Ecology output files, so recreate them based on a 1-hour\n",
    "# (1/4-day) interval\n",
    "timeVar = do_output.createVariable('time', \"f4\", ('time'))\n",
    "do_output['time'][:] = np.arange(0, exist_bottom_do.shape[0]/24, 1/24)\n",
    "existVar = do_output.createVariable('existing', \"f4\", ('time','node'))\n",
    "do_output['existing'][:] = exist_bottom_do\n",
    "refVar = do_output.createVariable('reference', \"f4\", ('time','node'))\n",
    "do_output['reference'][:] = ref_bottom_do\n",
    "do_output"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "do_output.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:ssm-analysis]",
   "language": "python",
   "name": "conda-env-ssm-analysis-py"
  },
  "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.10.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}