add section on mapping data between grids

csdms · Aug 16, 2023 · 575cfb3 · 575cfb3
1 parent cedd181
commit 575cfb3
Show file tree

Hide file tree

Showing 2 changed files with 187 additions and 0 deletions.
diff --git a/lessons/landlab/landlab/02-data.a.ipynb b/lessons/landlab/landlab/02-data.a.ipynb
@@ -18,6 +18,7 @@
     "# Table of Contents\n",
     "* [Attaching data to a grid](#Attaching-data-to-a-grid)\n",
     "  * [Load data from an ESRI ASCII file](#Load-data-from-an-ESRI-ASCII-file)\n",
+    "  * [Map data from one grid to another](#Map-data-from-one-grid-to-another)\n",
     "  * [Load data from OpenTopography](#Load-data-from-OpenTopography)"
    ]
   },
@@ -159,6 +160,98 @@
     "print(f\"Minimum elevation = {grid.at_node['topographic__elevation'].min()}\")"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Map data from one grid to another\n",
+    "\n",
+    "In our example we loaded data from an ESRI ASCII file into a ``RasterModelGrid`` but what if we wanted the data to be placed on, for example, a ``HexModelGrid``? *Landlab* doesn' have utilities to do this but we can use [interpolation utilities](https://docs.scipy.org/doc/scipy/reference/interpolate.html) from *scipy* to help us do this."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "First create a `HexModelGrid` with the same number of rows and columns, and grid spacing."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from landlab import HexModelGrid\n",
+    "\n",
+    "hex_grid = HexModelGrid((200, 200), spacing=30.0, node_layout=\"rect\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "To do the mapping, we can use *scipy*'s [LinearNDInterpolator](https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.LinearNDInterpolator.html#scipy.interpolate.LinearNDInterpolator).\n",
+    "\n",
+    "The ``LinearNDInterpolator`` takes two arguments: the first is an array of the *x* and *y* coordinates of the source grid's nodes, and the second the values at those nodes that we want to map. The return value is a ***function*** that maps the values."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from scipy.interpolate import LinearNDInterpolator\n",
+    "\n",
+    "interp = LinearNDInterpolator(grid.xy_of_node, grid.at_node[\"topographic__elevation\"])"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "To interpolate to the destination grid (our `HexModelGrid`) we call our interpolation function with the new grid's node values."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "z = interp(hex_grid.x_of_node, hex_grid.y_of_node)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "If we like, we can then add the values to the hex grid as a data field."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "hex_grid.at_node[\"topographic__elevation\"] = z"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "\n",
+    "hex_grid.imshow(\n",
+    "    \"topographic__elevation\", allow_colorbar=True, vmin=np.nanmin(z), vmax=np.nanmax(z)\n",
+    ")"
+   ]
+  },
   {
    "cell_type": "markdown",
    "metadata": {},

diff --git a/lessons/landlab/landlab/02-data.ipynb b/lessons/landlab/landlab/02-data.ipynb
@@ -18,6 +18,7 @@
     "# Table of Contents\n",
     "* [Attaching data to a grid](#Attaching-data-to-a-grid)\n",
     "  * [Load data from an ESRI ASCII file](#Load-data-from-an-ESRI-ASCII-file)\n",
+    "  * [Map data from one grid to another](#Map-data-from-one-grid-to-another)\n",
     "  * [Load data from OpenTopography](#Load-data-from-OpenTopography)"
    ]
   },
@@ -167,6 +168,98 @@
     "</details>"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Map data from one grid to another\n",
+    "\n",
+    "In our example we loaded data from an ESRI ASCII file into a ``RasterModelGrid`` but what if we wanted the data to be placed on, for example, a ``HexModelGrid``? *Landlab* doesn' have utilities to do this but we can use [interpolation utilities](https://docs.scipy.org/doc/scipy/reference/interpolate.html) from *scipy* to help us do this."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "First create a `HexModelGrid` with the same number of rows and columns, and grid spacing."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from landlab import HexModelGrid\n",
+    "\n",
+    "hex_grid = HexModelGrid((200, 200), spacing=30.0, node_layout=\"rect\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "To do the mapping, we can use *scipy*'s [LinearNDInterpolator](https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.LinearNDInterpolator.html#scipy.interpolate.LinearNDInterpolator).\n",
+    "\n",
+    "The ``LinearNDInterpolator`` takes two arguments: the first is an array of the *x* and *y* coordinates of the source grid's nodes, and the second the values at those nodes that we want to map. The return value is a ***function*** that maps the values."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from scipy.interpolate import LinearNDInterpolator\n",
+    "\n",
+    "interp = LinearNDInterpolator(grid.xy_of_node, grid.at_node[\"topographic__elevation\"])"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "To interpolate to the destination grid (our `HexModelGrid`) we call our interpolation function with the new grid's node values."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "z = interp(hex_grid.x_of_node, hex_grid.y_of_node)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "If we like, we can then add the values to the hex grid as a data field."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "hex_grid.at_node[\"topographic__elevation\"] = z"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "\n",
+    "hex_grid.imshow(\n",
+    "    \"topographic__elevation\", allow_colorbar=True, vmin=np.nanmin(z), vmax=np.nanmax(z)\n",
+    ")"
+   ]
+  },
   {
    "cell_type": "markdown",
    "metadata": {},
@@ -365,6 +458,7 @@
  ],
  "metadata": {
   "anaconda-cloud": {},
+  "celltoolbar": "Tags",
   "kernelspec": {
    "display_name": "Ivy",
    "language": "python",