diff --git a/.readthedocs.yaml b/.readthedocs.yaml
index f80269d6..12f27b85 100644
--- a/.readthedocs.yaml
+++ b/.readthedocs.yaml
@@ -4,13 +4,12 @@ build:
   os: "ubuntu-22.04"
   tools:
     python: "3.12"
-
-python:
-  install:
-    - method: pip
-      path: .
-      extra_requirements:
-        - docs
+  jobs:
+    # Avoid `pip install .[docs]`, which resolves geopandas, rasterio, dask, pvlib,
+    # etc. and often OOMs on RTD builders. Docs only need Sphinx + the source tree.
+    install:
+      - pip install -r docs/requirements.txt
+      - pip install --no-deps .
 
 sphinx:
   configuration: docs/source/conf.py
diff --git a/docs/jupyter_execute/legacy/mask_on_cutout.ipynb b/docs/jupyter_execute/legacy/mask_on_cutout.ipynb
new file mode 100644
index 00000000..b3d138ad
--- /dev/null
+++ b/docs/jupyter_execute/legacy/mask_on_cutout.ipynb
@@ -0,0 +1,444 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Tutorial: Incorporating Mask into Cutout Workflow\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Introduction"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Geodata is able to process geospatial data to extract cutouts over specified geographies. Built off the [rasterio library](https://rasterio.readthedocs.io/en/latest/quickstart.html), the **mask** module imports rasters and shapefiles, merges and flattens multiple layers together, and extracts subsetted cutout data from merged masks and shapefiles.\n",
+    "\n",
+    "After we create a mask, we can incorporate the suitability mask object/file into the Cutout. The cutouts are subsets of data based on specific time and geographic ranges. For more information on the creation of cutout, refer to these tutorials: [Creating Cutouts with MERRA2 Data](https://github.com/east-winds/geodata/blob/master/doc/merra2/merra2_createcutout.md), [Downloading and Creating Cutouts with ERA5 Data](https://github.com/east-winds/geodata/blob/master/doc/era5/era5_download.md)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Setup\n",
+    "\n",
+    "To start, import the geodata package and required libraries."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import matplotlib.pyplot as plt\n",
+    "import xarray as xr\n",
+    "\n",
+    "import geodata"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Download Data"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "We will use a Cutout object created from a downloaded dataset. **If you have already created a cutout, load it here and skip to step 3.**\n",
+    "\n",
+    "\n",
+    "We first download the dataset through `geodata.Dataset()`. In `get_data()`, if we specify `testing=True`, the program downloads only first file in download list (e.g., first day of month)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "dataset_test = geodata.Dataset(\n",
+    "    module=\"merra2\",\n",
+    "    years=slice(2011, 2011),\n",
+    "    months=slice(1, 1),\n",
+    "    weather_data_config=\"slv_radiation_hourly\",\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "if not dataset_test.prepared:\n",
+    "    dataset_test.get_data(testing=True)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Extract the cutout from the trimmed dataset."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "scrolled": true
+   },
+   "outputs": [],
+   "source": [
+    "cutout = geodata.Cutout(\n",
+    "    name=\"china-2011-slv-hourly-test\",\n",
+    "    module=\"merra2\",\n",
+    "    weather_data_config=\"slv_radiation_hourly\",\n",
+    "    xs=slice(73, 136),\n",
+    "    ys=slice(18, 54),\n",
+    "    years=slice(2011, 2011),\n",
+    "    months=slice(1, 1),\n",
+    ")\n",
+    "cutout.prepare()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Load Mask"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "In this tutorial, we use the `china` mask, created in this documentation: [mask_creation_workflow](mask_creation_workflow.ipynb)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# View the contents of the china mask\n",
+    "geodata.mask.load_mask(\"china\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Adding Mask Variables to a Cutout"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Adding Masking Variables"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The `add_mask` method will add attribute `merged_mask` and `shape_mask` from the Mask object to the Cutout object. Once the mask is added to the Cutout object, the `merged_mask` or `shape_mask` from the Mask object will be stored in the format of xarray.DataArray in the Cutout object, and their dimensions will be coarsened to the same dimension with the Cutout metadata.\n",
+    "\n",
+    "The `add_mask` method will look for both `merged_mask` and `shape_mask` attribute saved for the loaded mask, unless the user set the parameter `merged_mask=False`, or `shape_mask=False`.\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "scrolled": false
+   },
+   "outputs": [],
+   "source": [
+    "cutout.add_mask(\"china\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Plot the merged mask, coarsened to cutout resolution"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "cutout.merged_mask.plot()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Adding Area Variable\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "To calculate and add the variation of grid cell areas by latitude to the cutout, use the `add_grid_area` method. Keeping track of the area for each grid cell is necessary for analyses such as calculating the weighted sum of the grid cells based on their area."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "cutout.add_grid_area()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Creating PV Data Through Cutout Conversion"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The code block below will use the `geodata.convert.pv` method to generate `ds_cutout`, an xarray Dataset that contains the pv variable for the cutout.\n",
+    "\n",
+    "We transform the xarray DataArray into a xarray DataSet (which can contain multiple DataArray). "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "ds_cutout = geodata.convert.pv(cutout, panel=\"KANEKA\", orientation=\"latitude_optimal\").to_dataset(\n",
+    "    name=\"solar\"\n",
+    ")\n",
+    "len(ds_cutout.time)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "We also need to remove the time dimension by calculating daily means via `ds_cutout.coarsen(time=24, boundary=\"exact\").mean()`, which aggregates the values over its 24 timestamps."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "ds_cutout_mean = ds_cutout.coarsen(time=24, boundary=\"exact\").mean()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Combining PV Data with Mask"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The `mask` method for the Cutout will mask converted xarray.Dataset variable, such as `ds_cutout` and `ds_cutout_mean` created above, by combining it with merged_mask or shape_mask in the Cutout object. It will return a dictionary of xarray Dataset. Each key in the dictionary is one unique mask from either the merged_mask or shape_mask variable from the Cutout object, and each value is an xarray dataset containing the dataSet variable (`ds_cutout` or `ds_cutout_mean`) with the mask and area values.\n",
+    "\n",
+    "The program will automatically search for `merged_mask` and `shape_mask` to combine with the xarray.Dataset, unless the user specify `merged_mask=False` or `shape_mask=False`. The masks in `shape_mask` will have the same key as it has in the `shape_mask` attribute, and the mask for `merged_mask` will have the same key name `merged_mask`, as `merged_mask` is unique to each mask."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "#### Daily averaged PV values"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "scrolled": false
+   },
+   "outputs": [],
+   "source": [
+    "ds_mask_mean = cutout.mask(dataset=ds_cutout_mean)\n",
+    "ds_mask_mean.keys()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "From the output variable `ds_mask_mean`, check out the combined xarray.Dataset for the Jiangsu province, and plot each of its xarray.DataArray."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "ds_mask_mean[\"Jiangsu\"]"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Visualize the averaged PV value for each grid cell in the Cutout. Note that the data is the aggregated value for the date."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "ds_mask_mean[\"Jiangsu\"][\"solar\"].plot()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Visualize the masking value for each grid cell in the Cutout."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "ds_mask_mean[\"Jiangsu\"][\"mask\"].plot()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "#### Area and Mask-Weighted Hourly PV Values\n",
+    "\n",
+    "We use the raw hourly output generated by cutout to create time-series PV plots weighted by the mask and area. Note that we transposed ds_cutout so that time is set as the first dimension, which ease the following calculation since we want to aggregate the array spatially from each grid cell."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "ds_mask = cutout.mask(ds_cutout)\n",
+    "ds_mask.keys()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Calculate the aggregated mean solar PV for each provinces, at each time point. We will apply this equation below to calculate the area-weighted average. We save the result into a dictionary `PV_dict`, where its keys are the provinces, and the corresponding values are the PV series.\n",
+    "\n",
+    "$$\\text{Aggregated Solar Power For Each Region} = \\frac{\\sum_{}^\\text{For Each Grid Cell}\\text{Grid Cell Area} \\times \\text{Mask Value} \\times \\text{Solar Power}}{\\sum_{}^\\text{For Each Grid Cell}\\text{Grid Cell Area} \\times \\text{Mask Value}}$$"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "PV_dict = {}\n",
+    "\n",
+    "for prov_name in list(ds_mask)[1:]:\n",
+    "    PV_dict[prov_name] = (\n",
+    "        (ds_mask[prov_name][\"solar\"] * ds_mask[prov_name][\"mask\"] * ds_mask[prov_name][\"area\"])\n",
+    "        .sum(axis=1)\n",
+    "        .sum(axis=1)\n",
+    "    ) / (ds_mask[prov_name][\"mask\"] * ds_mask[prov_name][\"area\"]).sum()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The aggregated PV time-series for Zhejiang province."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "PV_dict[\"Zhejiang\"]"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Finally, for each province, plot the solar series weighted by mask * area."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "for prov_name, series in PV_dict.items():\n",
+    "    plt.plot(series, label=prov_name)\n",
+    "\n",
+    "    plt.title(f\"Solar series weighted by area for Chinese provinces.\")\n",
+    "    plt.grid()\n",
+    "    plt.legend()\n",
+    "    plt.xlabel(\"2011-01-01 Hour\")\n",
+    "    plt.ylabel(\"Aggregated weighted PV value for suitable area\")"
+   ]
+  }
+ ],
+ "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.10.11"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}
\ No newline at end of file
diff --git a/docs/jupyter_execute/legacy/merra2/merra2.ipynb b/docs/jupyter_execute/legacy/merra2/merra2.ipynb
new file mode 100644
index 00000000..92b1cf13
--- /dev/null
+++ b/docs/jupyter_execute/legacy/merra2/merra2.ipynb
@@ -0,0 +1,615 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# MERRA2 Analysis Process\n",
+    "\n",
+    "This Jupyter notebook provides a brief overview of how to use the **geodata** package to download MERRA2 climate data, create geographic-temporal subsets called cutouts, and use those cutouts to generate standalone datasets for separate analysis.\n",
+    "\n",
+    "*The following guide assumes you have installed and configured **geodata** and all required dependencies.*"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Step 1 - Setup\n",
+    "\n",
+    "Import the package first."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import geodata"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Notifications in **geodata** are implemented using `loggers` from the `logging` library.\n",
+    "It is recommended to always launch a logger to get information on what is going on. For debugging, you can use the more verbose `level=logging.DEBUG`:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import logging\n",
+    "\n",
+    "logging.basicConfig(level=logging.INFO)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Step 2 - Download\n",
+    "\n",
+    "Assuming you have previously created an Earthdata Login profile and approved the GES DISC app, you can download MERRA2 data from the source as follows.\n",
+    "\n",
+    "First, define a dataset object for the data you wish to download:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "DS = geodata.Dataset(\n",
+    "    module=\"merra2\",\n",
+    "    weather_data_config=\"surface_flux_monthly\",\n",
+    "    years=slice(2010, 2010),\n",
+    "    months=slice(1, 7),\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "* Use `module` to specify the data source. In this example, it is \"merra2\".\n",
+    "* Use `weather_data_config` to specifiy the dataset.  In this example, it is the [MERRA2 monthly mean, single-level surface flux diagnostics](https://disc.gsfc.nasa.gov/datasets/M2TMNXFLX_5.12.4/summary)\n",
+    "    * To download the [MERRA2 hourly, single-level surface flux diagnostics](https://disc.gsfc.nasa.gov/datasets/M2T1NXFLX_5.12.4/summary), specify `weather_data_config = \"surface_flux_hourly\"`.\n",
+    "* Use `years=slice()` and `months=slice()` to specify the years and months for download.  In each parameter, the first value indicates the start period, and the second value the end period.\n",
+    "\n",
+    "Use the code block below to begin the download."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "When a `dataset` object is created, **geodata** performs a check to see if the data specified has already been downloaded by checking for the existence of MERRA2 datafiles in the `merra2` directory configured in `src/geodata/config.py` (downloaded data is placed into subdirectories by year and then - for daily files - by month, ie `2011/01, 2011/02, 2012/01`, etc).  Monthly files are simply placed in the month's folder.  If downloaded data is found, the `prepared` attribute is set to `True` upon `dataset` object declaration.\n",
+    "\n",
+    "Accordingly, the snippet below saves you the trouble of accidentally redownloading data if it is already present in the correct subdirectories."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "if DS.prepared == False:\n",
+    "    DS.get_data()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Finally, in order to use the downloaded MERRA2 data with **geodata**, run:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "DS.trim_variables()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "`trim_variables()` subsets and resaves the downloaded files so that only those variables needed to generate **geodata** outputs are kept."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    " "
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Step 3 - Create Cutout\n",
+    "\n",
+    "A cutout is a subset of downloaded data based on specified time periods and geographic coordinates.  Cutouts are saved to the cutout directory specified in `src/geodata/config.py` and can be used to generate multiple outputs.\n",
+    "\n",
+    "*Note: 04/02/2020 - There is a known issue with MERRA2-based cutouts where running `cutout.prepare(overwrite=True)` on an existing cutout prevents the cutout from being used to generate outputs.  A workaround is to manually delete the problem cutout and recreate it from scratch.  A fix is planned pending investigation."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "To create a cutout, run the following:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "cutout = geodata.Cutout(\n",
+    "    name=\"tokyo-2010-test\",\n",
+    "    module=\"merra2\",\n",
+    "    weather_data_config=\"surface_flux_monthly\",\n",
+    "    xs=slice(138.5, 139.5),\n",
+    "    ys=slice(35, 36),\n",
+    "    years=slice(2010, 2010),\n",
+    "    months=slice(7, 7),\n",
+    ")\n",
+    "cutout.prepare()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The above code creates a cutout for July 2010 for a geographic area roughly corresponding to the Tokyo metropolitan area. Walking through the parameters:\n",
+    "\n",
+    "* `name` will be the name of the directory created in the cutouts folder where **geodata** will place the data files corresponding to the cutout.\n",
+    "* `module` indicates the source for the data from which the cutout is created.\n",
+    "* `weather_data_config` indicates the specific dataset from the source.  For MERRA2, the available options are `surface_flux_hourly` and `surface_flux_monthly`.\n",
+    "* Use `xs=slice()` and `ys=slice()` to define a geographical range for the cutout.\n",
+    "* Use `years=slice()` and `months=slice()` to define a temporal range for the cutout.  Naturally, the indicated time range must be present within the source data.\n",
+    "\n",
+    "`geodata.Cutout()` only defines the cutout object in memory.  To actually create the cutout files, run `prepare()`.  \n",
+    "As with `get_data()`, `prepare()` will first perform a check to see if a cutout has already been created at the same specified, and will exit the creation process if a cutout already exists.  To override this behavior and force a recalculation of the cutout, run `prepare(overwrite=True)`."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "To verify the results of the cutout, you can print some attributes to the console as follows.\n",
+    "\n",
+    "Basic information:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "cutout"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Name:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "cutout.name"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Coordinates:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "cutout.coords"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "All metadata:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "cutout.meta"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Information about the variable config used to download the data:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "cutout.dataset_module.weather_data_config"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "For Merra2, you can confirm variables downloaded this way:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "cutout.dataset_module.weather_data_config[\"surface_flux_monthly\"][\"variables\"]"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Step 4 - Generate Outputs\n",
+    "\n",
+    "**geodata** currently supports the following wind outputs using MERRA2 surface flux diagnostic data.\n",
+    "* Wind generation time-series (`wind`)\n",
+    "* Wind speed time-series (`windspd`)\n",
+    "* Wind power density time-series (`windpwd`)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Wind Generation Time-series\n",
+    "Convert wind speeds for turbine to wind energy generation using the following code:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "ds_wind = geodata.convert.wind(cutout, turbine=\"Suzlon_S82_1.5_MW\", smooth=True, var_height=\"lml\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Going over the parameters:\n",
+    "\n",
+    "* `cutout` - **string** -  A cutout created by `geodata.Cutout()`\n",
+    "* `turbine` - **string or dict** - Name of a turbine known by the reatlas client or a turbineconfig dictionary with the keys 'hub_height' for the hub height and 'V', 'POW' defining the power curve.  For a full list of currently supported turbines, see [the list of Turbines here.](https://github.com/east-winds/geodata/tree/master/geodata/resources/windturbine)\n",
+    "* `smooth` - **bool or dict** - If True smooth power curve with a gaussian kernel as determined for the Danish wind fleet to Delta_v = 1.27 and sigma = 2.29. A dict allows to tune these values.\n",
+    "\n",
+    "*Note* - \n",
+    "You can also specify all of the general conversion arguments documented in the `convert_and_aggregate` function (e.g. `var_height='lml'`)."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The convert function returns an xarray dataset, which is an in-memory representation of a NetCDF file."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "ds_wind"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "To convert this array to a more conventional dataframe, run:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "df_wind = ds_wind.to_dataframe(name=\"wind\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "which converts the xarray dataset into a pandas dataframe:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "df_wind"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "To output the data to a csv for separate analysis:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "df_wind.to_csv(\"merra2_wind_data.csv\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Extract wind speeds at given height (ms-1)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "ds_windspd = geodata.convert.windspd(cutout, turbine=\"Vestas_V66_1750kW\", var_height=\"lml\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Going over the parameters:\n",
+    "\n",
+    "* `cutout` - **string** -  A cutout created by `geodata.Cutout()`\n",
+    "* `**params` - Must have 1 of the following:\n",
+    "    - `turbine` - **string or dict** - Name of a turbine known by the reatlas client or a turbineconfig dictionary with the keys 'hub_height' for the hub height and 'V', 'POW' defining the power curve.  For a full list of currently supported turbines, see [the list of Turbines here.](https://github.com/east-winds/geodata/tree/master/geodata/resources/windturbine)\n",
+    "    - `hub-height` - **num** - Extrapolation height (m)\n",
+    "    \n",
+    "*Note* - \n",
+    "You can also specify all of the general conversion arguments documented in the `convert_and_aggregate` function (e.g. `var_height='lml'`)."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The convert function returns an xarray dataset, which is an in-memory representation of a NetCDF file."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "ds_windspd"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "To convert this array to a more conventional dataframe, run:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "df_windspd = ds_windspd.to_dataframe(name=\"windspd\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "which converts the xarray dataset into a pandas dataframe:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "df_windspd"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "To output the data to a csv for separate analysis:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "df_windspd.to_csv(\"merra2_windspd_data.csv\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Wind Power Density Time-series\n",
+    "\n",
+    "Extract wind power density at given height, according to:\n",
+    "**WPD = 0.5 * Density * Windspd^3**"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "ds_windwpd = geodata.convert.windwpd(cutout, turbine=\"Vestas_V66_1750kW\", var_height=\"lml\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Going over the parameters:\n",
+    "\n",
+    "* `cutout` - **string** -  A cutout created by `geodata.Cutout()`\n",
+    "* `**params` - Must have 1 of the following:\n",
+    "    - `turbine` - **string or dict** - Name of a turbine known by the reatlas client or a turbineconfig dictionary with the keys 'hub_height' for the hub height and 'V', 'POW' defining the power curve.  For a full list of currently supported turbines, see [the list of Turbines here.](https://github.com/east-winds/geodata/tree/master/geodata/resources/windturbine)\n",
+    "    - `hub-height` - **num** - Extrapolation height (m)\n",
+    "    \n",
+    "*Note* - \n",
+    "You can also specify all of the general conversion arguments documented in the `convert_and_aggregate` function (e.g. `var_height='lml'`)."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The convert function returns an xarray dataset, which is an in-memory representation of a NetCDF file."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "ds_windwpd"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "To convert this array to a more conventional dataframe, run:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "df_windwpd = ds_windwpd.to_dataframe(name=\"windwpd\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "which converts the xarray dataset into a pandas dataframe:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "df_windwpd"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "To output the data to a csv for separate analysis:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "df_windwpd.to_csv(\"merra2_windwpd_data.csv\")"
+   ]
+  }
+ ],
+ "metadata": {
+  "file_extension": ".py",
+  "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.4"
+  },
+  "mimetype": "text/x-python",
+  "name": "python",
+  "npconvert_exporter": "python",
+  "pygments_lexer": "ipython3",
+  "version": 3
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}
\ No newline at end of file
diff --git a/docs/jupyter_execute/mask/mask_creation_workflow.ipynb b/docs/jupyter_execute/mask/mask_creation_workflow.ipynb
new file mode 100644
index 00000000..9607f7a5
--- /dev/null
+++ b/docs/jupyter_execute/mask/mask_creation_workflow.ipynb
@@ -0,0 +1,1151 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Tutorial: Typical Mask Creation Workflow"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Introduction"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Geodata is able to process geospatial data to extract cutouts over specified geographies. Built off the [rasterio library](https://rasterio.readthedocs.io/en/latest/quickstart.html), the **mask** module imports rasters and shapefiles, merges and flattens multiple layers together, and extracts subsetted cutout data from merged masks and shapefiles.\n",
+    "\n",
+    "Functionalities explored in this notebook:\n",
+    "\n",
+    "- [Creating a mask object, adding and manipulating layers](#mask-creation--adding-and-manipulating-layers)\n",
+    "- [Opening a shapefile and adding shape features as layers](#adding-shape-features-as-a-layer)\n",
+    "- [Merging and flattening layers](#merging-and-flattening-layers)\n",
+    "- [Eliminate small contiguous areas](#eliminate-small-contiguous-areas)\n",
+    "- [Extracting shapes from mask](#extracting-shapes-from-masks)\n",
+    "- [Saving and loading masks](#saving-and-loading-masks)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Setup"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "To start, import the geodata package and required libraries. We can also import the `geodata.mask.show()` method for simplicity of its use."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import geopandas as gpd\n",
+    "import matplotlib.pyplot as plt\n",
+    "import numpy as np\n",
+    "import pandas as pd\n",
+    "\n",
+    "import geodata\n",
+    "from geodata.mask import show"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Additionally, we use [cartopy](https://scitools.org.uk/cartopy/docs/latest/tutorials/using_the_shapereader.html#cartopy.io.shapereader.Reader) to download some common administrative region shapes, but user-provided shapefiles will also work:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import cartopy.io.shapereader as shpreader"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Shapefiles and Rasters\n",
+    "\n",
+    "We will use the following geotiff and shape files for this demo:\n",
+    "\n",
+    "\n",
+    "- `china_modis.tif`\n",
+    "\n",
+    "    We downloaded the MODIS land cover data, which uses satellite remote sensing data to estimate the land use type on an annual basis. See: [EarthData_MCD12Q1](https://lpdaac.usgs.gov/products/mcd12q1v006/).\n",
+    "\n",
+    "    We will use the IGBP classification ('LC_Type1') which has 17 different land use characterizations (the corresponding data thus takes values from 1.0 to 17.0).\n",
+    "    All the \"Bands\" are listed here: [Google_earth_engine_MODIS_006_MCD12Q1](https://developers.google.com/earth-engine/datasets/catalog/MODIS_006_MCD12Q1#bands)\n",
+    "\n",
+    "- `china_elevation.tif` and `china_slope.tif`\n",
+    "\n",
+    "    These two rasters are based on the elevation map from: [Google_earth_engine_MODIS_CGIAR_SRTM90_V4](https://developers.google.com/earth-engine/datasets/catalog/CGIAR_SRTM90_V4?hl=en). Slope was computed in degrees using the 4-connected neighbors of each pixel. \n",
+    "\n",
+    "\n",
+    "- `UNEP_WDPA_China` Shapefiles\n",
+    "\n",
+    "    We downloaded the environmental protected area for China from: [ProtectedPlanet_China](https://www.protectedplanet.net/country/CHN). These shapefiles are distributed among 3 subfolders upon successful download and decompression due to the large size. We will create path variables for all three subfolders and we will only take the polygon shapes.\n",
+    "\n",
+    "    Alternatively, We can also retrieve the environmental protected area from Google Earth Engine: [Google_earth_engine_WCMC_WDPA](https://developers.google.com/earth-engine/datasets/catalog/WCMC_WDPA_current_polygons). The shapefile will contain the protected shapes from entire world (and the size is slightly over 1 GB), and additional data cleaning will be necessary if the user wants just the shapes within China. "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "modis_path = \"data/china_modis.tif\"\n",
+    "elevation_path = \"data/china_elevation.tif\"\n",
+    "slope_path = \"data/china_slope.tif\"\n",
+    "\n",
+    "wdpa_shape_path_0 = \"data/shapefiles/0/WDPA_WDOECM_Nov2021_Public_CHN_shp-polygons.shp\"\n",
+    "wdpa_shape_path_1 = \"data/shapefiles/1/WDPA_WDOECM_Nov2021_Public_CHN_shp-polygons.shp\"\n",
+    "wdpa_shape_path_2 = \"data/shapefiles/2/WDPA_WDOECM_Nov2021_Public_CHN_shp-polygons.shp\""
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Let us get province shapes from `cartopy` and save the path as `prov_path`. This can also be the path to user-supplied shape files."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "prov_path = shpreader.natural_earth(\n",
+    "    resolution=\"10m\", category=\"cultural\", name=\"admin_1_states_provinces\"\n",
+    ")\n",
+    "prov_path"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Load the shapes contained in path `prov_path` using the `geopandas` library."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "scrolled": true
+   },
+   "outputs": [],
+   "source": [
+    "all_shapes = gpd.read_file(prov_path, encoding=\"utf-8\")\n",
+    "all_shapes.head(2)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "GeoPandas data filtering with GeoDataFrame is identical to pandas. Let us select all the rows that contains shape within China."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "china_shapes = all_shapes[all_shapes[\"admin\"] == \"China\"]"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Next, to load the WDPA environmental protected shapefiles as a layer in the china mask, we will use the GeoPandas library. `gpd.read_file()` will return a GeoPandas dataframe including shape attributes and geometry given the file path. Like Pandas, we can read multiple dataframes and concat them together. In the code below, we will create one GeoPandas dataframe from three paths that we have for the Chinese environmental protected shapes."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "wdpa_shapes = pd.concat([\n",
+    "    gpd.read_file(wdpa_shape_path_0),\n",
+    "    gpd.read_file(wdpa_shape_path_1),\n",
+    "    gpd.read_file(wdpa_shape_path_2)\n",
+    "])\n",
+    "wdpa_shapes.head(2)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Mask Creation & Adding and Manipulating Layers"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The mask object consists of multiple layers and manipulations performed on them. To add a layer, the four methods below perform same functions. A user may add a layer to the mask by specifying paths when a new instance is created, or use the `add_layer` method. We will add the following two files: `china_elevation.tif`, and `china_modis.tif` to the `China` mask, and name them `elevation` and `modis` layers."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "scrolled": true
+   },
+   "outputs": [],
+   "source": [
+    "# Method 1: Initialize one layer, add one layer\n",
+    "china = geodata.Mask(\"China\", layer_path=elevation_path)\n",
+    "china.rename_layer(\"china_elevation\", \"elevation\")\n",
+    "china.add_layer(modis_path, layer_name=\"modis\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Method 2:  Initialize empty, add two layers using dict\n",
+    "china = geodata.Mask(\"China\")\n",
+    "china.add_layer(layer_path={\"elevation\": elevation_path, \"modis\": modis_path})"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Method 3:  Initalize with two layers passed as list\n",
+    "china = geodata.Mask(\n",
+    "    \"China\", layer_path=[elevation_path, modis_path], layer_name=[\"elevation\", \"modis\"]\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Method 4:  Initialize with two layers passed as dict\n",
+    "china = geodata.Mask(\n",
+    "    \"China\", layer_path={\"elevation\": elevation_path, \"modis\": modis_path}\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Display the mask object in the jupyter notebook:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "china"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Each mask object has several attributes:\n",
+    "- `layers`: a dictionary of name (key)  - rasterio file opener (values). The <\\open DatasetReader> can be the input for many other mask methods for the module.  \n",
+    "- `merged_mask`: the merged and flatten mask of its layers, the merged raster from `layers`\n",
+    "- `shape_mask`: similar to the `layers` attribute, but a dictionary of extracted shapes from the merged mask by default. Users may also extracted shape masks from specified layers in `self.layers`.\n",
+    "- `saved`: whether this mask object has been saved locally.\n",
+    "- `mask_dir`: the directory to save the mask object, by default it should be the mask dir in config.py."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Show the `slope` layer in mask `china`. The `show` method will always try to show the proper latitude and longitude, unless we call it `show(layer, lat_lon = False)`."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "china.layers[\"elevation\"]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "show(china.layers[\"elevation\"], title=\"Elevation of China in meters\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "china.layers"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "**Some useful methods to examine the layers**\n",
+    "\n",
+    "- `china.get_res()`: get resolution of each layer, in lat-lon coordinates\n",
+    "- `china.get_res(product = True)`: get grid cell size, in product of lat-lon coordinate differences\n",
+    "- `china.get_bounds()`: get bounds, in lat-lon coordinates"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "china.get_bounds()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Note that the modis layer has a very different bounding box then the slope layer in lat-lon coordinate system. This is because the modis layer was converted to the lat-lon CRS from a different CRS when it was added to the object. The following section will explore CRS conversion."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### CRS conversion, trimming, and cropping (Optional)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Method `open_tif` can open a layer without adding it to the layer, this allows us to visualize it before-hand. It is a good practice to close the raster after opening it to avoid writing permission conflict issues. Closing the raster below does not involve any layer operation associated with the mask object. "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "modis_opener = geodata.mask.open_tif(modis_path, show_raster=True)\n",
+    "modis_opener.close()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "We can use `remove_layer` method to remove a layer to mask `china`. This method will properly close the raster file, because the raster file would remain open after being added to the mask."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "china.remove_layer(\"modis\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The `add_layer` method incorporates coordinate reference system (CRS) conversion to lat-lon (EPSG:4326), if necessary. Note that this method will overwrite the layer by default, if it is in the object already, unless the user specifies `replace=False`. \n",
+    "\n",
+    "The method will automatically trim the all-zero columns/rows. By default, the paramater `trim` is set to `True`. If we do not set it to True, we might generate a converted raster with new CRS but many all-zero columns and rows."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "china.add_layer(modis_path, \"modis\", trim=False)\n",
+    "show(china.layers[\"modis\"], title=\"China Modis CRS converted (No trimming)\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "We can also crop a raster/layer with user-defined dimensions: method `crop_layer` can take either starting indices of top/left, ending indices of right/bottom, or coordinates values in lat/long to trim the raster.\n",
+    "\n",
+    "The difference between `crop_layer` and `trim_layer` is that `crop_layer` must take in user specified range to crop the raster, and `trim_layer` would remove the all zero rows and columns automatically for a raster. So that if the user do not know which index to start and end to remove the empty rows/columns, `trim_raster` is better.\n",
+    "\n",
+    "The method `crop_raster` (`geodata.mask.crop_raster`) is similar to `crop_layer` but can take a layer name as input, so that the user does not need to add a raster as a layer to call that method. (Similar method: `trim_layer`/`trim_raster`, `binarize_layer`/`binarize_raster`)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "china.crop_layer(\"modis\", bounds=(73, 17, 135, 54))\n",
+    "show(china.layers[\"modis\"], title=\"China Modis Layer Cropped\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "This performs the same function by passing the layer to `crop_raster`:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "china.layers[\"modis\"] = geodata.mask.crop_raster(\n",
+    "    china.layers[\"modis\"], (73, 17, 135, 54)\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Filter a layer"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The mask module also supports filtering a layer based on list of categorical values, a minimum (lower) boundary, or maximum (upper) boundary.\n",
+    "\n",
+    "In the `filter_raster` method, a user may specify any of the `value` (the list of numberic values in the raster array to be selected), `max_bound`, and `min_bound` parameters to selected desired values. If the parameter `binarize` is False (by default), the method will return the original values of the raster that satisfy the conditions, otherwise the method will return 1 for the values that satisfy the conditions and 0 elsewhere."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "#### Select Categorical Values from MODIS Layer\n",
+    "\n",
+    "Since the modis layer has 17 distinct values for different land use types, we want to create a layer of binary values, indicating unavailable land as 0, and available land as 1.\n",
+    "\n",
+    "We wish to create a mask where :\n",
+    "\n",
+    "- all forested areas (values 1-5) are 0 (i.e., unsuitable)\n",
+    "- all urban areas (13) are 0\n",
+    "- all others are 1\n",
+    "\n",
+    "\n",
+    "Let us use method `filter_raster` to create a layer of `modis_filtered` binary mask, where 1, 2, 3, 4, 5, and 13 will be unavailable land assigned 0 and the rest of the values will be 1 (available).\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "avail_values = list(set(range(1, 18)) - {1, 2, 3, 4, 5, 13})\n",
+    "avail_values"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "china.layers[\"modis_filtered\"] = geodata.mask.filter_raster(\n",
+    "    china.layers[\"modis\"], binarize=True, values=avail_values\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "scrolled": false
+   },
+   "outputs": [],
+   "source": [
+    "china.remove_layer(\"modis\")\n",
+    "show(china.layers[\"modis_filtered\"])"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "#### Filter elevation layer\n",
+    "\n",
+    "Because we cannot build renewable energy in areas with high elevation, let us set the constraint from the `elevation` layer, by using elevation < 4000m at 1 and other areas as 0. The result layer `elevation_filtered` will have only 1 and 0 as unique values."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "china.filter_layer(\n",
+    "    \"elevation\", dest_layer_name=\"elevation_filtered\", max_bound=4000, binarize=True\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "china.remove_layer(\"elevation\")\n",
+    "show(china.layers[\"elevation_filtered\"])"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "#### Filter Slope Layer\n",
+    "\n",
+    "We also cannot build renewable energy in area with large slopes, so let us set another constraint from the `slope` layer from the slope tif file, by using slope < 20 degree at 1 and else as 0. The result layer `slope_filtered` will have only 1 and 0 as unique values."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "First, add the slope raster to the china mask."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "china.add_layer(slope_path, layer_name=\"slope\")\n",
+    "show(china.layers[\"slope\"])"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Filter the raster, delete the old slope layer."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "china.filter_layer(\n",
+    "    \"slope\", dest_layer_name=\"slope_filtered\", max_bound=20, binarize=True\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "scrolled": false
+   },
+   "outputs": [],
+   "source": [
+    "china.remove_layer(\"slope\")\n",
+    "show(china.layers[\"slope_filtered\"])"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Additional Visualization Options"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "We can plot the provinces on a selected layer by taking `shape` input in the `show()` method. Here, we will use the `china_shapes` that we obtained from `all_shape`. Its `geometry` column is a Series of shapes (shapely.geometry or MultiPolygon) for Chinese provinces."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "show(china.layers[\"modis_filtered\"], shape=china_shapes[\"geometry\"])"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Adding Shape Features as a Layer"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Recall that we have previously loaded the environmental protected shapes of China in a GeoPandas dataframe."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "len(wdpa_shapes)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The three shapefiles have 78 features altogether, but we want to add all the features to one new layer instead of 78 new layers. The input shape should be a python dictionary, where there is a key for each unique shape. Also, in the `add_shape_layer` method, we will specify a `combine_name` to combine the features into one layer in this case, since we want the mask to have just one more layers, not 78 more layers.\n",
+    "\n",
+    "When adding a shapefile, we must specify the dimensions. We will also use `reference layer = 'slope_filtered'` so the new shape layer will have the same dimension with the `slope_filtered` layer. If the mask is empty and does not contain any layer, the user will have to specify the `resolution` parameter for the raster layer dimension.\n",
+    "\n",
+    "By default, this method will have paramater `exclude` that defaults to `False`. When it is true, area inside the shape is 0. When it is false, area inside the shape is 1. In this use case, however, we want 0 for area inside of the shape as they are environmental protected areas to exclude. We can just use the default method call."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "china.add_shape_layer(\n",
+    "    wdpa_shapes[\"geometry\"].to_dict(),\n",
+    "    reference_layer=\"slope_filtered\",\n",
+    "    combine_name=\"protected\",\n",
+    ")\n",
+    "show(\n",
+    "    china.layers[\"protected\"],\n",
+    "    title=\"WDPA Protected area shape features as a new layer\",\n",
+    "    grid=True,\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "We can also use the parameter `buffer` in `add_shape_layer` method to create an approximate representation of all locations within a given (perpindicular) distance of the shape object. The units for the buffer are given in kilometers.\n",
+    "\n",
+    "Note that since the units of the original shape are in lat-lon coordinates, when we add the buffer, we will need to have a CRS that has meter as unit. The program will convert the shapes to that CRS, add the buffer around shapes, then convert it back to the lat-lon CRS system. By default, we used \"EPSG:6933\", an equal area projection CRS to add buffer in kilometer. "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "km_buffer = 20\n",
+    "\n",
+    "china.add_shape_layer(\n",
+    "    wdpa_shapes[\"geometry\"].to_dict(),\n",
+    "    reference_layer=\"slope_filtered\",\n",
+    "    combine_name=\"protected_with_buffer\",\n",
+    "    buffer=km_buffer,\n",
+    ")\n",
+    "\n",
+    "show(\n",
+    "    china.layers[\"protected_with_buffer\"],\n",
+    "    title=f\"WDPA Protected area shape with {km_buffer}km buffer\",\n",
+    "    grid=True,\n",
+    ")\n",
+    "\n",
+    "china.remove_layer(\"protected_with_buffer\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Merging and Flattening Layers"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "In order to combine all layers into one, we use the `merge_layer` method which creates a new layer called `merged_mask`. This merges multiple layers together and flattens them using either **and** (default) or **sum** method, saving the result to `self.merged_mask` by default. Geospatial bounds and resolution of the output layer are in the units of the input file coordinate reference system, but by default, we will use the resolution of the layer with the best (finest) resolution for the output bounds/resolution, unless a reference layer is provided. In this case, the resolution of the merged_mask is the same with the `modis_filtered` layer. "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "china.get_res()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "china.merge_layer(attribute_save=False, show_raster=False).res"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Binary `AND` Method"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "By default, the `merge_layer` method will use a binary 'and' method: for each grid cell, if any of the n layers are 0, then the returned `self.merged_layer` will also have 0 at that location. In other words, if all the layers indicate that a land is available (!=0), the merged result will have value 1.\n",
+    "\n",
+    "`merge_layer` may also take in an optional parameter `layers`, which is a list of layer names stored in the object, if the user does not wish to merge all layers in the object. If the user does not want to save the result to the `merged_mask` attribute, the user can specify `attribute_save = False`."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# merge and plot only, do not save\n",
+    "china.merge_layer(attribute_save=False, layers=[\"slope_filtered\", \"modis_filtered\"])"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "china"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Try again with the `reference_layer` parameter:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "china.merge_layer(\n",
+    "    layers=[\"elevation_filtered\", \"modis_filtered\"],\n",
+    "    reference_layer=\"elevation_filtered\",\n",
+    "    show_raster=False,\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The result of the `merged_mask` method is saved to `china.merged_mask` with the same resolution as the reference layer, in this case `elevation_filtered`."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "china.merged_mask.res"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "For the purpose of this demonstration, we will select the `AND` method for the final merged_mask. We can also trim the border of the merged mask since the 4 layers have different boundaries. We can set the parameter `trim = True`."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "china.merge_layer(trim=True)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### `SUM` Method"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The sum method will add up the values from all the layers using weights. When there is no weight dict provided, all the layers for merging will have weights of 1 by default.\n",
+    "\n",
+    "Note: since we are not using the sum method to proceed to the following sections, we will keep `attribute_save = False` to prevent this method from overwriting the mask we have previously created above."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "china.merge_layer(method=\"sum\", attribute_save=False, trim=True)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "This distribution is completely arbitrary for the purpose of demonstration of the module: (Note: The weights do not need to have a total of 1)\n",
+    "- elevation_filtered: 0.15, slope_filtered: 0.1, modis_filtered: 0.3, protected: 0.45\n",
+    "\n",
+    "We will write the result to a new variable `customized_merged_layer` for continuing processing."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "customized_merged_layer = china.merge_layer(\n",
+    "    method=\"sum\",\n",
+    "    weights={\n",
+    "        \"elevation_filtered\": 0.15,\n",
+    "        \"slope_filtered\": 0.1,\n",
+    "        \"modis_filtered\": 0.3,\n",
+    "        \"protected\": 0.45,\n",
+    "    },\n",
+    "    attribute_save=False,\n",
+    "    trim=True,\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "If the continuous value created by `merged_mask` represents a suitability metric, we could set a minimum value of 0.8 to be considered \"suitable\" (or 1). We then apply the `filter_raster` method on the merged layer."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "customized_merged_layer = geodata.mask.filter_raster(\n",
+    "    customized_merged_layer, min_bound=0.8, binarize=True\n",
+    ")\n",
+    "show(customized_merged_layer)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Eliminate Small Contiguous Areas"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Using the above methods, we might end up with many small contiguous areas that are marked suitable but surrounded by an unsuitable region. We may want to exclude such regions from renewable energy development. The `filter_area` method will remove the small contiguous suitable regions by transforming the merged mask raster to polygons/shapes, calculating the area of each polygon, and filtering out polygons that are smaller than a given threshold. Units are given in kilometer-squared (km$^2$).\n",
+    "\n",
+    "By default, `filter_area` uses the merged mask raster and returns a new raster, unless input/output layers are specified by `layer_name` and `dest_layer_name`. \n",
+    "\n",
+    "By default, its `shape_value` parameter is 1, indicating that we are only interested in finding all groups of cells with value 1 (suitable) for elimination. We specify the threshold with the `min_area` parameter.\n",
+    "\n",
+    "Note: the `filter_area` method may take a long time (5 or more minutes depending on the complexity of your layer and your computational setup). The method relies upon `rasterio.rasterize`, see performance notes: https://rasterio.readthedocs.io/en/latest/api/rasterio.features.html#rasterio.features.rasterize\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "For example, if we focus on Guangdong province in Southern China from the merged mask, we notice that there are many small islands in the ocean that are marked as suitable areas. We want to exclude these small regions from our merged mask."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "plt.imshow(china.merged_mask.read(1)[4800:5300, 5700:6600], interpolation=\"none\")\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Call `filter_area` to remove all contiguous suitable region shapes smaller than 100 km$^2$:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "china.merged_mask = geodata.mask.filter_area(china, min_area=100)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "There shapes are removed in the new merged_mask."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "plt.imshow(china.merged_mask.read(1)[4800:5300, 5700:6600], interpolation=\"none\")\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Extracting Shapes from Masks"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Sometimes the user needs to generate masks and perform analysis for a collection of regions (e.g., at the state/province level). The purpose of shape extraction (`extract_shapes`) is to separate `merged_mask` values for each region, with the result a dictionary of name-mask pairs in the `shape_mask` attribute of the mask object. The values of `shape_mask` will be 0 outside of the shape, and will be `merged_mask` inside of the shape."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "For the purpose of this demonstration, we will only select the province of Jiangsu, Zhejiang, and Shanghai."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "scrolled": false
+   },
+   "outputs": [],
+   "source": [
+    "china_shapes_subset = china_shapes[\n",
+    "    china_shapes[\"name\"].isin([\"Jiangsu\", \"Zhejiang\", \"Shanghai\"])\n",
+    "]\n",
+    "china_shapes_subset"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Converting the filtered shape dictionary to a python dictionary as the input for `extract_shapes`, where the keys for the dictionary will be the names of the new extracted shape layers."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "china_shapes_subset = (\n",
+    "    china_shapes_subset[[\"name\", \"geometry\"]].set_index(\"name\")[\"geometry\"].to_dict()\n",
+    ")\n",
+    "china_shapes_subset"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Extract the shapes from the merged_mask. "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "china.extract_shapes(china_shapes_subset)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The resulting mask object contains the dictionary `shape_mask` with the extracted values:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "china"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Saving and Loading Masks"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "china.save_mask()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "With the mask saved, the user can now load the layers or shapes with `xarray` instead if preferred."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "shape_xr_lst = china.load_shape_xr()\n",
+    "shape_xr_lst[\"Zhejiang\"].plot()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Optional: closing all the files when saving the mask. This can avoid possible write permission error."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "china.save_mask(close_files=True)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Loading a previously saved mask."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "china_2 = geodata.mask.load_mask(\"china\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "china_2"
+   ]
+  }
+ ],
+ "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.10.11"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}
\ No newline at end of file
diff --git a/docs/jupyter_execute/mask/xarray_mask_tutorial.ipynb b/docs/jupyter_execute/mask/xarray_mask_tutorial.ipynb
new file mode 100644
index 00000000..17222420
--- /dev/null
+++ b/docs/jupyter_execute/mask/xarray_mask_tutorial.ipynb
@@ -0,0 +1,331 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "9572025b",
+   "metadata": {},
+   "source": [
+    "# Tutorial: Applying Saved Masks with `XarrayMask`\n",
+    "\n",
+    "This notebook shows how to apply a **saved** geographic mask to model or analysis\n",
+    "output represented as an `xarray.Dataset` or `xarray.DataArray` — without using\n",
+    "`Cutout.add_mask` or `Cutout.mask`.\n",
+    "\n",
+    "For contributor notes on the xarray masking design, see\n",
+    "[development/xarray_mask_workflow](../development/xarray_mask_workflow.rst).\n",
+    "To build masks from rasters and shapefiles, see\n",
+    "[mask creation workflow](mask_creation_workflow.ipynb)."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "d9005d5b",
+   "metadata": {},
+   "source": [
+    "## Overview\n",
+    "\n",
+    "| Step | API | Module |\n",
+    "|------|-----|--------|\n",
+    "| Create and save a mask | `geodata.Mask` | `src/geodata/mask.py` |\n",
+    "| Run a model (wind, pvlib, …) | model `estimate()` | `src/geodata/model/` |\n",
+    "| Align mask to your grid, attach or apply | `geodata.XarrayMask` | `src/geodata/mask/xarray_mask.py` |\n",
+    "\n",
+    "**`XarrayMask` does not replace mask creation.** It loads a saved mask and applies it\n",
+    "to xarray data on your target grid."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "1bfcabc6",
+   "metadata": {},
+   "source": [
+    "## Setup\n",
+    "\n",
+    "This tutorial runs **offline** using a small synthetic grid and a temporary mask\n",
+    "directory. The same API calls work for production masks saved under `GEODATA_ROOT`."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "66976b87",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import tempfile\n",
+    "from pathlib import Path\n",
+    "\n",
+    "import numpy as np\n",
+    "import rasterio as ras\n",
+    "import shapely.geometry\n",
+    "import xarray as xr\n",
+    "from rasterio.transform import from_bounds\n",
+    "\n",
+    "from geodata import Mask, XarrayMask"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "cacb7d20",
+   "metadata": {},
+   "source": [
+    "## Step 1: Stand in for model output\n",
+    "\n",
+    "Your analysis dataset can use `x`/`y` or `lat`/`lon`. `XarrayMask` normalizes\n",
+    "coordinates via `ds_reformat_index` before alignment.\n",
+    "\n",
+    "Below we use a small `(time, y, x)` dataset as if it came from a wind or PV model."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "04086b5d",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "y = np.array([30.75, 30.5, 30.25, 30.0])\n",
+    "x = np.array([100.0, 100.25, 100.5, 100.75])\n",
+    "time = np.array([\"2016-01-01T00:00:00\", \"2016-01-01T01:00:00\"], dtype=\"datetime64[ns]\")\n",
+    "\n",
+    "values = np.arange(len(time) * len(y) * len(x), dtype=np.float32).reshape(\n",
+    "    len(time), len(y), len(x)\n",
+    ")\n",
+    "model_ds = xr.Dataset(\n",
+    "    {\"signal\": ((\"time\", \"y\", \"x\"), values)},\n",
+    "    coords={\"time\": time, \"y\": y, \"x\": x},\n",
+    ")\n",
+    "model_ds"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "19a93781",
+   "metadata": {},
+   "source": [
+    "## Step 2: Create and save a mask (offline example)\n",
+    "\n",
+    "In practice you build masks with `Mask.add_layer`, `filter_layer`, `merge_layer`,\n",
+    "and `save_mask()` — see [mask creation workflow](mask_creation_workflow.ipynb).\n",
+    "\n",
+    "Mask rasters are often stored at **higher resolution** than model output.\n",
+    "`XarrayMask` coarsens them onto `grid` automatically.\n",
+    "\n",
+    "The helper below mirrors `tests/pr/mask/test_xarray_mask.py`."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "1ff7b32a",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "mask_dir = Path(tempfile.mkdtemp(prefix=\"geodata_xmask_tutorial_\"))\n",
+    "mask_name = \"tutorial_mask\"\n",
+    "\n",
+    "lon_step = float(np.abs(x[1] - x[0]))\n",
+    "lat_step = float(np.abs(y[1] - y[0]))\n",
+    "west = float(x.min() - lon_step / 2)\n",
+    "east = float(x.max() + lon_step / 2)\n",
+    "south = float(y.min() - lat_step / 2)\n",
+    "north = float(y.max() + lat_step / 2)\n",
+    "\n",
+    "nlon_hi = len(x) * 2\n",
+    "nlat_hi = len(y) * 2\n",
+    "transform = from_bounds(west, south, east, north, nlon_hi, nlat_hi)\n",
+    "\n",
+    "arr = np.zeros((nlat_hi, nlon_hi), dtype=np.uint8)\n",
+    "arr[nlat_hi // 4 : 3 * nlat_hi // 4, nlon_hi // 4 : 3 * nlat_hi // 4] = 1\n",
+    "\n",
+    "layer_path = mask_dir / \"source.tif\"\n",
+    "with ras.open(\n",
+    "    str(layer_path),\n",
+    "    \"w\",\n",
+    "    driver=\"GTiff\",\n",
+    "    height=arr.shape[0],\n",
+    "    width=arr.shape[1],\n",
+    "    count=1,\n",
+    "    dtype=arr.dtype,\n",
+    "    compress=\"lzw\",\n",
+    "    crs=\"+proj=latlong\",\n",
+    "    transform=transform,\n",
+    ") as dst:\n",
+    "    dst.write(arr, 1)\n",
+    "\n",
+    "mask = Mask(name=mask_name, mask_dir=str(mask_dir))\n",
+    "mask.add_layer(str(layer_path), layer_name=\"source\")\n",
+    "mask.merge_layer(show_raster=False)\n",
+    "\n",
+    "region = shapely.geometry.box(west, south, (west + east) / 2, (south + north) / 2)\n",
+    "mask.extract_shapes({\"region_a\": region}, show_raster=False)\n",
+    "mask.save_mask()\n",
+    "\n",
+    "print(f\"Saved mask '{mask_name}' under {mask_dir}\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "55774df8",
+   "metadata": {},
+   "source": [
+    "## Step 3: Load and align — `XarrayMask.from_name`\n",
+    "\n",
+    "Pass your model grid so the saved mask is coarsened and aligned to the same\n",
+    "`x`/`y` (or `lat`/`lon`) coordinates."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "77bd5f99",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "xmask = XarrayMask.from_name(mask_name, grid=model_ds, mask_dir=str(mask_dir))\n",
+    "xmask"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "dfda07e9",
+   "metadata": {},
+   "source": [
+    "You can also build from an in-memory `Mask` object:\n",
+    "\n",
+    "```python\n",
+    "loaded = Mask.from_name(mask_name, mask_dir=str(mask_dir))\n",
+    "xmask = XarrayMask.from_mask(loaded, grid=model_ds)\n",
+    "```"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "1614f7a9",
+   "metadata": {},
+   "source": [
+    "## Step 4: Attach — legacy-compatible output\n",
+    "\n",
+    "`attach()` returns a dict of datasets (keys: `merged_mask`, plus any shape masks).\n",
+    "Each dataset contains your original variables plus `mask` and optional `area` — the\n",
+    "same structure as `Cutout.mask()`."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "4d30517c",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "attached = xmask.attach(model_ds, include_area=True)\n",
+    "list(attached.keys())"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "299b93e6",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "merged = attached[\"merged_mask\"]\n",
+    "merged"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "b19d9eba",
+   "metadata": {},
+   "source": [
+    "## Step 5: Apply — filtered outputs\n",
+    "\n",
+    "- `mode=\"where\"` — set values outside the mask to NaN\n",
+    "- `mode=\"multiply\"` — set values outside the mask to zero"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "44b51430",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "where_out = xmask.apply(model_ds, mode=\"where\", include_area=True)[\"merged_mask\"]\n",
+    "multiply_out = xmask.apply(model_ds, mode=\"multiply\", include_area=False)[\"merged_mask\"]\n",
+    "\n",
+    "where_out[\"signal\"].isel(time=0)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "c9ce7c5f",
+   "metadata": {},
+   "source": [
+    "## Step 6: Area-weighted aggregation\n",
+    "\n",
+    "With `attach(..., include_area=True)` you can compute mask- and area-weighted\n",
+    "statistics over time — the same pattern as the legacy Cutout workflow."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "b5135761",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "ds = attached[\"merged_mask\"]\n",
+    "weighted_mean = (\n",
+    "    (ds[\"signal\"] * ds[\"mask\"] * ds[\"area\"]).sum(dim=[\"lat\", \"lon\"])\n",
+    "    / (ds[\"mask\"] * ds[\"area\"]).sum(dim=[\"lat\", \"lon\"])\n",
+    ")\n",
+    "weighted_mean"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "6905f522",
+   "metadata": {},
+   "source": [
+    "## Production usage\n",
+    "\n",
+    "When your mask is already saved under the default mask directory (`GEODATA_ROOT`):\n",
+    "\n",
+    "```python\n",
+    "xmask = XarrayMask.from_name(\"china\", grid=output_ds)  # uses geodata.config.MASK_DIR\n",
+    "masked = xmask.apply(output_ds, mode=\"where\")\n",
+    "```\n",
+    "\n",
+    "### Typical pipeline\n",
+    "\n",
+    "1. `output_ds = model.estimate(...)`\n",
+    "2. `xmask = XarrayMask.from_name(\"my_mask\", grid=output_ds, mask_dir=...)`\n",
+    "3. `xmask.attach(output_ds)` or `xmask.apply(output_ds, ...)`\n",
+    "\n",
+    "### See also\n",
+    "\n",
+    "| Topic | Page |\n",
+    "|-------|------|\n",
+    "| Create masks from GIS layers | [mask_creation_workflow](mask_creation_workflow.ipynb) |\n",
+    "| Legacy Cutout masking | [mask_on_cutout](../legacy/mask_on_cutout.ipynb) |\n",
+    "| Xarray masking design notes (contributors) | [xarray_mask_workflow](../development/xarray_mask_workflow.rst) |\n",
+    "| Migration plan (contributors) | [mask_xarray_migration_plan](../development/mask_xarray_migration_plan.md) |\n",
+    "| Automated examples | `tests/pr/mask/test_xarray_mask.py`, `tests/pr/test_wind_xarraymask_integration.py` |"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "name": "python",
+   "version": "3.11.0"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
\ No newline at end of file
diff --git a/docs/jupyter_execute/visualization/visualization.ipynb b/docs/jupyter_execute/visualization/visualization.ipynb
new file mode 100644
index 00000000..25423766
--- /dev/null
+++ b/docs/jupyter_execute/visualization/visualization.ipynb
@@ -0,0 +1,451 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Visualization Examples"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Geodata also provides the users with different methods to visualize outputs. \n",
+    "\n",
+    "To start, import the geodata package with a logger for detailed debugging."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import geodata"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "We also import the `geopandas` and `cartopy` libraries to retrieve and show geospatial [shapefiles](https://en.wikipedia.org/wiki/Shapefile) on the plot, and the `IPython` library to download generated animation as HTML file. These libaries are helpful, but not required to use geodata for visualization."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import cartopy.io.shapereader as shpreader\n",
+    "import geopandas as gpd\n",
+    "from IPython.display import HTML"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Download example datasets and create cutouts. We will get the hourly aerosol data and the hourly radiation data."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "scrolled": false
+   },
+   "outputs": [],
+   "source": [
+    "# Download aerosol hourly data\n",
+    "aerosol_hourly_data = geodata.Dataset(\n",
+    "    module=\"merra2\",\n",
+    "    years=slice(2020, 2020),\n",
+    "    months=slice(1, 12),\n",
+    "    weather_data_config=\"surface_aerosol_hourly\",\n",
+    ")\n",
+    "\n",
+    "# Download radiation hourly data\n",
+    "slv_hourly_data = geodata.Dataset(\n",
+    "    module=\"merra2\",\n",
+    "    years=slice(2011, 2011),\n",
+    "    months=slice(1, 1),\n",
+    "    weather_data_config=\"slv_radiation_hourly\",\n",
+    ")\n",
+    "\n",
+    "if aerosol_hourly_data.prepared == False:\n",
+    "    aerosol_hourly_data.get_data()\n",
+    "\n",
+    "# Download radiation hourly data only on 2011/01/01\n",
+    "if slv_hourly_data.prepared == False:\n",
+    "    slv_hourly_data.get_data(testing=True)\n",
+    "\n",
+    "# Create northern china aerosol Cutout\n",
+    "cutout_pm25 = geodata.Cutout(\n",
+    "    name=\"beijing19\",\n",
+    "    module=\"merra2\",\n",
+    "    weather_data_config=\"surface_aerosol_hourly\",\n",
+    "    xs=slice(105, 123),\n",
+    "    ys=slice(27, 43),\n",
+    "    years=slice(2020, 2020),\n",
+    "    months=slice(1, 12),\n",
+    ")\n",
+    "\n",
+    "# Create china solar Cutout\n",
+    "cutout_solar = geodata.Cutout(\n",
+    "    name=\"china-2011-slv-hourly-test\",\n",
+    "    module=\"merra2\",\n",
+    "    weather_data_config=\"slv_radiation_hourly\",\n",
+    "    xs=slice(73, 136),\n",
+    "    ys=slice(18, 54),\n",
+    "    years=slice(2011, 2011),\n",
+    "    months=slice(1, 1),\n",
+    ")\n",
+    "\n",
+    "cutout_solar.prepare()\n",
+    "cutout_pm25.prepare()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Generate PM2.5 and Solar PV Outputs."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "ds_pm25 = geodata.convert.pm25(cutout_pm25)\n",
+    "ds_solar = geodata.convert.pv(cutout_solar, panel=\"KANEKA\", orientation=\"latitude_optimal\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Time Series Visualization"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Default time series method call"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "We can use `geodata.plot.time_series` to visualize time series data from the output xarray DataArray, such as `ds_pm25` or `ds_solar`. Its minimal method call find the mean value of all grid cell for every time point in the dataset. For example, with `ds_solar`, we can visualize the spatially aggregated averages AC power over time."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "geodata.plot.time_series(ds_solar)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Spatial and temporal aggregation"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The `time_series` method can take in tuple parameters `lat_slice` and `lon_slice` to select grid cells within that range (inclusive). For example, if we want to find the aggregated value for all grid cells between latitude 35 degree and 36 degree, we set `lat_slice` to be (35, 36). The `agg_slice_method` parameter will specify the aggregation method for aggregating grid cells sliced by `lat_slice` or `lon_slice`. By default, `agg_slice_method` is set to mean aggregation. \n",
+    "\n",
+    "We use the latitude-sliced time-series visualization on the PM2.5 output below. Note that since we have hourly data for the year 2019, we will have 24 * 365 = 8760 timepoints for each hour. However, we can reduce the number of timepoints by taking in a `time_factor` parameter that tells the method how many timepoints to aggregate on. Here, we take 24 * 7 as the `time_factor` so that we will aggregate the data by week, as there are 24 * 7 hours in a week. The `agg_time_method` parameter will specify the aggregation method for time aggregation. By default, `agg_time_method` is set to mean aggregation. "
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "For example, below we visualize the weekly averages of sum of PM2.5 for region within latitude slice (35, 36)."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "geodata.plot.time_series(ds_pm25, lat_slice=(35, 36), agg_slice_method=\"sum\", time_factor=24 * 7)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "If we have `lat_slice` or `lon_slice` inputs, and want to plot the time series for every single grid cell without aggregating them, they can specify `agg_slice = False`. This will generate one line for each grid cell.\n",
+    "\n",
+    "The method also takes in user-defined title with the `title` parameter."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "geodata.plot.time_series(\n",
+    "    ds_pm25,\n",
+    "    lat_slice=(35, 36),\n",
+    "    lon_slice=(110, 111),\n",
+    "    agg_slice=False,\n",
+    "    time_factor=24 * 7,\n",
+    "    title=\"PM2.5 Time Series - lat(35-36) lon(110-111) weekly average\",\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Multiple coordinate points"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "You can also use a dictionary of name-coordinate pairs to plot different grid cells. The coordinates value of this `coord_dict` does not have to be exact, as the method can automatically find the grid cell containing the coordinate input."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "coord_d = {\"Beijing\": (30.9, 116.4), \"Shanghai\": (31.2, 121.47), \"Xi'an\": (34.2, 108.9)}\n",
+    "\n",
+    "geodata.plot.time_series(ds_pm25, coord_dict=coord_d, time_factor=24 * 7)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Heatmap Visualization"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Default Method Call"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Geodata can plot a spatial heatmap of output values. Since the output is a time-series containing more than 2 dimensions, this method will aggregate the values by mean at different timepoints for each grid cells by default. For example, to see the annual mean PM2.5 in our Cutout region, we use the following method call:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "geodata.plot.heatmap(ds_pm25)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Add shapefiles to Plot"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The `heatmap` method can also take in a `shape` parameter, which takes in a `geopandas` dataframe or series of shape objects. Let us use the province shapes from `cartopy` shape-reader and save the path as `prov_path`. This can also be the path to user-supplied shape files. "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "prov_path = shpreader.natural_earth(resolution=\"10m\", category=\"cultural\", name=\"admin_1_states_provinces\")\n",
+    "shapes = gpd.read_file(prov_path, encoding=\"utf-8\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "geodata.plot.heatmap(ds_pm25, shape=shapes)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Selecting Timepoint"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "If we do not want the temporally aggregated plot, we can specify the exact time point or its index in the dataArray. In the following method call, `t = 0` uses index to select the first time point in `ds_pm25`."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "geodata.plot.heatmap(ds_pm25, t=0, shape=shapes)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "We can also take in the exact time point from `ds_pm25` as a string. We can also change the map type from the default `colormesh` to `contour`, and customize the title text like the following:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "geodata.plot.heatmap(\n",
+    "    ds_pm25,\n",
+    "    t=\"2019-01-01T00:30:00\",\n",
+    "    map_type=\"contour\",\n",
+    "    shape=shapes,\n",
+    "    title=\"Contour plot\",\n",
+    "    title_size=20,\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Let's use the `heatmap` method on the solar PV output xarray `ds_solar`. Below we select the 7th time point for the `ds_solar` dataArray with the provincial shapes on the same plot.\n",
+    "\n",
+    "Note that the default map color of the method is `bone_r`, which is not ideal for visualizing solar PV. Therefore, we switch the `cmap` parameter to `Wistia`. You can view a complete list of matplotlib map color [here](https://matplotlib.org/stable/gallery/color/colormap_reference.html).\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "geodata.plot.heatmap(\n",
+    "    ds_solar,\n",
+    "    t=6,\n",
+    "    shape=shapes,\n",
+    "    shape_width=0.25,\n",
+    "    shape_color=\"navy\",\n",
+    "    map_type=\"contour\",\n",
+    "    cmap=\"Wistia\",\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Animation"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The drawback of plotting a static heatmap with `heatmap` is that we cannot see the changes over time like the `time_series` plots. However, the `heatmap_animation` method can create an animation of heatmap with time as another dimension in the plot.\n",
+    "\n",
+    "The parameters of the heatmap_animation is very similar to the ones for `heatmap`. You can use `time_factor` to find aggregated mean or sum. Here, we create the animation with averages for every two hours in the day."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "geodata.plot.heatmap_animation(\n",
+    "    ds_solar,\n",
+    "    cmap=\"Wistia\",\n",
+    "    time_factor=2,\n",
+    "    shape=shapes,\n",
+    "    shape_width=0.25,\n",
+    "    shape_color=\"navy\",\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The users can save the animation to a file, which requires the `HTML` method from the `IPython` package we imported earlier. It also requires the users to use the Jupyter Notebook in a browser, and have already generated the heatmap animation in the notebook, because `geodata.plot.save_animation` will extract the javascript content string from the animation in the Jupyter Notebook, and use HTML() method to enable the browser to download the file."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Save the animation above as a file named `solar_pv_2011_01_01_animation.html`."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "HTML(geodata.plot.save_animation(\"solar_pv_2011_01_01_animation.html\"))"
+   ]
+  }
+ ],
+ "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.10.11"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}
\ No newline at end of file
diff --git a/docs/requirements.txt b/docs/requirements.txt
new file mode 100644
index 00000000..4a487e8a
--- /dev/null
+++ b/docs/requirements.txt
@@ -0,0 +1,8 @@
+# Sphinx stack (matches pyproject [project.optional-dependencies] docs)
+sphinx>=8.0.0
+myst-nb>=1.1.2
+sphinx-book-theme>=1.1.3
+sphinx-autoapi==3.3.2
+
+# Install geodata itself without pulling the full runtime dependency tree
+# (see .readthedocs.yaml: pip install --no-deps .)
diff --git a/docs/source/conf.py b/docs/source/conf.py
index ea9d69e7..8262e956 100644
--- a/docs/source/conf.py
+++ b/docs/source/conf.py
@@ -6,12 +6,21 @@
 # -- Project information -----------------------------------------------------
 # https://www.sphinx-doc.org/en/master/usage/configuration.html#project-information
 
-from geodata import __version__
+import re
+from pathlib import Path
+
+# Read version without importing geodata (RTD installs with --no-deps).
+_version_file = Path(__file__).resolve().parents[2] / "src" / "geodata" / "_version.py"
+_release = re.search(
+    r'^__version__\s*=\s*["\']([^"\']+)["\']', _version_file.read_text(), re.M
+)
+if _release is None:
+    raise RuntimeError(f"Could not parse __version__ from {_version_file}")
+release = _release.group(1)
 
 project = "Geodata"
 copyright = "2025, Geodata Contributors"
 author = "Geodata Contributors"
-release = __version__
 
 # -- General configuration ---------------------------------------------------
 # https://www.sphinx-doc.org/en/master/usage/configuration.html#general-configuration
diff --git a/docs/source/datasets/era5.rst b/docs/source/datasets/era5.rst
index 2af0aee6..116cc9aa 100644
--- a/docs/source/datasets/era5.rst
+++ b/docs/source/datasets/era5.rst
@@ -1,81 +1,89 @@
 ERA5 Specific Instructions
 ==========================
 
-This page explains how you can set up access to ERA5 data from the `Copernicus Data Store <https://cds.climate.copernicus.eu/cdsapp#!/dataset/reanalysis-era5-single-levels?tab=overview>`_.
+This page covers **CDS account and API credential setup** for ERA5. Once credentials
+are in place, use the dataset classes — do not call ``cdsapi`` by hand for routine
+downloads.
 
-Creating a CDS account
-----------------------
-
-To download ERA5 data from the CDS, you'll need to create a free `CDS account here <https://cds.climate.copernicus.eu/>`_.
-
-Download data through CDS API
+Recommended download method
 -----------------------------
 
-Once your account has been created, set up access to the API by following these steps:
+The **recommended way** to fetch ERA5 data in Geodata is:
 
-1. Log into your CDS account and visit your `profile page <https://cds.climate.copernicus.eu/profile>`_.
-2. Install the API key. There will be a section called **Personal Access Token**.
-   Copy these two lines into a file called ``.cdsapirc`` in your user root folder.
+1. Complete the CDS setup below (one-time).
+2. Follow :ref:`downloading-era5-data` in :doc:`overview` — ``load_dataset``,
+   instantiate with ``years`` / ``months`` / optional ``bounds``, then ``download()``.
 
-- **macOS/Linux**: Open a terminal and run:
+Geodata's ERA5 classes (for example ``ERA5Wind3DHourlyDataset``) create a
+``cdsapi.Client`` internally and submit the correct product requests for each
+registered ``weather_config``.
 
-  .. code-block:: bash
+Creating a CDS account
+----------------------
 
-      touch ~/.cdsapirc
+To download ERA5 data from the CDS, create a free `CDS account here <https://cds.climate.copernicus.eu/>`_.
 
-  Then add the lines using:
+Configure CDS API credentials
+-----------------------------
 
-  .. code-block:: bash
+Once your account exists, install local API access:
 
-      echo [line 1 of the code] >> ~/.cdsapirc
-      echo [line 2 of the code] >> ~/.cdsapirc
+1. Log into your CDS account and visit your `profile page <https://cds.climate.copernicus.eu/profile>`_.
+2. Under **Personal Access Token**, copy the two lines for your ``.cdsapirc`` file
+   (URL and key).
 
+**macOS/Linux** — create ``~/.cdsapirc``:
 
-  - **Windows**: The process is slightly more complicated. Please refer to the in-depth guide at the Copernicus Knowledge Base `here <https://confluence.ecmwf.int/display/CKB/How+to+install+and+use+CDS+API+on+Windows>`_.
+.. code-block:: bash
 
-3. Install the CDS API client by opening a terminal/shell and running
+   touch ~/.cdsapirc
+   # Paste the two lines from your CDS profile into ~/.cdsapirc
 
-.. code-block:: bash
+**Windows** — see the Copernicus guide on
+`installing the CDS API on Windows <https://confluence.ecmwf.int/display/CKB/How+to+install+and+use+CDS+API+on+Windows>`_.
 
-     pip install ".[download]"
+Ensure ``cdsapi`` is available (it is a dependency of Geodata when you install the
+package). Then proceed to :ref:`downloading-era5-data` in :doc:`overview`.
 
-(Assuming you are in Geodata's *root directory*.)
+Verify CDS API access (optional)
+--------------------------------
 
-1. Once you've installed the API key and the API client, confirm access by running an
-   example in a Python script or a Jupyter notebook:
+You can confirm credentials with a minimal ``cdsapi`` script. This is **optional** —
+Geodata dataset downloads use the same client and credentials.
 
 .. code-block:: python
 
-     import cdsapi
-
-     c = cdsapi.Client()
-
-     c.retrieve(
-          "reanalysis-era5-single-levels",
-          {
-               "product_type": "reanalysis",
-               "format": "netcdf",
-               "variable": [
-                    "2m_dewpoint_temperature",
-                    "2m_temperature",
-               ],
-               "year": "2011",
-               "month": [
-                    "01",
-               ],
-               "day": ["01", "02", "03"],
-               "time": [
-                    "00:00",
-                    "12:00",
-               ],
-          },
-          "download.nc",
-     )
-
-The above example downloads 2m temperature and 2m dewpoint temperature with data points
-at 00:00 and 12:00 for each day, from January 1-3, 2011, in NetCDF format.
-
-If this works, you have successfully set up access to the ERA5 data through the CDS API.
-Please subsequently refer to the `general documentation on datasets <../overview.rst>`_
-for more information on how to download ERA5-based datasets using the ``geodata``
-package.
+   import cdsapi
+
+   c = cdsapi.Client()
+
+   c.retrieve(
+       "reanalysis-era5-single-levels",
+       {
+           "product_type": "reanalysis",
+           "format": "netcdf",
+           "variable": [
+               "2m_dewpoint_temperature",
+               "2m_temperature",
+           ],
+           "year": "2011",
+           "month": ["01"],
+           "day": ["01", "02", "03"],
+           "time": ["00:00", "12:00"],
+       },
+       "download.nc",
+   )
+
+This example fetches 2 m temperature and dewpoint at 00:00 and 12:00 UTC for
+2011-01-01 through 2011-01-03. If it succeeds, your ``.cdsapirc`` is valid.
+
+For production workflows, prefer :ref:`downloading-era5-data` in :doc:`overview` so
+Geodata requests the correct ERA5 products, paths, and post-processing for
+``wind_3d_hourly``, ``wind_solar_hourly``, and other registered configs.
+
+What's next
+-----------
+
+- :ref:`downloading-era5-data` in :doc:`overview` — **recommended** download workflow
+- :doc:`../development/offline-era5-fixture-datasets` — offline ``*_test`` configs for CI
+- :doc:`../modeling/wind/index` or :doc:`../modeling/pvlib/index` — run models on downloaded data
diff --git a/docs/source/datasets/overview.rst b/docs/source/datasets/overview.rst
index 57430059..d6622e48 100644
--- a/docs/source/datasets/overview.rst
+++ b/docs/source/datasets/overview.rst
@@ -9,85 +9,127 @@ data formats, handling metadata, and performing common geospatial operations.
 Key Features
 ------------
 
-- Supports the download and management of datasets from various sources, such as
-  `ERA5 <https://www.ecmwf.int/en/forecasts/dataset/ecmwf-reanalysis-v5>`_ and
-  `MERRA2 <https://gmao.gsfc.nasa.gov/reanalysis/MERRA-2/>`_.
+- Supports the download and management of **ERA5** datasets via ``load_dataset`` (see :doc:`era5` for CDS account setup).
 
-- Provides a consistent API for accessing geospatial data, regardless of the underlying
-  data source.
+- **MERRA2** remains in the codebase but is documented under :doc:`/legacy/index` (legacy ``Dataset`` / ``Cutout`` path, not part of the current tested workflow).
 
 Typical Usage
 -------------
 
-In the following example, we will demonstrate how to download a dataset containing wind
-and solar data from ECMWF's ERA5 dataset.
+Registered datasets are loaded by name, instantiated with a time range (and optional
+geographic bounds), then downloaded with ``download()`` if the files are not already
+on disk. The sections below use **ERA5** as the primary example; the same pattern
+applies to other configs returned by ``list_datasets()``.
+
+.. _downloading-era5-data:
+
+Downloading ERA5 data
+---------------------
+
+Geodata's ERA5 dataset classes wrap the `Copernicus CDS API <https://cds.climate.copernicus.eu/>`_.
+You configure credentials once (see :doc:`era5`), then download through Python — you do
+**not** need to call ``cdsapi`` directly for normal use.
+
+Files are stored under ``GEODATA_ROOT / era5 / <weather_config> / …`` (see
+:doc:`../quick_start/packagesetup` for ``GEODATA_ROOT``).
+
+**Example — 3D wind (for :doc:`../modeling/wind/index`):**
 
 .. code-block:: python
 
     from geodata.datasets import load_dataset
 
-    dataset_cls = load_dataset("wind_solar_hourly")
+    ds_cls = load_dataset("wind_3d_hourly")
+    ds = ds_cls(
+        years=slice(2016, 2016),
+        months=slice(1, 1),
+        bounds=[-10, 35, 10, 45],  # lon_min, lat_min, lon_max, lat_max
+    )
+
+    print(ds.downloaded)  # False until files exist locally
+
+    if not ds.downloaded:
+        ds.download()
+
+    print(ds.downloaded)  # True when the catalog is complete
+
+**Example — 2D wind and solar hourly (for :doc:`../modeling/pvlib/index`):**
+
+.. code-block:: python
+
+    from geodata.datasets import load_dataset
 
-    years = slice(2010, 2020)
-    months = slice(1, 13)
-    dataset = dataset_cls(years=years, months=months)
+    ds_cls = load_dataset("wind_solar_hourly")
+    ds = ds_cls(years=slice(2016, 2016), months=slice(1, 1))
 
-Here, we first create a dataset class using the `load_dataset` function, specifying the
-name of the dataset we want to load. We then instantiate the dataset class with the
-desired time range (years and months). Then, we can create a dataset instance with
-that class, which will handle the downloading and processing of the data.
+    if not ds.downloaded:
+        ds.download()
+
+``bounds`` is optional; omit it to use the full spatial extent allowed by the dataset
+class. With ``testing=True``, only a **small subset** of the catalog is requested (useful
+for trying a download before committing to a full month):
+
+.. code-block:: python
+
+    ds = ds_cls(
+        years=slice(2016, 2016),
+        months=slice(1, 1),
+        bounds=[50, 0, 48, 3],
+        testing=True,
+    )
+    ds.download()
+
+After ``downloaded`` is ``True``, pass ``ds`` to a model (for example
+``WindInterpolationModel(ds)`` or ``Pvlib(ds)``).
+
+Offline / CI without CDS
+~~~~~~~~~~~~~~~~~~~~~~~~
+
+For tests and local development without calling the CDS, use the committed fixture
+configs ``wind_3d_hourly_test`` and ``wind_solar_hourly_test`` — same API, no
+``download()`` required when fixture files are present. See
+:doc:`../development/offline-era5-fixture-datasets`.
 
 Dataset Classes
 -----------------
-The `geodata.datasets` module includes several dataset classes, each tailored for
+
+The ``geodata.datasets`` module includes several dataset classes, each tailored for
 specific datasets. These classes encapsulate the logic for downloading, processing, and
-accessing the data. Some of the available dataset classes
-(listed by `weather_data_config`) include:
+accessing the data. Some of the available ERA5 configs
+(listed by ``weather_config``) include:
 
-- `wind_solar_hourly`: A dataset containing hourly wind and solar data from ECMWF's
-  ERA5. It is important to note that the wind data are only recorded at
-  10 and 100 meters above ground level. Hence, this dataset is also referred to as
-  2D wind and solar dataset.
+- ``wind_solar_hourly``: Hourly wind (10 m and 100 m) and solar radiation from ERA5
+  single levels. Also referred to as the 2D wind and solar dataset.
 
-- `wind_3d_hourly`: A dataset containing hourly wind data from ECMWF's ERA5 at
-  multiple vertical levels, providing a more comprehensive view of the wind profile.
-  It can be used for wind speed estimation using and interpolation model built into
-  the geodata library.
+- ``wind_3d_hourly``: Hourly wind on ERA5 model levels (131–137), stored as **daily**
+  NetCDF files. Used by the wind interpolation model for hub-height wind speed.
 
-You can use the `list_datasets` function to see all available datasets in the
-`geodata.datasets` module. This function returns a list of dataset names that can be
-loaded using the `load_dataset` function. For example:
+You can use ``list_datasets()`` to see all registered names:
 
 .. code-block:: python
 
     from geodata.datasets import list_datasets
 
-    available_datasets = list_datasets()
-    print(available_datasets)  # Outputs a list of available dataset names.
-
-Check Preparedness of Datasets
-------------------------------------------------
-To check if a dataset is prepared and ready for use, you can use the `downloaded`
-property of the dataset instance. This property returns a boolean indicating whether the
-dataset is fully prepared. If the dataset is not prepared, you can call the `prepare`
-method to download and process the data. For example:
+    print(list_datasets())
 
-.. code-block:: python
+Check whether data is on disk
+------------------------------
 
-    print(dataset.downloaded)  # Check if the dataset is downloaded. Outputs False here.
+The ``downloaded`` property is ``True`` when every file in the dataset **catalog** exists
+under ``storage_root``. If any file is missing, call ``download()`` (or ``download(force=True)``
+to re-fetch):
 
-    if not dataset.downloaded:
-        dataset.download()
+.. code-block:: python
 
-    print(dataset.downloaded)  # Outputs True after downloading.
+    if not ds.downloaded:
+        ds.download()
 
 Dataset's Interoperability with Cutout
 ------------------------------------------------
 
-At the moment, the dataset classes are not interoperable with the `Cutout` class.
-In the future, we plan to consolidate the functionalities of the `Cutout` class into the
-dataset classes and the modeling module (see :doc:`here<../modeling/wind/index>`).
+At the moment, the dataset classes are not interoperable with the ``Cutout`` class.
+In the future, we plan to consolidate the functionalities of the ``Cutout`` class into the
+dataset classes and the modeling module (see :doc:`../modeling/wind/index`).
 
-For now, after downloading a dataset, a good point to move forward would be to use the
-:doc:`modeling module <../modeling/index>` to create a model that can do certain types
-of modeling with the dataset, such as wind speed estimation or solar PV generation.
+For now, after downloading a dataset, pass it to a modeling class — see
+:doc:`../modeling/wind/index` or :doc:`../modeling/pvlib/index`.
diff --git a/docs/source/datasets/weather_data_config.md b/docs/source/datasets/weather_data_config.md
index 098c64a4..8dcf9748 100644
--- a/docs/source/datasets/weather_data_config.md
+++ b/docs/source/datasets/weather_data_config.md
@@ -4,8 +4,10 @@ In Geodata, every downloadable dataset are associated with a unique `(module, we
 In this tuple, the `module` typicallly refers to the source of dataset, while the `weather_data_config` is 
 a dictionary that contains the information needed to download the specific form of the dataset. 
 
-As Geodata currently supports `ERA5` and `MERRA2` modules, you can find all relevant weather data configuration
-in each module's introduction pages here ([ERA5](era5/index.md), [MERRA2](merra2/index.md)). To find each config's actual definition, you can go to `src/geodata/datasets`. Within it, all available weather data configurations are located at the bottom of the file. 
+Geodata supports ERA5 through the modern `load_dataset` registry and MERRA2 through the legacy
+`Dataset(module="merra2", ...)` API. Introduction pages: [ERA5](era5.rst),
+[MERRA2 (legacy)](../legacy/merra2/index.md). To find each config's actual definition, go to
+`src/geodata/datasets` — weather data configurations are defined at the bottom of each module file. 
 
 In this tutorial, we will discuss the structure of each `weather_data_config` in more details.
 
diff --git a/docs/source/development/documentation-organization-plan.md b/docs/source/development/documentation-organization-plan.md
new file mode 100644
index 00000000..8d8ca91c
--- /dev/null
+++ b/docs/source/development/documentation-organization-plan.md
@@ -0,0 +1,396 @@
+# Geodata documentation organization plan
+
+This document defines how Geodata documentation is structured, how it maps to
+`src/geodata`, and how we keep prose, notebooks, and API reference aligned as
+the library evolves. It is intended for contributors working on the
+**documentation branch** and for anyone opening a PR that changes user-facing
+behavior.
+
+**Status:** living plan (update this file when conventions change).
+
+---
+
+## 1. Goals
+
+1. **One clear user journey** — readers should know whether to use the modern
+   (`load_dataset` → model → optional masking) or legacy (`Dataset` → `Cutout` →
+   `convert`) workflow without reading the entire site.
+2. **Docs follow code** — every public API change in `src/` has a defined doc
+   touchpoint (prose, notebook, or autoapi docstring).
+3. **No orphan pages** — every `.md`, `.rst`, and `.ipynb` under
+   `docs/source/` appears in a `toctree` or is explicitly marked as internal
+   (see [Section 5](#5-file-types-and-conventions)).
+4. **Reproducible examples** — tutorials should run offline where possible
+   (ERA5 `*_test` fixtures) so CI and local builds do not depend on CDS
+   credentials.
+5. **Separation of concerns** — migration plans and design notes stay in
+   `development/` or clearly labeled plan pages; user-facing tutorials stay
+   task-focused.
+
+---
+
+## 2. Current state (baseline)
+
+### 2.1 Two parallel workflows
+
+Geodata currently exposes two stacks. Both are valid; documentation must label
+them explicitly.
+
+| Aspect | Modern workflow | Legacy workflow |
+|--------|-----------------|-----------------|
+| Data access | `geodata.datasets.load_dataset(...)` | `geodata.Dataset(module=..., weather_data_config=...)` |
+| Subsetting | `BaseDataset` bounds / model `xs`/`ys` | `geodata.Cutout` + `prepare()` |
+| Transform | `geodata.model.wind`, `geodata.model.pvlib` | `geodata.convert.*` on Cutouts |
+| Masking (apply) | `geodata.XarrayMask` | `cutout.add_mask()` + `cutout.mask()` |
+| Masking (create) | `geodata.Mask` (same for both) | `geodata.Mask` (same for both) |
+| Primary docs | `datasets/`, `modeling/` | `intro.rst`, mask Cutout notebooks |
+
+**Canonical path for new features:** modern workflow. Legacy paths remain
+documented until explicitly deprecated.
+
+### 2.2 Documentation build stack
+
+| Piece | Location | Role |
+|-------|----------|------|
+| Sphinx config | `docs/source/conf.py` | MyST, notebooks, autoapi |
+| Site root | `docs/source/index.rst` | Top-level toctrees |
+| Landing narrative | `docs/source/intro.rst` | Overview (still legacy-heavy) |
+| API reference | autoapi → `src/geodata` | Generated from docstrings |
+| Notebooks | `myst_nb`, `nb_execution_mode = "off"` | Committed outputs; not executed on build |
+
+### 2.3 Known gaps (as of this plan)
+
+| Gap | Impact | Priority | Status |
+|-----|--------|----------|--------|
+| `intro.rst` teaches legacy Cutout/convert as the main story | New users miss models + `XarrayMask` | P0 | **Done** — modern intro on homepage; legacy moved to `legacy/workflow.rst` |
+| Modeling pages missing recent API options (`compact_output`, flexible `xs`/`ys`) | Docs diverge from `src` | P0 | **Done** — see modeling/wind/index and modeling/pvlib/index |
+| Wind capacity-factor internals not in wind toctree | Deep-dive exists only in source/comments | P1 | **Done (Option A)** — “Understanding the output” in interpolation/extrapolation Step 5 |
+| `xarray_mask_tutorial.ipynb` referenced by `xarray_mask_workflow.rst` but may be missing from tree | Broken `:doc:` link | P0 | **Done** — notebook added under `mask/` |
+| Mask section mixes user tutorials with `mask_xarray_migration_plan.md` | Hard to tell “how-to” vs “plan” | P1 | **Done** — plans moved to `development/` |
+| `development/offline-era5-fixture-datasets.md` not linked from modeling tutorials | Readers assume CDS required | P1 | **Done** — overview + era5 + intro link fixtures |
+| Example scripts in `docs/source/mask/*.py` not classified | Unclear if maintained or one-off | P2 | Open |
+| README points to placeholder doc URL | External discoverability | P2 | Open |
+
+---
+
+## 3. Target information architecture
+
+Organize the site by **user task**, not by file type. Recommended sidebar
+structure (matches `index.rst` with clearer intent):
+
+```
+Geodata docs
+├── Getting started
+│   ├── Package setup
+│   ├── Supported I/O formats
+│   └── Workflow chooser (NEW — short page: modern vs legacy)
+├── Datasets
+│   ├── Overview (load_dataset, list_datasets)
+│   ├── ERA5 (CDS setup + configs)
+│   ├── MERRA2
+│   └── Weather data config reference
+├── Modeling
+│   ├── Wind (index + interpolation + extrapolation; CF notes in Step 5)
+│   └── PVLib (index + future subpages)
+├── Masking
+│   ├── Create masks (mask_creation_workflow.ipynb)
+│   ├── Apply with XarrayMask (xarray_mask_tutorial.ipynb)
+│   ├── Troubleshoot (mask_troubleshoot.md)
+│   └── Cutout apply → legacy/mask_on_cutout
+├── Visualization
+├── Development (contributors)
+│   ├── Documentation organization (this file)
+│   ├── Offline ERA5 fixtures
+│   ├── mask_xarray_migration_plan
+│   └── xarray_mask_workflow (implementation notes)
+└── API reference (autoapi)
+```
+
+### 3.1 Page roles (Diátaxis)
+
+Use four doc types consistently:
+
+| Type | Purpose | Examples |
+|------|---------|----------|
+| **Tutorial** | Learning-oriented, step-by-step | Notebooks, `modeling/wind/interpolation.rst` |
+| **How-to guide** | Goal-oriented recipe | `xarray_mask_tutorial.ipynb`, ERA5 CDS setup |
+| **Reference** | Accurate, complete | autoapi, `weather_data_config.md`, turbine YAML lists |
+| **Explanation** | Concepts and design | migration plans, wind CF summary in interpolation Step 5 |
+
+Label migration/plan documents at the top:
+
+```markdown
+> **Audience:** contributors and maintainers. For usage, see [XarrayMask tutorial](../mask/xarray_mask_tutorial.ipynb).
+```
+
+---
+
+## 4. Source code ↔ documentation map
+
+Maintain this table when adding modules. **Primary doc** is the page that must
+be updated first when behavior changes.
+
+| `src/geodata` area | Primary doc | Secondary / API |
+|--------------------|-------------|-----------------|
+| `datasets/_base.py`, `datasets/era5/*` | `datasets/overview.rst`, `datasets/era5.rst` | autoapi |
+| `datasets/merra2/*` (legacy) | `legacy/merra2/*` | autoapi |
+| `datasets/era5/fixture.py` (`*_test`) | `development/offline-era5-fixture-datasets.md` | modeling tutorials (offline note) |
+| `model/wind/*` | `modeling/wind/index.rst`, `interpolation.rst`, `extrapolation.rst` | autoapi |
+| `model/pvlib/_base.py` | `modeling/pvlib/index.rst` | autoapi |
+| `model/_base.py` (slice sel, I/O) | modeling pages (bounding box sections) | autoapi |
+| `mask.py` (legacy Mask) | `mask/mask_creation_workflow.ipynb` | autoapi |
+| `mask/xarray_mask.py`, `mask/spatial.py` | `mask/xarray_mask_tutorial.ipynb` | `development/xarray_mask_workflow.rst` (contributors) |
+| `cutout.py`, `convert.py`, `preparation.py` | `legacy/workflow.rst` | autoapi |
+| Cutout-based masking | `legacy/mask_on_cutout.ipynb` | — |
+| `plot.py` | `visualization/visualization.ipynb` | autoapi |
+| `resource.py`, `resources/*` | modeling pages (turbine/panel names) | — |
+| `config.py` | `quick_start/packagesetup.md` | — |
+
+### 4.1 Public exports (`__init__.py`)
+
+When adding or removing symbols from `geodata.__all__`:
+
+1. Update docstrings (autoapi).
+2. Update `intro.rst` or the relevant tutorial if the symbol is part of a
+   documented workflow.
+3. Add a line to the [changelog section](#72-changelog-expectations) of the PR.
+
+---
+
+## 5. File types and conventions
+
+### 5.1 Where files live
+
+| Path | Use for |
+|------|---------|
+| `docs/source/quick_start/` | Install, env vars, I/O formats |
+| `docs/source/datasets/` | Download, configs, dataset-specific outputs |
+| `docs/source/modeling/<domain>/` | Model tutorials and domain index |
+| `docs/source/mask/` | Mask tutorials, workflows, troubleshooting |
+| `docs/source/visualization/` | Plotting notebooks |
+| `docs/source/development/` | Contributor docs, fixtures, **this plan**, internal design |
+| `docs/source/_static/` | Images referenced from rst/md |
+
+### 5.2 Format choice
+
+| Format | When to use |
+|--------|-------------|
+| `.rst` | Sphinx-native pages with toctrees (section indexes) |
+| `.md` (MyST) | Prose guides, plans, troubleshooting |
+| `.ipynb` | Executable narratives with plots; keep outputs committed |
+
+### 5.3 Naming
+
+- User-facing: `snake_case` or `kebab-case` descriptive names
+  (`xarray_mask_tutorial.ipynb`, `mask_troubleshoot.md`).
+- Plans: suffix or folder under `development/` (`*_plan.md`, `*_known_issues.md`).
+- Example scripts: `docs/source/<topic>/examples/` (proposed) — not mixed with
+  built pages unless listed in toctree.
+
+### 5.4 Internal vs published pages
+
+Pages under `development/` and `mask/*_plan.md` are **contributor-facing**.
+They stay in the toctree under **Development** or with an audience banner so
+users are not sent to migration checklists by mistake.
+
+Optional future convention: prefix internal-only files with `_` and exclude in
+`conf.py` `exclude_patterns` — not required if audience banners are used.
+
+---
+
+## 6. Keeping documentation up to date
+
+### 6.1 PR checklist (code changes)
+
+Every PR that changes `src/geodata` should answer:
+
+- [ ] Does this change **public API** or default behavior?
+- [ ] Which **primary doc** row in [Section 4](#4-source-code--documentation-map) applies?
+- [ ] Are **docstrings** updated for autoapi?
+- [ ] Is there a **minimal code snippet** in prose docs or a test that can be copied?
+- [ ] Do **notebooks** need re-run outputs (if affected)?
+- [ ] Does `intro.rst` need a **workflow label** (modern vs legacy) if touched?
+
+If the answer to the first question is yes and no doc file is updated, the PR
+should either include doc updates or link a follow-up issue.
+
+### 6.2 Documentation-only PRs (this branch)
+
+Recommended batching for the documentation branch:
+
+| Phase | Work | Outcome |
+|-------|------|---------|
+| **A — Structure** | Workflow chooser; relabel legacy in `intro.rst`; wire orphan pages into toctrees | Clear navigation |
+| **B — Sync with recent `src`** | `XarrayMask`, `compact_output`, slice/bounds notes, fixture offline path | Factual parity with code |
+| **C — Depth** | Wind CF explanation, mask troubleshooting, merge-layer known issues | Explanation layer |
+| **D — Hygiene** | Move example `.py` to `examples/`; README doc URL; trim stale “planned” notes | Lower maintenance cost |
+
+### 6.3 When to update which layer
+
+| Change in `src` | Update prose/notebook | Update docstrings only |
+|-----------------|----------------------|-------------------------|
+| New public class or method | Yes | Yes |
+| New optional parameter with non-obvious default | Yes (one example) | Yes |
+| Internal refactor, same API | No | Only if signatures changed |
+| Bug fix affecting coordinates, units, or outputs | Yes (note in troubleshooting or tutorial) | Yes |
+| New `*_test` fixture config | `development/offline-era5-fixture-datasets.md` | — |
+| Deprecation | Yes + migration plan | Yes |
+
+### 6.4 Single source of truth
+
+| Content | Source of truth | Docs should… |
+|---------|-----------------|--------------|
+| Function signatures | `src/` + autoapi | Not duplicate parameter lists |
+| End-to-end workflows | Notebooks + rst tutorials | Link to tests under `tests/pr/` |
+| Dataset registry names | `list_datasets()` / `datasets/registry` | Regenerate or manually sync lists in `overview.rst` when configs added |
+| Turbine/panel names | `resources/windturbine/`, `resources/solarpanel/` | Show representative examples, not full catalogs |
+
+Prefer **short examples from tests** over hand-written snippets that drift:
+
+```python
+# Pattern: tests/pr/test_xarray_mask.py → docs/source/mask/xarray_mask_tutorial.ipynb
+```
+
+### 6.5 Build and review
+
+Local build:
+
+```bash
+cd docs && make html
+# open _build/html/index.html
+```
+
+Before merging the documentation branch:
+
+1. `make html` completes without warnings for missing `:doc:` references.
+2. New pages appear in the correct toctree (sidebar).
+3. Notebooks render (committed outputs present; `nb_execution_mode` is `off`).
+4. autoapi pages generate for new modules.
+
+Future CI enhancements (optional):
+
+- Sphinx `-W` (warnings as errors) on PRs touching `docs/`.
+- Link check for internal `:doc:` and relative md links.
+- Script to diff `list_datasets()` output against `overview.rst` mentions.
+
+---
+
+## 7. Immediate backlog for the documentation branch
+
+Actionable items in recommended order.
+
+### P0 — Navigation and broken links
+
+1. ~~**Add workflow chooser**~~ — **Done:** homepage (`intro.rst`) is the modern workflow; legacy content lives under **Legacy workflow** (`legacy/workflow.rst`).
+2. ~~**Ensure `xarray_mask_tutorial.ipynb` exists**~~ — **Done** (`docs/source/mask/xarray_mask_tutorial.ipynb`, included via mask `*` toctree).
+3. ~~**Update `intro.rst` masking section**~~ — **Done:** modern intro uses `XarrayMask`; Cutout masking unchanged in `legacy/workflow.rst`.
+
+### P0 — Sync with recent source changes
+
+4. ~~**`modeling/pvlib/index.rst`** — document `compact_output`~~ — **Done**.
+5. ~~**`modeling/wind/index.rst` and interpolation.rst** — document flexible `xs`/`ys`~~ — **Done**.
+6. ~~**`datasets/era5.rst`** — clarify CDS download vs offline fixtures~~ — **Done:** overview has recommended download; era5.rst covers CDS setup + optional cdsapi verify.
+
+### P1 — Structure and depth
+
+7. ~~**Wind CF documentation**~~ — **Done (Option A):** expanded Step 5 in interpolation/extrapolation; no separate internals page.
+8. ~~**Reorganize mask toctree intent**~~ — **Done:** Mask = creation + apply tutorials + troubleshoot; plans under Development.
+9. **`mask/merge_layer_known_issues.md`** — publish under mask with troubleshooting cross-links.
+10. **Link fixture doc from modeling tutorials** — one paragraph + code using `load_dataset("wind_3d_hourly_test")`.
+
+### P2 — Hygiene
+
+11. Move `docs/source/mask/create_mask.py`, `split_china.py`, etc. to `docs/source/mask/examples/` (exclude from glob toctree or document as examples).
+12. Fix README documentation URL placeholder.
+13. Audit `input_output.md` “planned” notes against `list_datasets()`.
+14. Add **this plan** to `index.rst` Development toctree (done when this file is merged).
+
+---
+
+## 7.2 Changelog expectations
+
+Documentation PRs should summarize:
+
+- **User-visible:** what readers can now do or what corrected behavior is documented.
+- **Structural:** new pages, moved pages, deprecated paths.
+- **Not required:** typo fixes only.
+
+For paired code+doc releases, use a single changelog entry covering both.
+
+---
+
+## 8. Long-term governance
+
+### 8.1 Ownership (suggested)
+
+| Area | Default maintainer focus |
+|------|--------------------------|
+| Datasets / ERA5 fixtures | Whoever changes `datasets/era5/` |
+| Wind / PV modeling | Model module authors |
+| Mask / XarrayMask | Mask package authors |
+| Legacy Cutout/convert | Touch only when behavior changes; avoid new features here |
+
+### 8.2 Deprecation policy for docs
+
+When deprecating APIs:
+
+1. Mark in docstring + autoapi.
+2. Add “Deprecated” admonition in legacy tutorial.
+3. Record timeline in a `development/` plan or release notes.
+4. Remove legacy tutorial sections only after code removal or major version bump.
+
+### 8.3 Quarterly doc audit (lightweight)
+
+Every ~3 months or before a release:
+
+1. Run `list_datasets()` and compare to `datasets/overview.rst`.
+2. Scan `intro.rst` for legacy-only examples without modern pointers.
+3. Grep docs for `planned`, `TODO`, `FIXME`.
+4. Confirm `make html` clean build.
+5. Update [Section 2.3](#23-known-gaps-as-of-this-plan) gap table in this file.
+
+### 8.4 Relationship to API reference
+
+autoapi is the **reference layer**; tutorials should not duplicate every
+argument. Convention:
+
+- Tutorials: one worked example with common options.
+- Reference: full signatures via docstrings (NumPy style, `sphinx.ext.napoleon`).
+- Explanation pages: algorithms and data flow (e.g. wind CF pipeline).
+
+When adding a feature, **docstring first**, then **one tutorial paragraph** —
+not a third full copy in markdown.
+
+---
+
+## 9. Appendix: proposed `index.rst` Development toctree
+
+```rst
+.. toctree::
+   :maxdepth: 1
+   :caption: Development
+   :hidden:
+
+   development/documentation-organization-plan
+   development/offline-era5-fixture-datasets
+```
+
+Mask migration plan remains under `mask/` glob but should use a contributor
+banner (see [Section 3.1](#31-page-roles-diátaxis)).
+
+---
+
+## 10. Appendix: doc branch merge strategy
+
+1. **Land structure first** (toctrees, workflow chooser, intro labels) so follow-up edits have a home.
+2. **Land content sync** (modeling/mask/datasets factual updates) in the same branch or stacked PRs by area.
+3. **Avoid** mixing large narrative rewrites with unrelated code changes — keeps review focused.
+4. After merge, tag a docs release note listing: new XarrayMask path, fixture-based tutorials, deprecated/legacy labeling.
+
+---
+
+## Document history
+
+| Date | Change |
+|------|--------|
+| 2026-06-02 | Initial organization plan for documentation branch |
diff --git a/docs/source/mask/mask_xarray_migration_plan.md b/docs/source/development/mask_xarray_migration_plan.md
similarity index 97%
rename from docs/source/mask/mask_xarray_migration_plan.md
rename to docs/source/development/mask_xarray_migration_plan.md
index 56ee8247..7b3afb70 100644
--- a/docs/source/mask/mask_xarray_migration_plan.md
+++ b/docs/source/development/mask_xarray_migration_plan.md
@@ -1,5 +1,10 @@
 # Mask-Without-Cutout Migration Plan
 
+```{note}
+**Audience:** contributors and maintainers. For applying saved masks to model
+output, see [XarrayMask tutorial](../mask/xarray_mask_tutorial.ipynb).
+```
+
 ## Goal
 
 Replace Cutout-dependent masking with a direct xarray-based workflow:
diff --git a/docs/source/mask/xarray_mask_workflow.rst b/docs/source/development/xarray_mask_workflow.rst
similarity index 86%
rename from docs/source/mask/xarray_mask_workflow.rst
rename to docs/source/development/xarray_mask_workflow.rst
index c4a10c4d..929a115f 100644
--- a/docs/source/mask/xarray_mask_workflow.rst
+++ b/docs/source/development/xarray_mask_workflow.rst
@@ -1,6 +1,12 @@
 Xarray masking workflow
 =========================
 
+.. note::
+
+   **Audience:** contributors and maintainers. This page records the xarray-first
+   masking implementation phases. For usage, see
+   :doc:`/mask/xarray_mask_tutorial`.
+
 This page summarizes the **xarray-first masking** work added alongside the
 longer-term plan in :doc:`mask_xarray_migration_plan`. The legacy path based on
 ``Cutout`` (``add_mask``, ``add_grid_area``, ``mask``) is unchanged for now; the
@@ -57,7 +63,7 @@ APIs. The intended usage is:
 2. Build ``XarrayMask.from_name("my_mask", grid=output_ds, mask_dir=...)`` if needed.
 3. Call ``attach(output_ds)`` or ``apply(output_ds, ...)`` for analysis.
 
-See :doc:`xarray_mask_tutorial` for a step-by-step notebook, and the offline
+See :doc:`/mask/xarray_mask_tutorial` for a step-by-step notebook, and the offline
 tests under ``tests/pr/`` (e.g. ``test_xarray_mask.py``,
 ``test_wind_xarraymask_integration.py``) for concrete examples.
 
@@ -77,6 +83,7 @@ continue to work during this transition.
 See also
 --------
 
+* :doc:`/mask/xarray_mask_tutorial` — step-by-step notebook (offline runnable).
 * :doc:`mask_xarray_migration_plan` — full migration phases and deprecation plan.
-* :doc:`mask_on_cutout` — legacy notebook: masks via ``Cutout``.
-* :doc:`mask_creation_workflow` — building and saving ``Mask`` objects from rasters.
+* :doc:`/legacy/mask_on_cutout` — legacy notebook: masks via ``Cutout``.
+* :doc:`/mask/mask_creation_workflow` — building and saving ``Mask`` objects from rasters.
diff --git a/docs/source/index.rst b/docs/source/index.rst
index 9c19e82d..a72c54c0 100644
--- a/docs/source/index.rst
+++ b/docs/source/index.rst
@@ -17,13 +17,18 @@ Welcome to Geodata's documentation!
    quick_start/input_output
 
 .. toctree::
-   :caption: Dataset Specific Tutorials
+   :maxdepth: 1
+   :caption: Legacy workflow
+   :hidden:
+
+   legacy/index
+
+.. toctree::
+   :caption: Datasets
    :maxdepth: 1
    :glob:
    :hidden:
 
-   datasets/era5/index
-   datasets/merra2/index
    datasets/*
 
 .. toctree::
@@ -31,16 +36,19 @@ Welcome to Geodata's documentation!
    :caption: Modeling
    :hidden:
 
+   modeling/era5_outputs
    modeling/wind/index
    modeling/pvlib/index
 
 .. toctree::
    :maxdepth: 1
    :caption: Mask
-   :glob:
    :hidden:
 
-   mask/*
+   mask/mask_creation_workflow
+   mask/xarray_mask_tutorial
+   mask/mask_troubleshoot
+   mask/merge_layer_known_issues
 
 .. .. toctree::
 ..    :maxdepth: 1
@@ -70,7 +78,10 @@ Welcome to Geodata's documentation!
    :caption: Development
    :hidden:
 
+   development/documentation-organization-plan
    development/offline-era5-fixture-datasets
+   development/mask_xarray_migration_plan
+   development/xarray_mask_workflow
 
 .. toctree::
    :maxdepth: 1
diff --git a/docs/source/intro.rst b/docs/source/intro.rst
index 784349c9..f4fc47e4 100644
--- a/docs/source/intro.rst
+++ b/docs/source/intro.rst
@@ -7,7 +7,7 @@ engineering, and social science applications.
 .. figure:: _static/images/geodata_workflow_chart.png
    :alt: Geodata Workflow
 
-   A typical anaylsis workflow with Geodata
+   A typical analysis workflow with Geodata
 
 Motivation
 ----------
@@ -31,227 +31,130 @@ model inputs. Additionally, with a minimal amount of data consistency
 checks and metadata information, when one researcher goes through this
 exercise, everyone benefits.
 
-How To Use
+How to use
 ----------
 
-Download Datasets
-~~~~~~~~~~~~~~~~~
-
-Earth system datasets can be large (100+ MB / file with hundreds of
-files necessary for a single analysis) and their APIs and file
-structures (e.g., daily vs monthly) vary by source. Utilizing xarray and
-dask data parallelization, Geodata provides single call download with
-API credentials stored locally. Data requests are automatically trimmed
-to keep only required variables, significantly reducing bandwidth
-requirements and disk usage.
-
-Geodata currently supports MERRA-2 and ERA5 reanalysis products and
-various GIS file formats (see :doc:`here<quick_start/input_output>`).
-
-**Note**:
-If you are exploring ERA5 data with wind or pvlib model, please refer to the :doc:`modeling/wind/index` and :doc:`modeling/pvlib/index` pages for more details.
-As they are following the dataset module to download data, not the following legacy code.
+Overview
+~~~~~~~~
 
-For example, to evaluate solar PV availability using
-`MERRA2 <https://disc.gsfc.nasa.gov/datasets/M2T1NXRAD_5.12.4/summary&sa=D&source=docs&ust=1646032882397245&usg=AOvVaw1WJVA5CNhT-7x_XWIA58o5>`__
-on 01/01/2011, use the following method call:
+The recommended workflow follows four steps:
 
-.. code :: Python
+1. **Load and download** a registered dataset with ``load_dataset``.
+2. **Run a model** (wind or solar PV) to produce xarray outputs.
+3. **Apply a mask** (optional) with ``XarrayMask`` on model output.
+4. **Analyze or visualize** the results in xarray, pandas, or with
+   ``geodata.plot``.
 
-   from geodata import Dataset
+Geodata supports ERA5 reanalysis through ``load_dataset`` and common GIS
+formats (see :doc:`quick_start/input_output`). For dataset setup and configs,
+see :doc:`datasets/overview`. MERRA-2 cutout workflows are documented under
+:doc:`/legacy/index`.
 
-   solar = Dataset(
-      module="merra2", 
-      years= slice(2011, 2011),
-      months=slice(1,1),
-      weather_data_config="slv_radiation_hourly"
-   )
-   solar.get_data()
+.. note::
 
-Extract Cutouts
-~~~~~~~~~~~~~~~
+   If you rely on the older ``Dataset`` / ``Cutout`` / ``convert`` API,
+   see :doc:`legacy/workflow`.
 
-Most energy analyses (e.g., energy models, resource assessments,
-political economy studies) require time series on subsets of locations
-and time periods. Geodata can extract desired variables, time periods,
-and geographies from the dataset to a Cutout object. We then call various functions in
-``geodata.convert`` module to transform the raw data into analysis-ready
-variables with the option to export to CSV or combine with other GIS
-datasets through further masking analysis.
+Step 1: Load and download a dataset
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-After downloading the required
-`MERRA2 <https://disc.gsfc.nasa.gov/datasets/M2T1NXRAD_5.12.4/summary&sa=D&source=docs&ust=1646032882397245&usg=AOvVaw1WJVA5CNhT-7x_XWIA58o5>`__
-dataset, we create a Cutout object that contains solar irradiance over
-China.
+Earth system datasets can be large (100+ MB per file, with many files
+per analysis). The ``geodata.datasets`` module provides a unified
+interface: pick a registered config, instantiate the dataset class, and
+download only the variables and time range you need.
 
 .. code :: Python
 
-   from geodata import Cutout
+   from geodata.datasets import load_dataset
 
-   cutout = Cutout(
-      name="china-2011-slv-hourly-test",
-      module="merra2",
-      weather_data_config="slv_radiation_hourly",
-      xs=slice(73, 136),
-      ys=slice(18, 54),
-      years=slice(2011, 2011),
+   ds_cls = load_dataset("wind_3d_hourly")
+   ds = ds_cls(
+      years=slice(2016, 2016),
       months=slice(1, 1),
+      bounds=[-10, 35, 10, 45],  # optional bounding box
    )
-   cutout.prepare()
-
-
-Then, we can convert the downward-shortwave, upward-shortwave radiation
-flux, and ambient temperature variables from the Cutout data into a PV
-generation time-series using the geodata ``convert`` method. Geodata
-stores objects internally as an xarray DataArray, which can be easily
-converted to a Pandas DataFrame.
 
-.. code :: Python
-
-   from geodata import convert
+   if not ds.downloaded:
+      ds.download()
 
-   ds_solar = convert.pv(cutout, panel="KANEKA", orientation="latitude_optimal")
-   ds_solar.to_dataframe(name="pv")
+   print(ds.downloaded)
 
+Use ``list_datasets()`` to see all registered configs. For ERA5 CDS
+credentials and offline test fixtures, see :doc:`datasets/era5` and
+:doc:`development/offline-era5-fixture-datasets`.
 
-.. figure:: _static/images/example_output_dataframe.png
-   :alt: Output DataFrame
-   :scale: 50%
+Step 2: Run a model
+~~~~~~~~~~~~~~~~~~~
 
-   Output of the code above
+Models operate on downloaded datasets and return **xarray** objects.
+Import the model explicitly (models are not re-exported at the top-level
+``geodata`` namespace).
 
-We can plot a time series of average PV values for all grid cells on
-that day with geodata's visualization method:
+**Wind** — interpolate or extrapolate hub-height wind speed and capacity
+factor from ERA5 3D wind data:
 
 .. code :: Python
 
-   from geodata import plot
-
-   plot.time_series(ds_solar)
+   from geodata.model.wind import WindInterpolationModel
 
-.. figure:: _static/images/visualization/output_12_0.png
-   :alt: Time-Series Plot
+   model = WindInterpolationModel(ds)
+   model.prepare()
+   wind_speed = model.estimate(height=100.0)
 
-   Visualization of the average PV values over time
+See :doc:`modeling/wind/index` for wind interpolation and turbine capacity factor, and
+turbine capacity-factor details.
 
-We can also visualize the average solar PV for every two hours on this
-day through an animation:
+**Solar PV** — estimate AC power and capacity factor with pvlib-backed
+models on ERA5 wind/solar hourly data:
 
 .. code :: Python
 
-   import geopandas as gpdø
-
-   from geodata import plot
-
-   prov_shapes = gpd.read_file(prov_shapes_path)
-   geodata.plot.heatmap_animation(
-      ds_solar,
-      cmap="Wistia",
-      time_factor=2,
-      shape=prov_shapes,
-      shape_width=0.25,
-      shape_color="navy",
-   )
-
-
-.. figure:: _static/images/visualization/pv_animation.gif
-   :alt: animation
+   from geodata.datasets import load_dataset
+   from geodata.model.pvlib import Pvlib
 
-   Animated Result
+   solar_cls = load_dataset("wind_solar_hourly")
+   solar_ds = solar_cls(years=slice(2016, 2016), months=slice(1, 1))
+   if not solar_ds.downloaded:
+      solar_ds.download()
 
-Masking
-~~~~~~~
+   pv_model = Pvlib(solar_ds)
+   # configure pv_system and model config — see modeling/pvlib/index
+   cf = pv_model.estimate(years=slice(2016, 2016), months=slice(1, 1))
 
-Geographic masks help filter datasets for specific analyses. Geodata is
-able to process GIS datasets and extract cutouts over specified
-geographies. Built off the open-source binary libraries GDAL, GEOS, and
-PROJ, and Python libraries rasterio and shapely, the Mask module imports
-rasters and shapefiles, edits them as mask layers, merges and flattens
-multiple layers together, and extracts subsetted cutout data from merged
-masks and shapefiles.
+See :doc:`modeling/pvlib/index` for full PV system and ModelChain setup.
 
-For example, within Geodata the user can load the `MODIS land use
-dataset <https://developers.google.com/earth-engine/datasets/catalog/MODIS_006_MCD12Q1#bands>`__,
-the `elevation
-dataset <https://developers.google.com/earth-engine/datasets/catalog/CGIAR_SRTM90_V4?hl=en>`__,
-and `environmental protected
-shapes <https://www.protectedplanet.net/country/CHN>`__, filter these
-according to solar energy suitability criteria, and merge into a single
-binary siting mask, where values of 0 represent the unsuitable area, and
-values of 1 represent the suitable area. Masks can be saved locally for
-later use.
+Step 3: Apply a mask (optional)
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-Geodata automatically reprojects GIS data in different coordinate
-reference systems into degree coordinates for processing. Common
-manipulations include cropping, filtering on categorical values,
-filtering on thresholds, excluding small contiguous areas, and filtering
-by shape buffers. One multi-purpose plotting function (``mask.show``)
-supports visualizing the mask including relevant shape boundaries.
-
-For example, Geodata can create a binary mask of wind energy suitability
-in China based on the above GIS inputs.
+Mask **creation** uses ``geodata.Mask`` (see
+:doc:`mask/mask_creation_workflow`). To apply a saved mask to model
+output without a ``Cutout``, use ``XarrayMask``:
 
 .. code :: Python
 
-   import geopandas as gpd
-
-   from geodata import mask
-
-   china = mask.Mask("China")
-   china.add_layer(layer_path={"modis": modis_path, "elevation": elevation_path})
-
-   protected_area_shapes = gpd.read_file(protected_area_shapes_path)
-   china.add_shape_layer(
-      protected_area_shapes["geometry"].to_dict(),
-      reference_layer="elevation",
-      combine_name="protected",
-      buffer=20,
-   )
-
-   china.filter_layer(
-      "modis", binarize=True, values=[6, 7, 8, 9, 10, 11, 12, 14, 15, 16, 17]
-   )
-   china.filter_layer("elevation", binarize=True, max_bound=4000)
-   china.merge_layer(trim=True)
-
-   china_prov_shapes = gpd.read_file(china_prov_shapes_path)
-   mask.show(china.merged_mask, shape=china_prov_shapes["geometry"], title="Merged Mask")
+   from geodata import XarrayMask
 
-   china.save_mask()
+   xmask = XarrayMask.from_name("my_mask", grid=wind_speed)
+   masked = xmask.apply(wind_speed, mode="where")
 
-.. figure:: _static/images/mask_workflow.png
-   :alt: mask workflow
+See :doc:`mask/xarray_mask_tutorial` for ``attach``, ``apply``, and
+grid-area weighting.
 
-   Visualization of Mask Workflow
+Step 4: Visualize
+~~~~~~~~~~~~~~~~~
 
-In the final step, we apply the Mask object to the Cutout. Geodata
-automatically coarsens the (typically) high-resolution Mask into the
-same resolution as the Cutout, adding fractions of the coarse cells
-covered by the Mask and areas calculated via an equal-area projection.
+Plotting works on any xarray object returned by a model:
 
 .. code :: Python
 
-   ds_cutout = convert.pv(
-      cutout, panel="KANEKA", orientation="latitude_optimal"
-   ).to_dataset(name="solar")
-
-   cutout.add_mask("china")
-   cutout.add_grid_area()
-   ds_mask = cutout.mask(dataset=ds_cutout)["merged_mask"]
-
-   weighted_mean_pv_series = (
-      (ds_mask["solar"] * ds_mask["mask"] * ds_mask["area"]).sum(axis=1).sum(axis=1)
-   ) / (ds_mask["mask"] * ds_mask["area"]).sum()
-
-   plt.plot(weighted_mean_pv_series)
-
+   from geodata import plot
 
-.. figure:: _static/images/mask_cutout_workflow.png
-   :alt: Mask-Cutout Workflow
+   plot.time_series(wind_speed)
 
-   Mask-Cutout Workflow
+See :doc:`visualization/visualization` for heatmaps and animations.
 
 What's next?
 ============
 
-To further explore the capabilities of Geodata, check out the table of contents on the left!
+Use the table of contents on the left to go deeper into datasets,
+modeling, masking, and the API reference.
diff --git a/docs/source/legacy/index.rst b/docs/source/legacy/index.rst
new file mode 100644
index 00000000..1a57f52a
--- /dev/null
+++ b/docs/source/legacy/index.rst
@@ -0,0 +1,25 @@
+Legacy workflow
+===============
+
+The pages below document the original Geodata API built around
+``Dataset``, ``Cutout``, ``geodata.convert``, and Cutout-based masking,
+including MERRA2 download and cutout tutorials.
+
+.. note::
+
+   This path is **not** part of the current tested workflow
+   (``load_dataset`` → models → ``XarrayMask``). It remains available for
+   existing analyses and reference.
+
+For the recommended path, see the :doc:`documentation homepage </index>`.
+
+.. toctree::
+   :maxdepth: 1
+
+   workflow
+   mask_on_cutout
+   merra2/index
+   merra2/merra2_download
+   merra2/merra2_outputs
+   merra2/merra2
+   wind_extrapolation
diff --git a/docs/source/mask/mask_on_cutout.ipynb b/docs/source/legacy/mask_on_cutout.ipynb
similarity index 100%
rename from docs/source/mask/mask_on_cutout.ipynb
rename to docs/source/legacy/mask_on_cutout.ipynb
diff --git a/docs/source/datasets/merra2/index.md b/docs/source/legacy/merra2/index.md
similarity index 77%
rename from docs/source/datasets/merra2/index.md
rename to docs/source/legacy/merra2/index.md
index fb020ff4..d571067a 100644
--- a/docs/source/datasets/merra2/index.md
+++ b/docs/source/legacy/merra2/index.md
@@ -1,5 +1,11 @@
 # MERRA2 Related Tutorials
 
+```{note}
+**Legacy documentation.** These tutorials use the older ``Dataset`` / ``Cutout`` API and are
+not part of the current tested workflow. For the recommended ERA5 path, see
+[Dataset module overview](../../datasets/overview.rst) and [ERA5 setup](../../datasets/era5.rst).
+```
+
 This page explains how you can setup access MERRA2 data from NASA's [GES DISC](https://disc.gsfc.nasa.gov/).
 
 ## Creating an Earthdata Login Profile and Approving the GES DISC App
@@ -41,9 +47,12 @@ For Windows, open Notepad and enter the following line in a new document, making
 
 Save the file to `C:\Users\<username>\.netrc`
 
-## What' next?
+## What's next?
 
 Now that you have configured your Earthdata Login credentials, you have successfully set up access to the MERRA-2 data.
-Please subsequently refer to the [general documentation on datasets](../overview.rst)
-for more information on how to download ERA5-based datasets using the `geodata`
-package.
+
+* [Download MERRA2 data and create cutouts](merra2_download.md)
+* [MERRA2 outputs via `convert`](merra2_outputs.md)
+* [MERRA2 workflow notebook](merra2.ipynb)
+
+For the current ERA5 + `load_dataset` workflow, see [Dataset module overview](../../datasets/overview.rst).
diff --git a/docs/source/datasets/merra2/merra2.ipynb b/docs/source/legacy/merra2/merra2.ipynb
similarity index 100%
rename from docs/source/datasets/merra2/merra2.ipynb
rename to docs/source/legacy/merra2/merra2.ipynb
diff --git a/docs/source/datasets/merra2/merra2_download.md b/docs/source/legacy/merra2/merra2_download.md
similarity index 100%
rename from docs/source/datasets/merra2/merra2_download.md
rename to docs/source/legacy/merra2/merra2_download.md
diff --git a/docs/source/datasets/merra2/merra2_outputs.md b/docs/source/legacy/merra2/merra2_outputs.md
similarity index 100%
rename from docs/source/datasets/merra2/merra2_outputs.md
rename to docs/source/legacy/merra2/merra2_outputs.md
diff --git a/docs/source/legacy/wind_extrapolation.rst b/docs/source/legacy/wind_extrapolation.rst
new file mode 100644
index 00000000..81dab3be
--- /dev/null
+++ b/docs/source/legacy/wind_extrapolation.rst
@@ -0,0 +1,102 @@
+Wind extrapolation (legacy)
+===========================
+
+.. note::
+
+   **Legacy / untested in CI.** ``WindExtrapolationModel`` only supports the
+   ``slv_flux_hourly`` weather config (MERRA-2 via ``load_dataset``). It is **not**
+   part of the current ERA5 workflow documented on the homepage. For ERA5 wind, use
+   :doc:`/modeling/wind/interpolation` instead.
+
+For the recommended modern path, see the :doc:`documentation homepage </index>`.
+
+Tutorial: Estimate Wind Speed with Extrapolation
+------------------------------------------------
+
+In this tutorial, we will learn how to estimate wind speed using the extrapolation model
+from the geodata library.
+
+.. warning::
+
+   Extrapolation requires a dataset with wind speed at **multiple** heights. In Geodata,
+   only ``slv_flux_hourly`` (MERRA-2) is registered for ``WindExtrapolationModel``.
+   Using any other ``weather_config`` raises ``ValueError``.
+
+Step 1: Import the necessary libraries
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+.. code:: Python
+
+    import xarray as xr
+
+    from geodata.datasets import load_dataset
+    from geodata.model.wind import WindExtrapolationModel
+
+
+Step 2: Load the dataset
+~~~~~~~~~~~~~~~~~~~~~~~~
+
+Use the MERRA-2 ``slv_flux_hourly`` config (not ERA5):
+
+.. code:: Python
+
+    ds_cls = load_dataset("slv_flux_hourly")
+    ds = ds_cls(
+        years=slice(2006, 2006),
+        months=slice(1, 1),
+        bounds=[-10, 35, 10, 45],
+    )
+
+    if not ds.downloaded:
+        ds.download()
+
+    print(ds.downloaded)
+
+
+Step 3: Compute extrapolation parameters
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+.. code:: Python
+
+    model = WindExtrapolationModel(ds)
+    model.prepare()
+
+Prepared coefficients are stored under ``GEODATA_ROOT/models/`` (see
+:doc:`/modeling/wind/index` — **Preparing the model** for ``prepare`` / ``prepared`` /
+``force``).
+
+Step 4: Estimate using the extrapolation model
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+.. code:: Python
+
+    estimated_wind_speed = model.estimate(
+        height=60,
+        years=slice(2006, 2006),
+        months=slice(1, 1),
+    )
+
+Step 5: Estimate wind turbine capacity factor
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+See :doc:`/modeling/wind/interpolation` Step 5 — **Understanding the output** for what
+``cf`` means. Example:
+
+.. code:: Python
+
+    estimated_cf = model.estimate(
+        turbine="Vestas_V112_3MW",
+        years=slice(2006, 2006),
+        months=slice(1, 1),
+    )
+
+
+How the Extrapolation Model Works
+---------------------------------
+
+The model calculates hub height wind speed from MERRA2 surface and low-level winds,
+extrapolating variables in MERRA's tavg1_2d_slv_Nx collection (2 m, 10 m, 50 m winds,
+displacement height, lowest model level winds, etc.).
+
+The hub height wind speed uses a log-profile fit (see the original tutorial in the
+repository history for the full equations).
diff --git a/docs/source/legacy/workflow.rst b/docs/source/legacy/workflow.rst
new file mode 100644
index 00000000..476fe909
--- /dev/null
+++ b/docs/source/legacy/workflow.rst
@@ -0,0 +1,227 @@
+Legacy workflow
+===============
+
+.. note::
+
+   This page documents the original ``Dataset`` → ``Cutout`` → ``convert`` workflow.
+   For the current recommended path, see the :doc:`documentation homepage </index>`.
+
+How To Use
+----------
+
+Download Datasets
+~~~~~~~~~~~~~~~~~
+
+Earth system datasets can be large (100+ MB / file with hundreds of
+files necessary for a single analysis) and their APIs and file
+structures (e.g., daily vs monthly) vary by source. Utilizing xarray and
+dask data parallelization, Geodata provides single call download with
+API credentials stored locally. Data requests are automatically trimmed
+to keep only required variables, significantly reducing bandwidth
+requirements and disk usage.
+
+Geodata currently supports MERRA-2 and ERA5 reanalysis products and
+various GIS file formats (see :doc:`here </quick_start/input_output>`).
+
+**Note**:
+If you are exploring ERA5 data with wind or pvlib model, please refer to the :doc:`/modeling/wind/index` and :doc:`/modeling/pvlib/index` pages for more details.
+As they are following the dataset module to download data, not the following legacy code.
+
+For example, to evaluate solar PV availability using
+`MERRA2 <https://disc.gsfc.nasa.gov/datasets/M2T1NXRAD_5.12.4/summary&sa=D&source=docs&ust=1646032882397245&usg=AOvVaw1WJVA5CNhT-7x_XWIA58o5>`__
+on 01/01/2011, use the following method call:
+
+.. code :: Python
+
+   from geodata import Dataset
+
+   solar = Dataset(
+      module="merra2", 
+      years= slice(2011, 2011),
+      months=slice(1,1),
+      weather_data_config="slv_radiation_hourly"
+   )
+   solar.get_data()
+
+Extract Cutouts
+~~~~~~~~~~~~~~~
+
+Most energy analyses (e.g., energy models, resource assessments,
+political economy studies) require time series on subsets of locations
+and time periods. Geodata can extract desired variables, time periods,
+and geographies from the dataset to a Cutout object. We then call various functions in
+``geodata.convert`` module to transform the raw data into analysis-ready
+variables with the option to export to CSV or combine with other GIS
+datasets through further masking analysis.
+
+After downloading the required
+`MERRA2 <https://disc.gsfc.nasa.gov/datasets/M2T1NXRAD_5.12.4/summary&sa=D&source=docs&ust=1646032882397245&usg=AOvVaw1WJVA5CNhT-7x_XWIA58o5>`__
+dataset, we create a Cutout object that contains solar irradiance over
+China.
+
+.. code :: Python
+
+   from geodata import Cutout
+
+   cutout = Cutout(
+      name="china-2011-slv-hourly-test",
+      module="merra2",
+      weather_data_config="slv_radiation_hourly",
+      xs=slice(73, 136),
+      ys=slice(18, 54),
+      years=slice(2011, 2011),
+      months=slice(1, 1),
+   )
+   cutout.prepare()
+
+
+Then, we can convert the downward-shortwave, upward-shortwave radiation
+flux, and ambient temperature variables from the Cutout data into a PV
+generation time-series using the geodata ``convert`` method. Geodata
+stores objects internally as an xarray DataArray, which can be easily
+converted to a Pandas DataFrame.
+
+.. code :: Python
+
+   from geodata import convert
+
+   ds_solar = convert.pv(cutout, panel="KANEKA", orientation="latitude_optimal")
+   ds_solar.to_dataframe(name="pv")
+
+
+.. figure:: ../_static/images/example_output_dataframe.png
+   :alt: Output DataFrame
+   :scale: 50%
+
+   Output of the code above
+
+We can plot a time series of average PV values for all grid cells on
+that day with geodata's visualization method:
+
+.. code :: Python
+
+   from geodata import plot
+
+   plot.time_series(ds_solar)
+
+.. figure:: ../_static/images/visualization/output_12_0.png
+   :alt: Time-Series Plot
+
+   Visualization of the average PV values over time
+
+We can also visualize the average solar PV for every two hours on this
+day through an animation:
+
+.. code :: Python
+
+   import geopandas as gpdø
+
+   from geodata import plot
+
+   prov_shapes = gpd.read_file(prov_shapes_path)
+   geodata.plot.heatmap_animation(
+      ds_solar,
+      cmap="Wistia",
+      time_factor=2,
+      shape=prov_shapes,
+      shape_width=0.25,
+      shape_color="navy",
+   )
+
+
+.. figure:: ../_static/images/visualization/pv_animation.gif
+   :alt: animation
+
+   Animated Result
+
+Masking
+~~~~~~~
+
+Geographic masks help filter datasets for specific analyses. Geodata is
+able to process GIS datasets and extract cutouts over specified
+geographies. Built off the open-source binary libraries GDAL, GEOS, and
+PROJ, and Python libraries rasterio and shapely, the Mask module imports
+rasters and shapefiles, edits them as mask layers, merges and flattens
+multiple layers together, and extracts subsetted cutout data from merged
+masks and shapefiles.
+
+For example, within Geodata the user can load the `MODIS land use
+dataset <https://developers.google.com/earth-engine/datasets/catalog/MODIS_006_MCD12Q1#bands>`__,
+the `elevation
+dataset <https://developers.google.com/earth-engine/datasets/catalog/CGIAR_SRTM90_V4?hl=en>`__,
+and `environmental protected
+shapes <https://www.protectedplanet.net/country/CHN>`__, filter these
+according to solar energy suitability criteria, and merge into a single
+binary siting mask, where values of 0 represent the unsuitable area, and
+values of 1 represent the suitable area. Masks can be saved locally for
+later use.
+
+Geodata automatically reprojects GIS data in different coordinate
+reference systems into degree coordinates for processing. Common
+manipulations include cropping, filtering on categorical values,
+filtering on thresholds, excluding small contiguous areas, and filtering
+by shape buffers. One multi-purpose plotting function (``mask.show``)
+supports visualizing the mask including relevant shape boundaries.
+
+For example, Geodata can create a binary mask of wind energy suitability
+in China based on the above GIS inputs.
+
+.. code :: Python
+
+   import geopandas as gpd
+
+   from geodata import mask
+
+   china = mask.Mask("China")
+   china.add_layer(layer_path={"modis": modis_path, "elevation": elevation_path})
+
+   protected_area_shapes = gpd.read_file(protected_area_shapes_path)
+   china.add_shape_layer(
+      protected_area_shapes["geometry"].to_dict(),
+      reference_layer="elevation",
+      combine_name="protected",
+      buffer=20,
+   )
+
+   china.filter_layer(
+      "modis", binarize=True, values=[6, 7, 8, 9, 10, 11, 12, 14, 15, 16, 17]
+   )
+   china.filter_layer("elevation", binarize=True, max_bound=4000)
+   china.merge_layer(trim=True)
+
+   china_prov_shapes = gpd.read_file(china_prov_shapes_path)
+   mask.show(china.merged_mask, shape=china_prov_shapes["geometry"], title="Merged Mask")
+
+   china.save_mask()
+
+.. figure:: ../_static/images/mask_workflow.png
+   :alt: mask workflow
+
+   Visualization of Mask Workflow
+
+In the final step, we apply the Mask object to the Cutout. Geodata
+automatically coarsens the (typically) high-resolution Mask into the
+same resolution as the Cutout, adding fractions of the coarse cells
+covered by the Mask and areas calculated via an equal-area projection.
+
+.. code :: Python
+
+   ds_cutout = convert.pv(
+      cutout, panel="KANEKA", orientation="latitude_optimal"
+   ).to_dataset(name="solar")
+
+   cutout.add_mask("china")
+   cutout.add_grid_area()
+   ds_mask = cutout.mask(dataset=ds_cutout)["merged_mask"]
+
+   weighted_mean_pv_series = (
+      (ds_mask["solar"] * ds_mask["mask"] * ds_mask["area"]).sum(axis=1).sum(axis=1)
+   ) / (ds_mask["mask"] * ds_mask["area"]).sum()
+
+   plt.plot(weighted_mean_pv_series)
+
+
+.. figure:: ../_static/images/mask_cutout_workflow.png
+   :alt: Mask-Cutout Workflow
+
+   Mask-Cutout Workflow
diff --git a/docs/source/mask/mask_creation_workflow.ipynb b/docs/source/mask/mask_creation_workflow.ipynb
index 636f6a1c..7d532403 100644
--- a/docs/source/mask/mask_creation_workflow.ipynb
+++ b/docs/source/mask/mask_creation_workflow.ipynb
@@ -1,1151 +1,1151 @@
 {
- "cells": [
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "# Tutorial: Typical Mask Workflow"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "## Introduction"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "Geodata is able to process geospatial data to extract cutouts over specified geographies. Built off the [rasterio library](https://rasterio.readthedocs.io/en/latest/quickstart.html), the **mask** module imports rasters and shapefiles, merges and flattens multiple layers together, and extracts subsetted cutout data from merged masks and shapefiles.\n",
-    "\n",
-    "Functionalities explored in this notebook:\n",
-    "\n",
-    "- [Creating a mask object, adding and manipulating layers](#mask-creation--adding-and-manipulating-layers)\n",
-    "- [Opening a shapefile and adding shape features as layers](#adding-shape-features-as-a-layer)\n",
-    "- [Merging and flattening layers](#merging-and-flattening-layers)\n",
-    "- [Eliminate small contiguous areas](#eliminate-small-contiguous-areas)\n",
-    "- [Extracting shapes from mask](#extracting-shapes-from-masks)\n",
-    "- [Saving and loading masks](#saving-and-loading-masks)"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "## Setup"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "To start, import the geodata package and required libraries. We can also import the `geodata.mask.show()` method for simplicity of its use."
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "import geopandas as gpd\n",
-    "import matplotlib.pyplot as plt\n",
-    "import numpy as np\n",
-    "import pandas as pd\n",
-    "\n",
-    "import geodata\n",
-    "from geodata.mask import show"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "Additionally, we use [cartopy](https://scitools.org.uk/cartopy/docs/latest/tutorials/using_the_shapereader.html#cartopy.io.shapereader.Reader) to download some common administrative region shapes, but user-provided shapefiles will also work:"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "import cartopy.io.shapereader as shpreader"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "### Shapefiles and Rasters\n",
-    "\n",
-    "We will use the following geotiff and shape files for this demo:\n",
-    "\n",
-    "\n",
-    "- `china_modis.tif`\n",
-    "\n",
-    "    We downloaded the MODIS land cover data, which uses satellite remote sensing data to estimate the land use type on an annual basis. See: [EarthData_MCD12Q1](https://lpdaac.usgs.gov/products/mcd12q1v006/).\n",
-    "\n",
-    "    We will use the IGBP classification ('LC_Type1') which has 17 different land use characterizations (the corresponding data thus takes values from 1.0 to 17.0).\n",
-    "    All the \"Bands\" are listed here: [Google_earth_engine_MODIS_006_MCD12Q1](https://developers.google.com/earth-engine/datasets/catalog/MODIS_006_MCD12Q1#bands)\n",
-    "\n",
-    "- `china_elevation.tif` and `china_slope.tif`\n",
-    "\n",
-    "    These two rasters are based on the elevation map from: [Google_earth_engine_MODIS_CGIAR_SRTM90_V4](https://developers.google.com/earth-engine/datasets/catalog/CGIAR_SRTM90_V4?hl=en). Slope was computed in degrees using the 4-connected neighbors of each pixel. \n",
-    "\n",
-    "\n",
-    "- `UNEP_WDPA_China` Shapefiles\n",
-    "\n",
-    "    We downloaded the environmental protected area for China from: [ProtectedPlanet_China](https://www.protectedplanet.net/country/CHN). These shapefiles are distributed among 3 subfolders upon successful download and decompression due to the large size. We will create path variables for all three subfolders and we will only take the polygon shapes.\n",
-    "\n",
-    "    Alternatively, We can also retrieve the environmental protected area from Google Earth Engine: [Google_earth_engine_WCMC_WDPA](https://developers.google.com/earth-engine/datasets/catalog/WCMC_WDPA_current_polygons). The shapefile will contain the protected shapes from entire world (and the size is slightly over 1 GB), and additional data cleaning will be necessary if the user wants just the shapes within China. "
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "modis_path = \"data/china_modis.tif\"\n",
-    "elevation_path = \"data/china_elevation.tif\"\n",
-    "slope_path = \"data/china_slope.tif\"\n",
-    "\n",
-    "wdpa_shape_path_0 = \"data/shapefiles/0/WDPA_WDOECM_Nov2021_Public_CHN_shp-polygons.shp\"\n",
-    "wdpa_shape_path_1 = \"data/shapefiles/1/WDPA_WDOECM_Nov2021_Public_CHN_shp-polygons.shp\"\n",
-    "wdpa_shape_path_2 = \"data/shapefiles/2/WDPA_WDOECM_Nov2021_Public_CHN_shp-polygons.shp\""
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "Let us get province shapes from `cartopy` and save the path as `prov_path`. This can also be the path to user-supplied shape files."
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "prov_path = shpreader.natural_earth(\n",
-    "    resolution=\"10m\", category=\"cultural\", name=\"admin_1_states_provinces\"\n",
-    ")\n",
-    "prov_path"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "Load the shapes contained in path `prov_path` using the `geopandas` library."
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "scrolled": true
-   },
-   "outputs": [],
-   "source": [
-    "all_shapes = gpd.read_file(prov_path, encoding=\"utf-8\")\n",
-    "all_shapes.head(2)"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "GeoPandas data filtering with GeoDataFrame is identical to pandas. Let us select all the rows that contains shape within China."
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "china_shapes = all_shapes[all_shapes[\"admin\"] == \"China\"]"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "Next, to load the WDPA environmental protected shapefiles as a layer in the china mask, we will use the GeoPandas library. `gpd.read_file()` will return a GeoPandas dataframe including shape attributes and geometry given the file path. Like Pandas, we can read multiple dataframes and concat them together. In the code below, we will create one GeoPandas dataframe from three paths that we have for the Chinese environmental protected shapes."
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "wdpa_shapes = pd.concat([\n",
-    "    gpd.read_file(wdpa_shape_path_0),\n",
-    "    gpd.read_file(wdpa_shape_path_1),\n",
-    "    gpd.read_file(wdpa_shape_path_2)\n",
-    "])\n",
-    "wdpa_shapes.head(2)"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "## Mask Creation & Adding and Manipulating Layers"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "The mask object consists of multiple layers and manipulations performed on them. To add a layer, the four methods below perform same functions. A user may add a layer to the mask by specifying paths when a new instance is created, or use the `add_layer` method. We will add the following two files: `china_elevation.tif`, and `china_modis.tif` to the `China` mask, and name them `elevation` and `modis` layers."
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "scrolled": true
-   },
-   "outputs": [],
-   "source": [
-    "# Method 1: Initialize one layer, add one layer\n",
-    "china = geodata.Mask(\"China\", layer_path=elevation_path)\n",
-    "china.rename_layer(\"china_elevation\", \"elevation\")\n",
-    "china.add_layer(modis_path, layer_name=\"modis\")"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "# Method 2:  Initialize empty, add two layers using dict\n",
-    "china = geodata.Mask(\"China\")\n",
-    "china.add_layer(layer_path={\"elevation\": elevation_path, \"modis\": modis_path})"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "# Method 3:  Initalize with two layers passed as list\n",
-    "china = geodata.Mask(\n",
-    "    \"China\", layer_path=[elevation_path, modis_path], layer_name=[\"elevation\", \"modis\"]\n",
-    ")"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "# Method 4:  Initialize with two layers passed as dict\n",
-    "china = geodata.Mask(\n",
-    "    \"China\", layer_path={\"elevation\": elevation_path, \"modis\": modis_path}\n",
-    ")"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "Display the mask object in the jupyter notebook:"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "china"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "Each mask object has several attributes:\n",
-    "- `layers`: a dictionary of name (key)  - rasterio file opener (values). The <\\open DatasetReader> can be the input for many other mask methods for the module.  \n",
-    "- `merged_mask`: the merged and flatten mask of its layers, the merged raster from `layers`\n",
-    "- `shape_mask`: similar to the `layers` attribute, but a dictionary of extracted shapes from the merged mask by default. Users may also extracted shape masks from specified layers in `self.layers`.\n",
-    "- `saved`: whether this mask object has been saved locally.\n",
-    "- `mask_dir`: the directory to save the mask object, by default it should be the mask dir in config.py."
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "Show the `slope` layer in mask `china`. The `show` method will always try to show the proper latitude and longitude, unless we call it `show(layer, lat_lon = False)`."
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "china.layers[\"elevation\"]"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "show(china.layers[\"elevation\"], title=\"Elevation of China in meters\")"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "china.layers"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "**Some useful methods to examine the layers**\n",
-    "\n",
-    "- `china.get_res()`: get resolution of each layer, in lat-lon coordinates\n",
-    "- `china.get_res(product = True)`: get grid cell size, in product of lat-lon coordinate differences\n",
-    "- `china.get_bounds()`: get bounds, in lat-lon coordinates"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "china.get_bounds()"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "Note that the modis layer has a very different bounding box then the slope layer in lat-lon coordinate system. This is because the modis layer was converted to the lat-lon CRS from a different CRS when it was added to the object. The following section will explore CRS conversion."
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "### CRS conversion, trimming, and cropping (Optional)"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "Method `open_tif` can open a layer without adding it to the layer, this allows us to visualize it before-hand. It is a good practice to close the raster after opening it to avoid writing permission conflict issues. Closing the raster below does not involve any layer operation associated with the mask object. "
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "modis_opener = geodata.mask.open_tif(modis_path, show_raster=True)\n",
-    "modis_opener.close()"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "We can use `remove_layer` method to remove a layer to mask `china`. This method will properly close the raster file, because the raster file would remain open after being added to the mask."
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "china.remove_layer(\"modis\")"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "The `add_layer` method incorporates coordinate reference system (CRS) conversion to lat-lon (EPSG:4326), if necessary. Note that this method will overwrite the layer by default, if it is in the object already, unless the user specifies `replace=False`. \n",
-    "\n",
-    "The method will automatically trim the all-zero columns/rows. By default, the paramater `trim` is set to `True`. If we do not set it to True, we might generate a converted raster with new CRS but many all-zero columns and rows."
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "china.add_layer(modis_path, \"modis\", trim=False)\n",
-    "show(china.layers[\"modis\"], title=\"China Modis CRS converted (No trimming)\")"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "We can also crop a raster/layer with user-defined dimensions: method `crop_layer` can take either starting indices of top/left, ending indices of right/bottom, or coordinates values in lat/long to trim the raster.\n",
-    "\n",
-    "The difference between `crop_layer` and `trim_layer` is that `crop_layer` must take in user specified range to crop the raster, and `trim_layer` would remove the all zero rows and columns automatically for a raster. So that if the user do not know which index to start and end to remove the empty rows/columns, `trim_raster` is better.\n",
-    "\n",
-    "The method `crop_raster` (`geodata.mask.crop_raster`) is similar to `crop_layer` but can take a layer name as input, so that the user does not need to add a raster as a layer to call that method. (Similar method: `trim_layer`/`trim_raster`, `binarize_layer`/`binarize_raster`)"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "china.crop_layer(\"modis\", bounds=(73, 17, 135, 54))\n",
-    "show(china.layers[\"modis\"], title=\"China Modis Layer Cropped\")"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "This performs the same function by passing the layer to `crop_raster`:"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "china.layers[\"modis\"] = geodata.mask.crop_raster(\n",
-    "    china.layers[\"modis\"], (73, 17, 135, 54)\n",
-    ")"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "### Filter a layer"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "The mask module also supports filtering a layer based on list of categorical values, a minimum (lower) boundary, or maximum (upper) boundary.\n",
-    "\n",
-    "In the `filter_raster` method, a user may specify any of the `value` (the list of numberic values in the raster array to be selected), `max_bound`, and `min_bound` parameters to selected desired values. If the parameter `binarize` is False (by default), the method will return the original values of the raster that satisfy the conditions, otherwise the method will return 1 for the values that satisfy the conditions and 0 elsewhere."
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "#### Select Categorical Values from MODIS Layer\n",
-    "\n",
-    "Since the modis layer has 17 distinct values for different land use types, we want to create a layer of binary values, indicating unavailable land as 0, and available land as 1.\n",
-    "\n",
-    "We wish to create a mask where :\n",
-    "\n",
-    "- all forested areas (values 1-5) are 0 (i.e., unsuitable)\n",
-    "- all urban areas (13) are 0\n",
-    "- all others are 1\n",
-    "\n",
-    "\n",
-    "Let us use method `filter_raster` to create a layer of `modis_filtered` binary mask, where 1, 2, 3, 4, 5, and 13 will be unavailable land assigned 0 and the rest of the values will be 1 (available).\n"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "avail_values = list(set(range(1, 18)) - {1, 2, 3, 4, 5, 13})\n",
-    "avail_values"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "china.layers[\"modis_filtered\"] = geodata.mask.filter_raster(\n",
-    "    china.layers[\"modis\"], binarize=True, values=avail_values\n",
-    ")"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "scrolled": false
-   },
-   "outputs": [],
-   "source": [
-    "china.remove_layer(\"modis\")\n",
-    "show(china.layers[\"modis_filtered\"])"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "#### Filter elevation layer\n",
-    "\n",
-    "Because we cannot build renewable energy in areas with high elevation, let us set the constraint from the `elevation` layer, by using elevation < 4000m at 1 and other areas as 0. The result layer `elevation_filtered` will have only 1 and 0 as unique values."
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "china.filter_layer(\n",
-    "    \"elevation\", dest_layer_name=\"elevation_filtered\", max_bound=4000, binarize=True\n",
-    ")"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "china.remove_layer(\"elevation\")\n",
-    "show(china.layers[\"elevation_filtered\"])"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "#### Filter Slope Layer\n",
-    "\n",
-    "We also cannot build renewable energy in area with large slopes, so let us set another constraint from the `slope` layer from the slope tif file, by using slope < 20 degree at 1 and else as 0. The result layer `slope_filtered` will have only 1 and 0 as unique values."
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "First, add the slope raster to the china mask."
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "china.add_layer(slope_path, layer_name=\"slope\")\n",
-    "show(china.layers[\"slope\"])"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "Filter the raster, delete the old slope layer."
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "china.filter_layer(\n",
-    "    \"slope\", dest_layer_name=\"slope_filtered\", max_bound=20, binarize=True\n",
-    ")"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "scrolled": false
-   },
-   "outputs": [],
-   "source": [
-    "china.remove_layer(\"slope\")\n",
-    "show(china.layers[\"slope_filtered\"])"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "### Additional Visualization Options"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "We can plot the provinces on a selected layer by taking `shape` input in the `show()` method. Here, we will use the `china_shapes` that we obtained from `all_shape`. Its `geometry` column is a Series of shapes (shapely.geometry or MultiPolygon) for Chinese provinces."
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "show(china.layers[\"modis_filtered\"], shape=china_shapes[\"geometry\"])"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "## Adding Shape Features as a Layer"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "Recall that we have previously loaded the environmental protected shapes of China in a GeoPandas dataframe."
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "len(wdpa_shapes)"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "The three shapefiles have 78 features altogether, but we want to add all the features to one new layer instead of 78 new layers. The input shape should be a python dictionary, where there is a key for each unique shape. Also, in the `add_shape_layer` method, we will specify a `combine_name` to combine the features into one layer in this case, since we want the mask to have just one more layers, not 78 more layers.\n",
-    "\n",
-    "When adding a shapefile, we must specify the dimensions. We will also use `reference layer = 'slope_filtered'` so the new shape layer will have the same dimension with the `slope_filtered` layer. If the mask is empty and does not contain any layer, the user will have to specify the `resolution` parameter for the raster layer dimension.\n",
-    "\n",
-    "By default, this method will have paramater `exclude` that defaults to `False`. When it is true, area inside the shape is 0. When it is false, area inside the shape is 1. In this use case, however, we want 0 for area inside of the shape as they are environmental protected areas to exclude. We can just use the default method call."
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "china.add_shape_layer(\n",
-    "    wdpa_shapes[\"geometry\"].to_dict(),\n",
-    "    reference_layer=\"slope_filtered\",\n",
-    "    combine_name=\"protected\",\n",
-    ")\n",
-    "show(\n",
-    "    china.layers[\"protected\"],\n",
-    "    title=\"WDPA Protected area shape features as a new layer\",\n",
-    "    grid=True,\n",
-    ")"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "We can also use the parameter `buffer` in `add_shape_layer` method to create an approximate representation of all locations within a given (perpindicular) distance of the shape object. The units for the buffer are given in kilometers.\n",
-    "\n",
-    "Note that since the units of the original shape are in lat-lon coordinates, when we add the buffer, we will need to have a CRS that has meter as unit. The program will convert the shapes to that CRS, add the buffer around shapes, then convert it back to the lat-lon CRS system. By default, we used \"EPSG:6933\", an equal area projection CRS to add buffer in kilometer. "
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "km_buffer = 20\n",
-    "\n",
-    "china.add_shape_layer(\n",
-    "    wdpa_shapes[\"geometry\"].to_dict(),\n",
-    "    reference_layer=\"slope_filtered\",\n",
-    "    combine_name=\"protected_with_buffer\",\n",
-    "    buffer=km_buffer,\n",
-    ")\n",
-    "\n",
-    "show(\n",
-    "    china.layers[\"protected_with_buffer\"],\n",
-    "    title=f\"WDPA Protected area shape with {km_buffer}km buffer\",\n",
-    "    grid=True,\n",
-    ")\n",
-    "\n",
-    "china.remove_layer(\"protected_with_buffer\")"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "## Merging and Flattening Layers"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "In order to combine all layers into one, we use the `merge_layer` method which creates a new layer called `merged_mask`. This merges multiple layers together and flattens them using either **and** (default) or **sum** method, saving the result to `self.merged_mask` by default. Geospatial bounds and resolution of the output layer are in the units of the input file coordinate reference system, but by default, we will use the resolution of the layer with the best (finest) resolution for the output bounds/resolution, unless a reference layer is provided. In this case, the resolution of the merged_mask is the same with the `modis_filtered` layer. "
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "china.get_res()"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "china.merge_layer(attribute_save=False, show_raster=False).res"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "### Binary `AND` Method"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "By default, the `merge_layer` method will use a binary 'and' method: for each grid cell, if any of the n layers are 0, then the returned `self.merged_layer` will also have 0 at that location. In other words, if all the layers indicate that a land is available (!=0), the merged result will have value 1.\n",
-    "\n",
-    "`merge_layer` may also take in an optional parameter `layers`, which is a list of layer names stored in the object, if the user does not wish to merge all layers in the object. If the user does not want to save the result to the `merged_mask` attribute, the user can specify `attribute_save = False`."
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "# merge and plot only, do not save\n",
-    "china.merge_layer(attribute_save=False, layers=[\"slope_filtered\", \"modis_filtered\"])"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "china"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "Try again with the `reference_layer` parameter:"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "china.merge_layer(\n",
-    "    layers=[\"elevation_filtered\", \"modis_filtered\"],\n",
-    "    reference_layer=\"elevation_filtered\",\n",
-    "    show_raster=False,\n",
-    ")"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "The result of the `merged_mask` method is saved to `china.merged_mask` with the same resolution as the reference layer, in this case `elevation_filtered`."
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "china.merged_mask.res"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "For the purpose of this demonstration, we will select the `AND` method for the final merged_mask. We can also trim the border of the merged mask since the 4 layers have different boundaries. We can set the parameter `trim = True`."
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "china.merge_layer(trim=True)"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "### `SUM` Method"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "The sum method will add up the values from all the layers using weights. When there is no weight dict provided, all the layers for merging will have weights of 1 by default.\n",
-    "\n",
-    "Note: since we are not using the sum method to proceed to the following sections, we will keep `attribute_save = False` to prevent this method from overwriting the mask we have previously created above."
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "china.merge_layer(method=\"sum\", attribute_save=False, trim=True)"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "This distribution is completely arbitrary for the purpose of demonstration of the module: (Note: The weights do not need to have a total of 1)\n",
-    "- elevation_filtered: 0.15, slope_filtered: 0.1, modis_filtered: 0.3, protected: 0.45\n",
-    "\n",
-    "We will write the result to a new variable `customized_merged_layer` for continuing processing."
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "customized_merged_layer = china.merge_layer(\n",
-    "    method=\"sum\",\n",
-    "    weights={\n",
-    "        \"elevation_filtered\": 0.15,\n",
-    "        \"slope_filtered\": 0.1,\n",
-    "        \"modis_filtered\": 0.3,\n",
-    "        \"protected\": 0.45,\n",
-    "    },\n",
-    "    attribute_save=False,\n",
-    "    trim=True,\n",
-    ")"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "If the continuous value created by `merged_mask` represents a suitability metric, we could set a minimum value of 0.8 to be considered \"suitable\" (or 1). We then apply the `filter_raster` method on the merged layer."
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "customized_merged_layer = geodata.mask.filter_raster(\n",
-    "    customized_merged_layer, min_bound=0.8, binarize=True\n",
-    ")\n",
-    "show(customized_merged_layer)"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "## Eliminate Small Contiguous Areas"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "Using the above methods, we might end up with many small contiguous areas that are marked suitable but surrounded by an unsuitable region. We may want to exclude such regions from renewable energy development. The `filter_area` method will remove the small contiguous suitable regions by transforming the merged mask raster to polygons/shapes, calculating the area of each polygon, and filtering out polygons that are smaller than a given threshold. Units are given in kilometer-squared (km$^2$).\n",
-    "\n",
-    "By default, `filter_area` uses the merged mask raster and returns a new raster, unless input/output layers are specified by `layer_name` and `dest_layer_name`. \n",
-    "\n",
-    "By default, its `shape_value` parameter is 1, indicating that we are only interested in finding all groups of cells with value 1 (suitable) for elimination. We specify the threshold with the `min_area` parameter.\n",
-    "\n",
-    "Note: the `filter_area` method may take a long time (5 or more minutes depending on the complexity of your layer and your computational setup). The method relies upon `rasterio.rasterize`, see performance notes: https://rasterio.readthedocs.io/en/latest/api/rasterio.features.html#rasterio.features.rasterize\n"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "For example, if we focus on Guangdong province in Southern China from the merged mask, we notice that there are many small islands in the ocean that are marked as suitable areas. We want to exclude these small regions from our merged mask."
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "plt.imshow(china.merged_mask.read(1)[4800:5300, 5700:6600], interpolation=\"none\")\n",
-    "plt.show()"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "Call `filter_area` to remove all contiguous suitable region shapes smaller than 100 km$^2$:"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "china.merged_mask = geodata.mask.filter_area(china, min_area=100)"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "There shapes are removed in the new merged_mask."
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "plt.imshow(china.merged_mask.read(1)[4800:5300, 5700:6600], interpolation=\"none\")\n",
-    "plt.show()"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "## Extracting Shapes from Masks"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "Sometimes the user needs to generate masks and perform analysis for a collection of regions (e.g., at the state/province level). The purpose of shape extraction (`extract_shapes`) is to separate `merged_mask` values for each region, with the result a dictionary of name-mask pairs in the `shape_mask` attribute of the mask object. The values of `shape_mask` will be 0 outside of the shape, and will be `merged_mask` inside of the shape."
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "For the purpose of this demonstration, we will only select the province of Jiangsu, Zhejiang, and Shanghai."
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "scrolled": false
-   },
-   "outputs": [],
-   "source": [
-    "china_shapes_subset = china_shapes[\n",
-    "    china_shapes[\"name\"].isin([\"Jiangsu\", \"Zhejiang\", \"Shanghai\"])\n",
-    "]\n",
-    "china_shapes_subset"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "Converting the filtered shape dictionary to a python dictionary as the input for `extract_shapes`, where the keys for the dictionary will be the names of the new extracted shape layers."
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "china_shapes_subset = (\n",
-    "    china_shapes_subset[[\"name\", \"geometry\"]].set_index(\"name\")[\"geometry\"].to_dict()\n",
-    ")\n",
-    "china_shapes_subset"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "Extract the shapes from the merged_mask. "
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "china.extract_shapes(china_shapes_subset)"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "The resulting mask object contains the dictionary `shape_mask` with the extracted values:"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "china"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "## Saving and Loading Masks"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "china.save_mask()"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "With the mask saved, the user can now load the layers or shapes with `xarray` instead if preferred."
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "shape_xr_lst = china.load_shape_xr()\n",
-    "shape_xr_lst[\"Zhejiang\"].plot()"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "Optional: closing all the files when saving the mask. This can avoid possible write permission error."
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "china.save_mask(close_files=True)"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "Loading a previously saved mask."
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "china_2 = geodata.mask.load_mask(\"china\")"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "china_2"
-   ]
-  }
- ],
- "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.10.11"
-  }
- },
- "nbformat": 4,
- "nbformat_minor": 4
-}
+  "cells": [
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "# Tutorial: Typical Mask Creation Workflow"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "## Introduction"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "Geodata is able to process geospatial data to extract cutouts over specified geographies. Built off the [rasterio library](https://rasterio.readthedocs.io/en/latest/quickstart.html), the **mask** module imports rasters and shapefiles, merges and flattens multiple layers together, and extracts subsetted cutout data from merged masks and shapefiles.\n",
+        "\n",
+        "Functionalities explored in this notebook:\n",
+        "\n",
+        "- [Creating a mask object, adding and manipulating layers](#mask-creation--adding-and-manipulating-layers)\n",
+        "- [Opening a shapefile and adding shape features as layers](#adding-shape-features-as-a-layer)\n",
+        "- [Merging and flattening layers](#merging-and-flattening-layers)\n",
+        "- [Eliminate small contiguous areas](#eliminate-small-contiguous-areas)\n",
+        "- [Extracting shapes from mask](#extracting-shapes-from-masks)\n",
+        "- [Saving and loading masks](#saving-and-loading-masks)"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "## Setup"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "To start, import the geodata package and required libraries. We can also import the `geodata.mask.show()` method for simplicity of its use."
+      ]
+    },
+    {
+      "cell_type": "code",
+      "metadata": {},
+      "source": [
+        "import geopandas as gpd\n",
+        "import matplotlib.pyplot as plt\n",
+        "import numpy as np\n",
+        "import pandas as pd\n",
+        "\n",
+        "import geodata\n",
+        "from geodata.mask import show"
+      ],
+      "execution_count": null,
+      "outputs": []
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "Additionally, we use [cartopy](https://scitools.org.uk/cartopy/docs/latest/tutorials/using_the_shapereader.html#cartopy.io.shapereader.Reader) to download some common administrative region shapes, but user-provided shapefiles will also work:"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "metadata": {},
+      "source": [
+        "import cartopy.io.shapereader as shpreader"
+      ],
+      "execution_count": null,
+      "outputs": []
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "### Shapefiles and Rasters\n",
+        "\n",
+        "We will use the following geotiff and shape files for this demo:\n",
+        "\n",
+        "\n",
+        "- `china_modis.tif`\n",
+        "\n",
+        "    We downloaded the MODIS land cover data, which uses satellite remote sensing data to estimate the land use type on an annual basis. See: [EarthData_MCD12Q1](https://lpdaac.usgs.gov/products/mcd12q1v006/).\n",
+        "\n",
+        "    We will use the IGBP classification ('LC_Type1') which has 17 different land use characterizations (the corresponding data thus takes values from 1.0 to 17.0).\n",
+        "    All the \"Bands\" are listed here: [Google_earth_engine_MODIS_006_MCD12Q1](https://developers.google.com/earth-engine/datasets/catalog/MODIS_006_MCD12Q1#bands)\n",
+        "\n",
+        "- `china_elevation.tif` and `china_slope.tif`\n",
+        "\n",
+        "    These two rasters are based on the elevation map from: [Google_earth_engine_MODIS_CGIAR_SRTM90_V4](https://developers.google.com/earth-engine/datasets/catalog/CGIAR_SRTM90_V4?hl=en). Slope was computed in degrees using the 4-connected neighbors of each pixel. \n",
+        "\n",
+        "\n",
+        "- `UNEP_WDPA_China` Shapefiles\n",
+        "\n",
+        "    We downloaded the environmental protected area for China from: [ProtectedPlanet_China](https://www.protectedplanet.net/country/CHN). These shapefiles are distributed among 3 subfolders upon successful download and decompression due to the large size. We will create path variables for all three subfolders and we will only take the polygon shapes.\n",
+        "\n",
+        "    Alternatively, We can also retrieve the environmental protected area from Google Earth Engine: [Google_earth_engine_WCMC_WDPA](https://developers.google.com/earth-engine/datasets/catalog/WCMC_WDPA_current_polygons). The shapefile will contain the protected shapes from entire world (and the size is slightly over 1 GB), and additional data cleaning will be necessary if the user wants just the shapes within China. "
+      ]
+    },
+    {
+      "cell_type": "code",
+      "metadata": {},
+      "source": [
+        "modis_path = \"data/china_modis.tif\"\n",
+        "elevation_path = \"data/china_elevation.tif\"\n",
+        "slope_path = \"data/china_slope.tif\"\n",
+        "\n",
+        "wdpa_shape_path_0 = \"data/shapefiles/0/WDPA_WDOECM_Nov2021_Public_CHN_shp-polygons.shp\"\n",
+        "wdpa_shape_path_1 = \"data/shapefiles/1/WDPA_WDOECM_Nov2021_Public_CHN_shp-polygons.shp\"\n",
+        "wdpa_shape_path_2 = \"data/shapefiles/2/WDPA_WDOECM_Nov2021_Public_CHN_shp-polygons.shp\""
+      ],
+      "execution_count": null,
+      "outputs": []
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "Let us get province shapes from `cartopy` and save the path as `prov_path`. This can also be the path to user-supplied shape files."
+      ]
+    },
+    {
+      "cell_type": "code",
+      "metadata": {},
+      "source": [
+        "prov_path = shpreader.natural_earth(\n",
+        "    resolution=\"10m\", category=\"cultural\", name=\"admin_1_states_provinces\"\n",
+        ")\n",
+        "prov_path"
+      ],
+      "execution_count": null,
+      "outputs": []
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "Load the shapes contained in path `prov_path` using the `geopandas` library."
+      ]
+    },
+    {
+      "cell_type": "code",
+      "metadata": {
+        "scrolled": true
+      },
+      "source": [
+        "all_shapes = gpd.read_file(prov_path, encoding=\"utf-8\")\n",
+        "all_shapes.head(2)"
+      ],
+      "execution_count": null,
+      "outputs": []
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "GeoPandas data filtering with GeoDataFrame is identical to pandas. Let us select all the rows that contains shape within China."
+      ]
+    },
+    {
+      "cell_type": "code",
+      "metadata": {},
+      "source": [
+        "china_shapes = all_shapes[all_shapes[\"admin\"] == \"China\"]"
+      ],
+      "execution_count": null,
+      "outputs": []
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "Next, to load the WDPA environmental protected shapefiles as a layer in the china mask, we will use the GeoPandas library. `gpd.read_file()` will return a GeoPandas dataframe including shape attributes and geometry given the file path. Like Pandas, we can read multiple dataframes and concat them together. In the code below, we will create one GeoPandas dataframe from three paths that we have for the Chinese environmental protected shapes."
+      ]
+    },
+    {
+      "cell_type": "code",
+      "metadata": {},
+      "source": [
+        "wdpa_shapes = pd.concat([\n",
+        "    gpd.read_file(wdpa_shape_path_0),\n",
+        "    gpd.read_file(wdpa_shape_path_1),\n",
+        "    gpd.read_file(wdpa_shape_path_2)\n",
+        "])\n",
+        "wdpa_shapes.head(2)"
+      ],
+      "execution_count": null,
+      "outputs": []
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "## Mask Creation & Adding and Manipulating Layers"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "The mask object consists of multiple layers and manipulations performed on them. To add a layer, the four methods below perform same functions. A user may add a layer to the mask by specifying paths when a new instance is created, or use the `add_layer` method. We will add the following two files: `china_elevation.tif`, and `china_modis.tif` to the `China` mask, and name them `elevation` and `modis` layers."
+      ]
+    },
+    {
+      "cell_type": "code",
+      "metadata": {
+        "scrolled": true
+      },
+      "source": [
+        "# Method 1: Initialize one layer, add one layer\n",
+        "china = geodata.Mask(\"China\", layer_path=elevation_path)\n",
+        "china.rename_layer(\"china_elevation\", \"elevation\")\n",
+        "china.add_layer(modis_path, layer_name=\"modis\")"
+      ],
+      "execution_count": null,
+      "outputs": []
+    },
+    {
+      "cell_type": "code",
+      "metadata": {},
+      "source": [
+        "# Method 2:  Initialize empty, add two layers using dict\n",
+        "china = geodata.Mask(\"China\")\n",
+        "china.add_layer(layer_path={\"elevation\": elevation_path, \"modis\": modis_path})"
+      ],
+      "execution_count": null,
+      "outputs": []
+    },
+    {
+      "cell_type": "code",
+      "metadata": {},
+      "source": [
+        "# Method 3:  Initalize with two layers passed as list\n",
+        "china = geodata.Mask(\n",
+        "    \"China\", layer_path=[elevation_path, modis_path], layer_name=[\"elevation\", \"modis\"]\n",
+        ")"
+      ],
+      "execution_count": null,
+      "outputs": []
+    },
+    {
+      "cell_type": "code",
+      "metadata": {},
+      "source": [
+        "# Method 4:  Initialize with two layers passed as dict\n",
+        "china = geodata.Mask(\n",
+        "    \"China\", layer_path={\"elevation\": elevation_path, \"modis\": modis_path}\n",
+        ")"
+      ],
+      "execution_count": null,
+      "outputs": []
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "Display the mask object in the jupyter notebook:"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "metadata": {},
+      "source": [
+        "china"
+      ],
+      "execution_count": null,
+      "outputs": []
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "Each mask object has several attributes:\n",
+        "- `layers`: a dictionary of name (key)  - rasterio file opener (values). The <\\open DatasetReader> can be the input for many other mask methods for the module.  \n",
+        "- `merged_mask`: the merged and flatten mask of its layers, the merged raster from `layers`\n",
+        "- `shape_mask`: similar to the `layers` attribute, but a dictionary of extracted shapes from the merged mask by default. Users may also extracted shape masks from specified layers in `self.layers`.\n",
+        "- `saved`: whether this mask object has been saved locally.\n",
+        "- `mask_dir`: the directory to save the mask object, by default it should be the mask dir in config.py."
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "Show the `slope` layer in mask `china`. The `show` method will always try to show the proper latitude and longitude, unless we call it `show(layer, lat_lon = False)`."
+      ]
+    },
+    {
+      "cell_type": "code",
+      "metadata": {},
+      "source": [
+        "china.layers[\"elevation\"]"
+      ],
+      "execution_count": null,
+      "outputs": []
+    },
+    {
+      "cell_type": "code",
+      "metadata": {},
+      "source": [
+        "show(china.layers[\"elevation\"], title=\"Elevation of China in meters\")"
+      ],
+      "execution_count": null,
+      "outputs": []
+    },
+    {
+      "cell_type": "code",
+      "metadata": {},
+      "source": [
+        "china.layers"
+      ],
+      "execution_count": null,
+      "outputs": []
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "**Some useful methods to examine the layers**\n",
+        "\n",
+        "- `china.get_res()`: get resolution of each layer, in lat-lon coordinates\n",
+        "- `china.get_res(product = True)`: get grid cell size, in product of lat-lon coordinate differences\n",
+        "- `china.get_bounds()`: get bounds, in lat-lon coordinates"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "metadata": {},
+      "source": [
+        "china.get_bounds()"
+      ],
+      "execution_count": null,
+      "outputs": []
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "Note that the modis layer has a very different bounding box then the slope layer in lat-lon coordinate system. This is because the modis layer was converted to the lat-lon CRS from a different CRS when it was added to the object. The following section will explore CRS conversion."
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "### CRS conversion, trimming, and cropping (Optional)"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "Method `open_tif` can open a layer without adding it to the layer, this allows us to visualize it before-hand. It is a good practice to close the raster after opening it to avoid writing permission conflict issues. Closing the raster below does not involve any layer operation associated with the mask object. "
+      ]
+    },
+    {
+      "cell_type": "code",
+      "metadata": {},
+      "source": [
+        "modis_opener = geodata.mask.open_tif(modis_path, show_raster=True)\n",
+        "modis_opener.close()"
+      ],
+      "execution_count": null,
+      "outputs": []
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "We can use `remove_layer` method to remove a layer to mask `china`. This method will properly close the raster file, because the raster file would remain open after being added to the mask."
+      ]
+    },
+    {
+      "cell_type": "code",
+      "metadata": {},
+      "source": [
+        "china.remove_layer(\"modis\")"
+      ],
+      "execution_count": null,
+      "outputs": []
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "The `add_layer` method incorporates coordinate reference system (CRS) conversion to lat-lon (EPSG:4326), if necessary. Note that this method will overwrite the layer by default, if it is in the object already, unless the user specifies `replace=False`. \n",
+        "\n",
+        "The method will automatically trim the all-zero columns/rows. By default, the paramater `trim` is set to `True`. If we do not set it to True, we might generate a converted raster with new CRS but many all-zero columns and rows."
+      ]
+    },
+    {
+      "cell_type": "code",
+      "metadata": {},
+      "source": [
+        "china.add_layer(modis_path, \"modis\", trim=False)\n",
+        "show(china.layers[\"modis\"], title=\"China Modis CRS converted (No trimming)\")"
+      ],
+      "execution_count": null,
+      "outputs": []
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "We can also crop a raster/layer with user-defined dimensions: method `crop_layer` can take either starting indices of top/left, ending indices of right/bottom, or coordinates values in lat/long to trim the raster.\n",
+        "\n",
+        "The difference between `crop_layer` and `trim_layer` is that `crop_layer` must take in user specified range to crop the raster, and `trim_layer` would remove the all zero rows and columns automatically for a raster. So that if the user do not know which index to start and end to remove the empty rows/columns, `trim_raster` is better.\n",
+        "\n",
+        "The method `crop_raster` (`geodata.mask.crop_raster`) is similar to `crop_layer` but can take a layer name as input, so that the user does not need to add a raster as a layer to call that method. (Similar method: `trim_layer`/`trim_raster`, `binarize_layer`/`binarize_raster`)"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "metadata": {},
+      "source": [
+        "china.crop_layer(\"modis\", bounds=(73, 17, 135, 54))\n",
+        "show(china.layers[\"modis\"], title=\"China Modis Layer Cropped\")"
+      ],
+      "execution_count": null,
+      "outputs": []
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "This performs the same function by passing the layer to `crop_raster`:"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "metadata": {},
+      "source": [
+        "china.layers[\"modis\"] = geodata.mask.crop_raster(\n",
+        "    china.layers[\"modis\"], (73, 17, 135, 54)\n",
+        ")"
+      ],
+      "execution_count": null,
+      "outputs": []
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "### Filter a layer"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "The mask module also supports filtering a layer based on list of categorical values, a minimum (lower) boundary, or maximum (upper) boundary.\n",
+        "\n",
+        "In the `filter_raster` method, a user may specify any of the `value` (the list of numberic values in the raster array to be selected), `max_bound`, and `min_bound` parameters to selected desired values. If the parameter `binarize` is False (by default), the method will return the original values of the raster that satisfy the conditions, otherwise the method will return 1 for the values that satisfy the conditions and 0 elsewhere."
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "#### Select Categorical Values from MODIS Layer\n",
+        "\n",
+        "Since the modis layer has 17 distinct values for different land use types, we want to create a layer of binary values, indicating unavailable land as 0, and available land as 1.\n",
+        "\n",
+        "We wish to create a mask where :\n",
+        "\n",
+        "- all forested areas (values 1-5) are 0 (i.e., unsuitable)\n",
+        "- all urban areas (13) are 0\n",
+        "- all others are 1\n",
+        "\n",
+        "\n",
+        "Let us use method `filter_raster` to create a layer of `modis_filtered` binary mask, where 1, 2, 3, 4, 5, and 13 will be unavailable land assigned 0 and the rest of the values will be 1 (available).\n"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "metadata": {},
+      "source": [
+        "avail_values = list(set(range(1, 18)) - {1, 2, 3, 4, 5, 13})\n",
+        "avail_values"
+      ],
+      "execution_count": null,
+      "outputs": []
+    },
+    {
+      "cell_type": "code",
+      "metadata": {},
+      "source": [
+        "china.layers[\"modis_filtered\"] = geodata.mask.filter_raster(\n",
+        "    china.layers[\"modis\"], binarize=True, values=avail_values\n",
+        ")"
+      ],
+      "execution_count": null,
+      "outputs": []
+    },
+    {
+      "cell_type": "code",
+      "metadata": {
+        "scrolled": false
+      },
+      "source": [
+        "china.remove_layer(\"modis\")\n",
+        "show(china.layers[\"modis_filtered\"])"
+      ],
+      "execution_count": null,
+      "outputs": []
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "#### Filter elevation layer\n",
+        "\n",
+        "Because we cannot build renewable energy in areas with high elevation, let us set the constraint from the `elevation` layer, by using elevation < 4000m at 1 and other areas as 0. The result layer `elevation_filtered` will have only 1 and 0 as unique values."
+      ]
+    },
+    {
+      "cell_type": "code",
+      "metadata": {},
+      "source": [
+        "china.filter_layer(\n",
+        "    \"elevation\", dest_layer_name=\"elevation_filtered\", max_bound=4000, binarize=True\n",
+        ")"
+      ],
+      "execution_count": null,
+      "outputs": []
+    },
+    {
+      "cell_type": "code",
+      "metadata": {},
+      "source": [
+        "china.remove_layer(\"elevation\")\n",
+        "show(china.layers[\"elevation_filtered\"])"
+      ],
+      "execution_count": null,
+      "outputs": []
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "#### Filter Slope Layer\n",
+        "\n",
+        "We also cannot build renewable energy in area with large slopes, so let us set another constraint from the `slope` layer from the slope tif file, by using slope < 20 degree at 1 and else as 0. The result layer `slope_filtered` will have only 1 and 0 as unique values."
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "First, add the slope raster to the china mask."
+      ]
+    },
+    {
+      "cell_type": "code",
+      "metadata": {},
+      "source": [
+        "china.add_layer(slope_path, layer_name=\"slope\")\n",
+        "show(china.layers[\"slope\"])"
+      ],
+      "execution_count": null,
+      "outputs": []
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "Filter the raster, delete the old slope layer."
+      ]
+    },
+    {
+      "cell_type": "code",
+      "metadata": {},
+      "source": [
+        "china.filter_layer(\n",
+        "    \"slope\", dest_layer_name=\"slope_filtered\", max_bound=20, binarize=True\n",
+        ")"
+      ],
+      "execution_count": null,
+      "outputs": []
+    },
+    {
+      "cell_type": "code",
+      "metadata": {
+        "scrolled": false
+      },
+      "source": [
+        "china.remove_layer(\"slope\")\n",
+        "show(china.layers[\"slope_filtered\"])"
+      ],
+      "execution_count": null,
+      "outputs": []
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "### Additional Visualization Options"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "We can plot the provinces on a selected layer by taking `shape` input in the `show()` method. Here, we will use the `china_shapes` that we obtained from `all_shape`. Its `geometry` column is a Series of shapes (shapely.geometry or MultiPolygon) for Chinese provinces."
+      ]
+    },
+    {
+      "cell_type": "code",
+      "metadata": {},
+      "source": [
+        "show(china.layers[\"modis_filtered\"], shape=china_shapes[\"geometry\"])"
+      ],
+      "execution_count": null,
+      "outputs": []
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "## Adding Shape Features as a Layer"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "Recall that we have previously loaded the environmental protected shapes of China in a GeoPandas dataframe."
+      ]
+    },
+    {
+      "cell_type": "code",
+      "metadata": {},
+      "source": [
+        "len(wdpa_shapes)"
+      ],
+      "execution_count": null,
+      "outputs": []
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "The three shapefiles have 78 features altogether, but we want to add all the features to one new layer instead of 78 new layers. The input shape should be a python dictionary, where there is a key for each unique shape. Also, in the `add_shape_layer` method, we will specify a `combine_name` to combine the features into one layer in this case, since we want the mask to have just one more layers, not 78 more layers.\n",
+        "\n",
+        "When adding a shapefile, we must specify the dimensions. We will also use `reference layer = 'slope_filtered'` so the new shape layer will have the same dimension with the `slope_filtered` layer. If the mask is empty and does not contain any layer, the user will have to specify the `resolution` parameter for the raster layer dimension.\n",
+        "\n",
+        "By default, this method will have paramater `exclude` that defaults to `False`. When it is true, area inside the shape is 0. When it is false, area inside the shape is 1. In this use case, however, we want 0 for area inside of the shape as they are environmental protected areas to exclude. We can just use the default method call."
+      ]
+    },
+    {
+      "cell_type": "code",
+      "metadata": {},
+      "source": [
+        "china.add_shape_layer(\n",
+        "    wdpa_shapes[\"geometry\"].to_dict(),\n",
+        "    reference_layer=\"slope_filtered\",\n",
+        "    combine_name=\"protected\",\n",
+        ")\n",
+        "show(\n",
+        "    china.layers[\"protected\"],\n",
+        "    title=\"WDPA Protected area shape features as a new layer\",\n",
+        "    grid=True,\n",
+        ")"
+      ],
+      "execution_count": null,
+      "outputs": []
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "We can also use the parameter `buffer` in `add_shape_layer` method to create an approximate representation of all locations within a given (perpindicular) distance of the shape object. The units for the buffer are given in kilometers.\n",
+        "\n",
+        "Note that since the units of the original shape are in lat-lon coordinates, when we add the buffer, we will need to have a CRS that has meter as unit. The program will convert the shapes to that CRS, add the buffer around shapes, then convert it back to the lat-lon CRS system. By default, we used \"EPSG:6933\", an equal area projection CRS to add buffer in kilometer. "
+      ]
+    },
+    {
+      "cell_type": "code",
+      "metadata": {},
+      "source": [
+        "km_buffer = 20\n",
+        "\n",
+        "china.add_shape_layer(\n",
+        "    wdpa_shapes[\"geometry\"].to_dict(),\n",
+        "    reference_layer=\"slope_filtered\",\n",
+        "    combine_name=\"protected_with_buffer\",\n",
+        "    buffer=km_buffer,\n",
+        ")\n",
+        "\n",
+        "show(\n",
+        "    china.layers[\"protected_with_buffer\"],\n",
+        "    title=f\"WDPA Protected area shape with {km_buffer}km buffer\",\n",
+        "    grid=True,\n",
+        ")\n",
+        "\n",
+        "china.remove_layer(\"protected_with_buffer\")"
+      ],
+      "execution_count": null,
+      "outputs": []
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "## Merging and Flattening Layers"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "In order to combine all layers into one, we use the `merge_layer` method which creates a new layer called `merged_mask`. This merges multiple layers together and flattens them using either **and** (default) or **sum** method, saving the result to `self.merged_mask` by default. Geospatial bounds and resolution of the output layer are in the units of the input file coordinate reference system, but by default, we will use the resolution of the layer with the best (finest) resolution for the output bounds/resolution, unless a reference layer is provided. In this case, the resolution of the merged_mask is the same with the `modis_filtered` layer. "
+      ]
+    },
+    {
+      "cell_type": "code",
+      "metadata": {},
+      "source": [
+        "china.get_res()"
+      ],
+      "execution_count": null,
+      "outputs": []
+    },
+    {
+      "cell_type": "code",
+      "metadata": {},
+      "source": [
+        "china.merge_layer(attribute_save=False, show_raster=False).res"
+      ],
+      "execution_count": null,
+      "outputs": []
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "### Binary `AND` Method"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "By default, the `merge_layer` method will use a binary 'and' method: for each grid cell, if any of the n layers are 0, then the returned `self.merged_layer` will also have 0 at that location. In other words, if all the layers indicate that a land is available (!=0), the merged result will have value 1.\n",
+        "\n",
+        "`merge_layer` may also take in an optional parameter `layers`, which is a list of layer names stored in the object, if the user does not wish to merge all layers in the object. If the user does not want to save the result to the `merged_mask` attribute, the user can specify `attribute_save = False`."
+      ]
+    },
+    {
+      "cell_type": "code",
+      "metadata": {},
+      "source": [
+        "# merge and plot only, do not save\n",
+        "china.merge_layer(attribute_save=False, layers=[\"slope_filtered\", \"modis_filtered\"])"
+      ],
+      "execution_count": null,
+      "outputs": []
+    },
+    {
+      "cell_type": "code",
+      "metadata": {},
+      "source": [
+        "china"
+      ],
+      "execution_count": null,
+      "outputs": []
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "Try again with the `reference_layer` parameter:"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "metadata": {},
+      "source": [
+        "china.merge_layer(\n",
+        "    layers=[\"elevation_filtered\", \"modis_filtered\"],\n",
+        "    reference_layer=\"elevation_filtered\",\n",
+        "    show_raster=False,\n",
+        ")"
+      ],
+      "execution_count": null,
+      "outputs": []
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "The result of the `merged_mask` method is saved to `china.merged_mask` with the same resolution as the reference layer, in this case `elevation_filtered`."
+      ]
+    },
+    {
+      "cell_type": "code",
+      "metadata": {},
+      "source": [
+        "china.merged_mask.res"
+      ],
+      "execution_count": null,
+      "outputs": []
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "For the purpose of this demonstration, we will select the `AND` method for the final merged_mask. We can also trim the border of the merged mask since the 4 layers have different boundaries. We can set the parameter `trim = True`."
+      ]
+    },
+    {
+      "cell_type": "code",
+      "metadata": {},
+      "source": [
+        "china.merge_layer(trim=True)"
+      ],
+      "execution_count": null,
+      "outputs": []
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "### `SUM` Method"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "The sum method will add up the values from all the layers using weights. When there is no weight dict provided, all the layers for merging will have weights of 1 by default.\n",
+        "\n",
+        "Note: since we are not using the sum method to proceed to the following sections, we will keep `attribute_save = False` to prevent this method from overwriting the mask we have previously created above."
+      ]
+    },
+    {
+      "cell_type": "code",
+      "metadata": {},
+      "source": [
+        "china.merge_layer(method=\"sum\", attribute_save=False, trim=True)"
+      ],
+      "execution_count": null,
+      "outputs": []
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "This distribution is completely arbitrary for the purpose of demonstration of the module: (Note: The weights do not need to have a total of 1)\n",
+        "- elevation_filtered: 0.15, slope_filtered: 0.1, modis_filtered: 0.3, protected: 0.45\n",
+        "\n",
+        "We will write the result to a new variable `customized_merged_layer` for continuing processing."
+      ]
+    },
+    {
+      "cell_type": "code",
+      "metadata": {},
+      "source": [
+        "customized_merged_layer = china.merge_layer(\n",
+        "    method=\"sum\",\n",
+        "    weights={\n",
+        "        \"elevation_filtered\": 0.15,\n",
+        "        \"slope_filtered\": 0.1,\n",
+        "        \"modis_filtered\": 0.3,\n",
+        "        \"protected\": 0.45,\n",
+        "    },\n",
+        "    attribute_save=False,\n",
+        "    trim=True,\n",
+        ")"
+      ],
+      "execution_count": null,
+      "outputs": []
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "If the continuous value created by `merged_mask` represents a suitability metric, we could set a minimum value of 0.8 to be considered \"suitable\" (or 1). We then apply the `filter_raster` method on the merged layer."
+      ]
+    },
+    {
+      "cell_type": "code",
+      "metadata": {},
+      "source": [
+        "customized_merged_layer = geodata.mask.filter_raster(\n",
+        "    customized_merged_layer, min_bound=0.8, binarize=True\n",
+        ")\n",
+        "show(customized_merged_layer)"
+      ],
+      "execution_count": null,
+      "outputs": []
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "## Eliminate Small Contiguous Areas"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "Using the above methods, we might end up with many small contiguous areas that are marked suitable but surrounded by an unsuitable region. We may want to exclude such regions from renewable energy development. The `filter_area` method will remove the small contiguous suitable regions by transforming the merged mask raster to polygons/shapes, calculating the area of each polygon, and filtering out polygons that are smaller than a given threshold. Units are given in kilometer-squared (km$^2$).\n",
+        "\n",
+        "By default, `filter_area` uses the merged mask raster and returns a new raster, unless input/output layers are specified by `layer_name` and `dest_layer_name`. \n",
+        "\n",
+        "By default, its `shape_value` parameter is 1, indicating that we are only interested in finding all groups of cells with value 1 (suitable) for elimination. We specify the threshold with the `min_area` parameter.\n",
+        "\n",
+        "Note: the `filter_area` method may take a long time (5 or more minutes depending on the complexity of your layer and your computational setup). The method relies upon `rasterio.rasterize`, see performance notes: https://rasterio.readthedocs.io/en/latest/api/rasterio.features.html#rasterio.features.rasterize\n"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "For example, if we focus on Guangdong province in Southern China from the merged mask, we notice that there are many small islands in the ocean that are marked as suitable areas. We want to exclude these small regions from our merged mask."
+      ]
+    },
+    {
+      "cell_type": "code",
+      "metadata": {},
+      "source": [
+        "plt.imshow(china.merged_mask.read(1)[4800:5300, 5700:6600], interpolation=\"none\")\n",
+        "plt.show()"
+      ],
+      "execution_count": null,
+      "outputs": []
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "Call `filter_area` to remove all contiguous suitable region shapes smaller than 100 km$^2$:"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "metadata": {},
+      "source": [
+        "china.merged_mask = geodata.mask.filter_area(china, min_area=100)"
+      ],
+      "execution_count": null,
+      "outputs": []
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "There shapes are removed in the new merged_mask."
+      ]
+    },
+    {
+      "cell_type": "code",
+      "metadata": {},
+      "source": [
+        "plt.imshow(china.merged_mask.read(1)[4800:5300, 5700:6600], interpolation=\"none\")\n",
+        "plt.show()"
+      ],
+      "execution_count": null,
+      "outputs": []
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "## Extracting Shapes from Masks"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "Sometimes the user needs to generate masks and perform analysis for a collection of regions (e.g., at the state/province level). The purpose of shape extraction (`extract_shapes`) is to separate `merged_mask` values for each region, with the result a dictionary of name-mask pairs in the `shape_mask` attribute of the mask object. The values of `shape_mask` will be 0 outside of the shape, and will be `merged_mask` inside of the shape."
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "For the purpose of this demonstration, we will only select the province of Jiangsu, Zhejiang, and Shanghai."
+      ]
+    },
+    {
+      "cell_type": "code",
+      "metadata": {
+        "scrolled": false
+      },
+      "source": [
+        "china_shapes_subset = china_shapes[\n",
+        "    china_shapes[\"name\"].isin([\"Jiangsu\", \"Zhejiang\", \"Shanghai\"])\n",
+        "]\n",
+        "china_shapes_subset"
+      ],
+      "execution_count": null,
+      "outputs": []
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "Converting the filtered shape dictionary to a python dictionary as the input for `extract_shapes`, where the keys for the dictionary will be the names of the new extracted shape layers."
+      ]
+    },
+    {
+      "cell_type": "code",
+      "metadata": {},
+      "source": [
+        "china_shapes_subset = (\n",
+        "    china_shapes_subset[[\"name\", \"geometry\"]].set_index(\"name\")[\"geometry\"].to_dict()\n",
+        ")\n",
+        "china_shapes_subset"
+      ],
+      "execution_count": null,
+      "outputs": []
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "Extract the shapes from the merged_mask. "
+      ]
+    },
+    {
+      "cell_type": "code",
+      "metadata": {},
+      "source": [
+        "china.extract_shapes(china_shapes_subset)"
+      ],
+      "execution_count": null,
+      "outputs": []
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "The resulting mask object contains the dictionary `shape_mask` with the extracted values:"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "metadata": {},
+      "source": [
+        "china"
+      ],
+      "execution_count": null,
+      "outputs": []
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "## Saving and Loading Masks"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "metadata": {},
+      "source": [
+        "china.save_mask()"
+      ],
+      "execution_count": null,
+      "outputs": []
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "With the mask saved, the user can now load the layers or shapes with `xarray` instead if preferred."
+      ]
+    },
+    {
+      "cell_type": "code",
+      "metadata": {},
+      "source": [
+        "shape_xr_lst = china.load_shape_xr()\n",
+        "shape_xr_lst[\"Zhejiang\"].plot()"
+      ],
+      "execution_count": null,
+      "outputs": []
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "Optional: closing all the files when saving the mask. This can avoid possible write permission error."
+      ]
+    },
+    {
+      "cell_type": "code",
+      "metadata": {},
+      "source": [
+        "china.save_mask(close_files=True)"
+      ],
+      "execution_count": null,
+      "outputs": []
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "Loading a previously saved mask."
+      ]
+    },
+    {
+      "cell_type": "code",
+      "metadata": {},
+      "source": [
+        "china_2 = geodata.mask.load_mask(\"china\")"
+      ],
+      "execution_count": null,
+      "outputs": []
+    },
+    {
+      "cell_type": "code",
+      "metadata": {},
+      "source": [
+        "china_2"
+      ],
+      "execution_count": null,
+      "outputs": []
+    }
+  ],
+  "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.10.11"
+    }
+  },
+  "nbformat": 4,
+  "nbformat_minor": 4
+}
\ No newline at end of file
diff --git a/docs/source/mask/merge_layer_known_issues.md b/docs/source/mask/merge_layer_known_issues.md
new file mode 100644
index 00000000..f5502895
--- /dev/null
+++ b/docs/source/mask/merge_layer_known_issues.md
@@ -0,0 +1,34 @@
+# `merge_layer` known issues (historical)
+
+```{note}
+**Historical context.** Older geodata versions could raise ``RasterioIOError: No such
+file or directory`` when merging **in-memory** (``/vsimem``) mask layers. Current code
+pins memory files for the lifetime of each layer reader so ``filter_layer`` →
+``merge_layer`` normally works without pre-saving layers to disk.
+```
+
+## Symptom
+
+``merge_layer`` (or ``merge_and`` / ``merge_sum`` after filters) fails with an error
+referring to a missing path under ``/vsimem/``.
+
+## Cause (legacy behavior)
+
+Raster layers stored in GDAL memory files were sometimes closed before merge read them
+back, so the virtual path was no longer valid.
+
+## Current behavior
+
+The mask module keeps layer readers alive while a ``Mask`` object uses in-memory
+layers. If you still see this error on an old install, upgrade geodata or save
+intermediate layers to disk before merging.
+
+## Workaround (older versions)
+
+1. Save filtered layers to GeoTIFF before ``merge_layer``.
+2. Call ``save_mask(close_files=True)`` when finished, and avoid two ``Mask`` objects
+   opening the same files simultaneously (see [mask troubleshooting](mask_troubleshoot.md)).
+
+## Tests
+
+Regression coverage lives under ``tests/pr/mask/test_mask_merge_inmemory.py``.
diff --git a/docs/source/mask/xarray_mask_tutorial.ipynb b/docs/source/mask/xarray_mask_tutorial.ipynb
new file mode 100644
index 00000000..fdc86f89
--- /dev/null
+++ b/docs/source/mask/xarray_mask_tutorial.ipynb
@@ -0,0 +1,312 @@
+{
+  "cells": [
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "# Tutorial: Applying Saved Masks with `XarrayMask`\n",
+        "\n",
+        "This notebook shows how to apply a **saved** geographic mask to model or analysis\n",
+        "output represented as an `xarray.Dataset` or `xarray.DataArray` — without using\n",
+        "`Cutout.add_mask` or `Cutout.mask`.\n",
+        "\n",
+        "For contributor notes on the xarray masking design, see\n",
+        "[development/xarray_mask_workflow](../development/xarray_mask_workflow.rst).\n",
+        "To build masks from rasters and shapefiles, see\n",
+        "[mask creation workflow](mask_creation_workflow.ipynb)."
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "## Overview\n",
+        "\n",
+        "| Step | API | Module |\n",
+        "|------|-----|--------|\n",
+        "| Create and save a mask | `geodata.Mask` | `src/geodata/mask.py` |\n",
+        "| Run a model (wind, pvlib, …) | model `estimate()` | `src/geodata/model/` |\n",
+        "| Align mask to your grid, attach or apply | `geodata.XarrayMask` | `src/geodata/mask/xarray_mask.py` |\n",
+        "\n",
+        "**`XarrayMask` does not replace mask creation.** It loads a saved mask and applies it\n",
+        "to xarray data on your target grid."
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "## Setup\n",
+        "\n",
+        "This tutorial runs **offline** using a small synthetic grid and a temporary mask\n",
+        "directory. The same API calls work for production masks saved under `GEODATA_ROOT`."
+      ]
+    },
+    {
+      "cell_type": "code",
+      "metadata": {},
+      "source": [
+        "import tempfile\n",
+        "from pathlib import Path\n",
+        "\n",
+        "import numpy as np\n",
+        "import rasterio as ras\n",
+        "import shapely.geometry\n",
+        "import xarray as xr\n",
+        "from rasterio.transform import from_bounds\n",
+        "\n",
+        "from geodata import Mask, XarrayMask"
+      ],
+      "execution_count": null,
+      "outputs": []
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "## Step 1: Stand in for model output\n",
+        "\n",
+        "Your analysis dataset can use `x`/`y` or `lat`/`lon`. `XarrayMask` normalizes\n",
+        "coordinates via `ds_reformat_index` before alignment.\n",
+        "\n",
+        "Below we use a small `(time, y, x)` dataset as if it came from a wind or PV model."
+      ]
+    },
+    {
+      "cell_type": "code",
+      "metadata": {},
+      "source": [
+        "y = np.array([30.75, 30.5, 30.25, 30.0])\n",
+        "x = np.array([100.0, 100.25, 100.5, 100.75])\n",
+        "time = np.array([\"2016-01-01T00:00:00\", \"2016-01-01T01:00:00\"], dtype=\"datetime64[ns]\")\n",
+        "\n",
+        "values = np.arange(len(time) * len(y) * len(x), dtype=np.float32).reshape(\n",
+        "    len(time), len(y), len(x)\n",
+        ")\n",
+        "model_ds = xr.Dataset(\n",
+        "    {\"signal\": ((\"time\", \"y\", \"x\"), values)},\n",
+        "    coords={\"time\": time, \"y\": y, \"x\": x},\n",
+        ")\n",
+        "model_ds"
+      ],
+      "execution_count": null,
+      "outputs": []
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "## Step 2: Create and save a mask (offline example)\n",
+        "\n",
+        "In practice you build masks with `Mask.add_layer`, `filter_layer`, `merge_layer`,\n",
+        "and `save_mask()` — see [mask creation workflow](mask_creation_workflow.ipynb).\n",
+        "\n",
+        "Mask rasters are often stored at **higher resolution** than model output.\n",
+        "`XarrayMask` coarsens them onto `grid` automatically.\n",
+        "\n",
+        "The helper below mirrors `tests/pr/mask/test_xarray_mask.py`."
+      ]
+    },
+    {
+      "cell_type": "code",
+      "metadata": {},
+      "source": [
+        "mask_dir = Path(tempfile.mkdtemp(prefix=\"geodata_xmask_tutorial_\"))\n",
+        "mask_name = \"tutorial_mask\"\n",
+        "\n",
+        "lon_step = float(np.abs(x[1] - x[0]))\n",
+        "lat_step = float(np.abs(y[1] - y[0]))\n",
+        "west = float(x.min() - lon_step / 2)\n",
+        "east = float(x.max() + lon_step / 2)\n",
+        "south = float(y.min() - lat_step / 2)\n",
+        "north = float(y.max() + lat_step / 2)\n",
+        "\n",
+        "nlon_hi = len(x) * 2\n",
+        "nlat_hi = len(y) * 2\n",
+        "transform = from_bounds(west, south, east, north, nlon_hi, nlat_hi)\n",
+        "\n",
+        "arr = np.zeros((nlat_hi, nlon_hi), dtype=np.uint8)\n",
+        "arr[nlat_hi // 4 : 3 * nlat_hi // 4, nlon_hi // 4 : 3 * nlat_hi // 4] = 1\n",
+        "\n",
+        "layer_path = mask_dir / \"source.tif\"\n",
+        "with ras.open(\n",
+        "    str(layer_path),\n",
+        "    \"w\",\n",
+        "    driver=\"GTiff\",\n",
+        "    height=arr.shape[0],\n",
+        "    width=arr.shape[1],\n",
+        "    count=1,\n",
+        "    dtype=arr.dtype,\n",
+        "    compress=\"lzw\",\n",
+        "    crs=\"+proj=latlong\",\n",
+        "    transform=transform,\n",
+        ") as dst:\n",
+        "    dst.write(arr, 1)\n",
+        "\n",
+        "mask = Mask(name=mask_name, mask_dir=str(mask_dir))\n",
+        "mask.add_layer(str(layer_path), layer_name=\"source\")\n",
+        "mask.merge_layer(show_raster=False)\n",
+        "\n",
+        "region = shapely.geometry.box(west, south, (west + east) / 2, (south + north) / 2)\n",
+        "mask.extract_shapes({\"region_a\": region}, show_raster=False)\n",
+        "mask.save_mask()\n",
+        "\n",
+        "print(f\"Saved mask '{mask_name}' under {mask_dir}\")"
+      ],
+      "execution_count": null,
+      "outputs": []
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "## Step 3: Load and align — `XarrayMask.from_name`\n",
+        "\n",
+        "Pass your model grid so the saved mask is coarsened and aligned to the same\n",
+        "`x`/`y` (or `lat`/`lon`) coordinates."
+      ]
+    },
+    {
+      "cell_type": "code",
+      "metadata": {},
+      "source": [
+        "xmask = XarrayMask.from_name(mask_name, grid=model_ds, mask_dir=str(mask_dir))\n",
+        "xmask"
+      ],
+      "execution_count": null,
+      "outputs": []
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "You can also build from an in-memory `Mask` object:\n",
+        "\n",
+        "```python\n",
+        "loaded = Mask.from_name(mask_name, mask_dir=str(mask_dir))\n",
+        "xmask = XarrayMask.from_mask(loaded, grid=model_ds)\n",
+        "```"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "## Step 4: Attach — legacy-compatible output\n",
+        "\n",
+        "`attach()` returns a dict of datasets (keys: `merged_mask`, plus any shape masks).\n",
+        "Each dataset contains your original variables plus `mask` and optional `area` — the\n",
+        "same structure as `Cutout.mask()`."
+      ]
+    },
+    {
+      "cell_type": "code",
+      "metadata": {},
+      "source": [
+        "attached = xmask.attach(model_ds, include_area=True)\n",
+        "list(attached.keys())"
+      ],
+      "execution_count": null,
+      "outputs": []
+    },
+    {
+      "cell_type": "code",
+      "metadata": {},
+      "source": [
+        "merged = attached[\"merged_mask\"]\n",
+        "merged"
+      ],
+      "execution_count": null,
+      "outputs": []
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "## Step 5: Apply — filtered outputs\n",
+        "\n",
+        "- `mode=\"where\"` — set values outside the mask to NaN\n",
+        "- `mode=\"multiply\"` — set values outside the mask to zero"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "metadata": {},
+      "source": [
+        "where_out = xmask.apply(model_ds, mode=\"where\", include_area=True)[\"merged_mask\"]\n",
+        "multiply_out = xmask.apply(model_ds, mode=\"multiply\", include_area=False)[\"merged_mask\"]\n",
+        "\n",
+        "where_out[\"signal\"].isel(time=0)"
+      ],
+      "execution_count": null,
+      "outputs": []
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "## Step 6: Area-weighted aggregation\n",
+        "\n",
+        "With `attach(..., include_area=True)` you can compute mask- and area-weighted\n",
+        "statistics over time — the same pattern as the legacy Cutout workflow."
+      ]
+    },
+    {
+      "cell_type": "code",
+      "metadata": {},
+      "source": [
+        "ds = attached[\"merged_mask\"]\n",
+        "weighted_mean = (\n",
+        "    (ds[\"signal\"] * ds[\"mask\"] * ds[\"area\"]).sum(dim=[\"lat\", \"lon\"])\n",
+        "    / (ds[\"mask\"] * ds[\"area\"]).sum(dim=[\"lat\", \"lon\"])\n",
+        ")\n",
+        "weighted_mean"
+      ],
+      "execution_count": null,
+      "outputs": []
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "## Production usage\n",
+        "\n",
+        "When your mask is already saved under the default mask directory (`GEODATA_ROOT`):\n",
+        "\n",
+        "```python\n",
+        "xmask = XarrayMask.from_name(\"china\", grid=output_ds)  # uses geodata.config.MASK_DIR\n",
+        "masked = xmask.apply(output_ds, mode=\"where\")\n",
+        "```\n",
+        "\n",
+        "### Typical pipeline\n",
+        "\n",
+        "1. `output_ds = model.estimate(...)`\n",
+        "2. `xmask = XarrayMask.from_name(\"my_mask\", grid=output_ds, mask_dir=...)`\n",
+        "3. `xmask.attach(output_ds)` or `xmask.apply(output_ds, ...)`\n",
+        "\n",
+        "### See also\n",
+        "\n",
+        "| Topic | Page |\n",
+        "|-------|------|\n",
+        "| Create masks from GIS layers | [mask_creation_workflow](mask_creation_workflow.ipynb) |\n",
+        "| Legacy Cutout masking | [mask_on_cutout](../legacy/mask_on_cutout.ipynb) |\n",
+        "| Xarray masking design notes (contributors) | [xarray_mask_workflow](../development/xarray_mask_workflow.rst) |\n",
+        "| Migration plan (contributors) | [mask_xarray_migration_plan](../development/mask_xarray_migration_plan.md) |\n",
+        "| Automated examples | `tests/pr/mask/test_xarray_mask.py`, `tests/pr/test_wind_xarraymask_integration.py` |"
+      ]
+    }
+  ],
+  "metadata": {
+    "kernelspec": {
+      "display_name": "Python 3",
+      "language": "python",
+      "name": "python3"
+    },
+    "language_info": {
+      "name": "python",
+      "version": "3.11.0"
+    }
+  },
+  "nbformat": 4,
+  "nbformat_minor": 5
+}
\ No newline at end of file
diff --git a/docs/source/modeling/era5_outputs.md b/docs/source/modeling/era5_outputs.md
new file mode 100644
index 00000000..834042c2
--- /dev/null
+++ b/docs/source/modeling/era5_outputs.md
@@ -0,0 +1,58 @@
+# ERA5 model outputs
+
+After you download ERA5 data (see :ref:`downloading-era5-data` in the
+[Dataset module overview](../datasets/overview.rst)), **models** turn raw files into
+analysis-ready time series. Masking applies afterward on model results (see
+[Mask tutorials](../mask/xarray_mask_tutorial.ipynb)).
+
+This page is a short catalog of common **model** outputs on the current tested path. It
+does not describe legacy ``Cutout`` / ``geodata.convert`` products (see
+[Legacy MERRA2 outputs](../legacy/merra2/merra2_outputs.md)).
+
+## Wind generation time-series
+
+Hub-height **capacity factor** (``cf``) from ERA5 3D wind:
+
+| Step | Component |
+|------|-----------|
+| Dataset | ``wind_3d_hourly`` (or ``wind_3d_hourly_test`` for offline fixtures) |
+| Model | ``WindInterpolationModel`` — [Wind modeling](wind/index.rst), [interpolation tutorial](wind/interpolation.rst) |
+
+```python
+from geodata.datasets import load_dataset
+from geodata.model.wind import WindInterpolationModel
+
+ds = load_dataset("wind_3d_hourly")(years=slice(2016, 2016), months=slice(1, 1))
+if not ds.downloaded:
+    ds.download()
+
+model = WindInterpolationModel(ds)
+model.prepare()
+cf = model.estimate(turbine="Vestas_V112_3MW", years=slice(2016, 2016), months=slice(1, 1))
+```
+
+## Wind speed time-series
+
+Same dataset and model; pass ``height=<meters>`` instead of ``turbine=``:
+
+```python
+wind_speed = model.estimate(height=100.0, years=slice(2016, 2016), months=slice(1, 1))
+```
+
+## Solar photovoltaic generation time-series
+
+Hourly **AC power** (``ac``) and **capacity factor** (``pv``):
+
+| Step | Component |
+|------|-----------|
+| Dataset | ``wind_solar_hourly`` (or ``wind_solar_hourly_test`` for offline fixtures) |
+| Model | ``Pvlib`` — [PVLib modeling](pvlib/index.rst) |
+
+After ``init_pv_system()`` and ``init_model_config()``, call ``estimate()`` (see the
+pvlib docs for ``compact_output`` and spatial subsetting).
+
+## See also
+
+- [ERA5 CDS setup](../datasets/era5.rst)
+- [Offline ERA5 fixtures](../development/offline-era5-fixture-datasets.md)
+- [Supported input/output formats](../quick_start/input_output.md)
diff --git a/docs/source/modeling/pvlib/index.rst b/docs/source/modeling/pvlib/index.rst
index 211cc0ed..68f4a41d 100644
--- a/docs/source/modeling/pvlib/index.rst
+++ b/docs/source/modeling/pvlib/index.rst
@@ -97,12 +97,76 @@ Next, we can estimate the AC Power and PV capacity using the model.
     cf = model.estimate(
         years = slice(2016, 2016), 
         months = slice(1, 1),
-        xs = slice(8, 10), # Optional: specify the bounding box
-        ys = slice(48, 46), # here is an example bounding box for central europe
+        xs = slice(8, 10), # Optional: longitude subset
+        ys = slice(48, 46), # Optional: latitude subset (see below)
     )
     print(cf)
 
-The output will be an xarray Dataset containing the estimated AC Power and PV capacity values for the specified region and time period.
+The output will be an xarray Dataset containing the estimated AC power (``ac``) and
+capacity factor (``pv``) for the specified region and time period.
+
+Estimate options
+----------------
+
+All models inherit a common pattern for **time** and **space** subsetting via
+``estimate()``. The PVLib model adds one extra output option.
+
+Temporal subsetting
+~~~~~~~~~~~~~~~~~~~
+
+Pass ``years`` and ``months`` as ``slice`` objects to limit the period processed.
+Omit either argument to use the prepared model's full range (subject to what was
+available when ``prepare()`` ran).
+
+Spatial subsetting (``xs``, ``ys``)
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+Pass ``xs`` and ``ys`` as ``slice(start, stop)`` to restrict longitude (``x``) and
+latitude (``y``). Omit either argument to keep the full horizontal extent of the
+prepared dataset.
+
+Geodata **normalizes slice bounds** before calling xarray's ``.sel()``. You can pass
+bounds in either order (for example ``ys=slice(48, 46)`` for a band in central
+Europe) and still get a non-empty selection. This matters for ERA5-style grids where
+latitude is often stored in **descending** order: a naive ``slice(46, 48)`` would
+return no points without normalization.
+
+.. code:: Python
+
+    # Equivalent selections on a descending-latitude grid:
+    cf_a = model.estimate(years=slice(2016, 2016), months=slice(1, 1), ys=slice(48, 46))
+    cf_b = model.estimate(years=slice(2016, 2016), months=slice(1, 1), ys=slice(46, 48))
+
+Compact output (``compact_output``)
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+By default, ``estimate()`` returns a compact dataset with only two data variables:
+
+- ``ac`` — AC power (W)
+- ``pv`` — capacity factor (AC output normalized by module nameplate)
+
+Set ``compact_output=False`` to retain **all intermediate weather and ModelChain
+columns** per grid cell (irradiance components, temperature, wind, and other inputs
+used along the chain). Use this for debugging or when you need columns beyond
+``ac`` and ``pv``; the result is larger and slower to write.
+
+.. code:: Python
+
+    # Default: only ac and pv
+    cf = model.estimate(
+        years=slice(2016, 2016),
+        months=slice(1, 1),
+        compact_output=True,
+    )
+    list(cf.data_vars)  # ['ac', 'pv']
+
+    # Full per-coordinate table (debugging / downstream analysis)
+    full = model.estimate(
+        years=slice(2016, 2016),
+        months=slice(1, 1),
+        compact_output=False,
+    )
+    list(full.data_vars)  # ac, pv, plus weather and intermediate columns
 
 .. toctree::
    :maxdepth: 1
diff --git a/docs/source/modeling/wind/extrapolation.rst b/docs/source/modeling/wind/extrapolation.rst
deleted file mode 100644
index a983806e..00000000
--- a/docs/source/modeling/wind/extrapolation.rst
+++ /dev/null
@@ -1,166 +0,0 @@
-Tutorial: Estimate Wind Speed with Extrapolation
-================================================
-
-In this tutorial, we will learn how to estimate wind speed using the extrapolation model
- from the geodata library.
-
-.. warning::
-   Performing wind speed estimation using extrapolation requires a dataset with known
-   wind speed values at **multiple** locations.
-
-   Currently, only the :code:`weather_data_config` :code:`slv_flux_hourly` from the MERRA2 dataset
-   contains the necessary wind speed data for extrapolation.
-
-   Therefore, all of the information below only applies with :code:`slv_flux_hourly` or cutouts
-   derived from it. Using any other dataset will lead to a :code:`ValueError`.
-
-Step 1: Import the necessary libraries
-----------------------------------------
-
-To get started, we need to import the required libraries. We will import the `WindExtrapolationModel` from the `geodata` library, as well as any other libraries needed for data handling and visualization.
-
-.. code:: Python
-
-    import xarray as xr
-
-    from geodata.datasets import load_dataset
-    from geodata.model.wind import WindExtrapolationModel
-
-
-Step 2: Load the dataset
-------------------------
-
-Next, we need to load the dataset that contains the wind speed data. We will use the `slv_flux_hourly` dataset from the ERA5 dataset.
-
-.. code:: Python
-
-    # Load the dataset
-    ds_cls = load_dataset("slv_flux_hourly")
-    ds = ds_cls(
-        years=slice(2006, 2006),
-        months=slice(1, 1),
-        bounds=[-10, 35, 10, 45]  # Optional: specify the bounding box
-    )
-
-    if not ds.downloaded:
-        ds.download()  # Download the data if we don't have it locally
-
-    print(ds.downloaded)  # Check if the dataset is downloaded. Should return True.
-
-
-Step 3: Compute extrapolation parameters
---------------------------------------------
-The extrapolation is separated into two steps, first estimating extrapolation parameters
-using linear regression, and second extrapolating to desired heights.
-First, we compute the extrapolation parameters.
-For more information on the model, see the section below: `How the Extrapolation Model Works`_.
-
-.. code:: Python
-
-    # Create a model based on the above dataset. The model will be associated with
-    # the dataset forever. If you wish to use a different dataset, you will need to
-    # create a new model.
-
-    model = WindExtrapolationModel(ds)
-    model.prepare()
-
-If you have already prepared a cutout with the config :code:`slv_flux_hourly`, you
-can also pass
-that into the model as well. The model treats dataset and cutouts indifferently.
-Simply replace :code:`ds` with your cutout variable.
-
-.. note::
-   The `prepare` method computes the necessary parameters for the extrapolation model
-   based on the loaded dataset. Everything will be saved under the :code:`models`
-   directory under the path :code:`GEODATA_ROOT`.
-
-.. note::
-    It is not necessary to call the `prepare` method every time you want to perform
-    extrapolation. You only need to call it once after loading the dataset. From that
-    point on, you can load and use the model directly without re-preparing it.
-
-Step 4: Estimate using the extrapolation model
-----------------------------------------------
-
-Now that we have prepared the model, we can perform the extrapolation to estimate wind
-speed at the desired locations. Suppose we want to estimate the wind speed at a height
-of 60 above ground during January of 2006 for the entire region covered by the original
-dataset, we can do this as follows:
-
-.. code:: Python
-
-    estimated_wind_speed = model.estimate(
-        height=60,
-        years=slice(2006, 2006),
-        months=slice(1, 1),
-    )
-
-This will return an xarray DataArray containing the estimated wind speed values. Note
-that you can also select a subset area by passing in :code:`xs=slice(start, end)`
-and/or :code:`ys=slice(start, end)` parameters to the `estimate` method.
-
-
-Step 5: Estimate Wind Turbine Capacity Factor (CF) using the interpolation model
---------------------------------------------------------------------------------
-
-Geodata also supports a limited set of wind turbine models to estimate the capacity
-factor (CF) of a wind turbine directly. To get a list of available wind turbine models,
-you can use the `get_available_windturbines` function:
-
-.. code:: Python
-
-    from geodata.resource import get_available_windturbines
-
-    turbines = get_available_windturbines()
-    print(turbines)  # List of available wind turbine configurations
-
-
-To estimate the capacity factor of a wind turbine, you can use the `estimate` method
-and passign in the `turbine` parameter with the name of the wind turbine model.
-
-.. code:: Python
-
-    # Estimate the capacity factor for a specific wind turbine model
-    estimated_cf: xr.Dataset = model.estimate(
-        turbine="Vestas_V112_3MW",  # Example wind turbine model
-        years=slice(2006, 2006),
-        months=slice(1, 1),
-    )
-
-    print(estimated_cf)  # Display the estimated capacity factor
-
-
-The output will be an xarray Dataset containing the estimated capacity factor values
-for the specified wind turbine model over the given time period and region.
-
-
-How the Extrapolation Model Works
----------------------------------
-
-The model calculates hub height wind speed from MERRA2, extrapolating the variables in
-MERRA's tavg1_2d_slv_Nx data collection, which is a set of the time-averaged
-single-layer diagnostics.
-
-Specifically, the variables we use for extrapolation are: 2-m wind (U2M, V2M, in m/s),
-10-m wind (U10M, V10M), 50-m wind (U50M, V50M), and the zero-plane displacement
-height (DISPH, in meters). Additionally, we also use the wind speed at MERRA2's lowest
-model level (ULML, VLML, in m/s), the height of the lowest model level
-(HLML, in meters), may vary depending on the location. We can obtain the wind speed at
-any given location and height by computing the norm of the vector sum of the U and V
-components.
-
-
-The hub height wind speed can be calculated as
-
-.. math::
-   \nu  = \alpha \ln\left(\frac{H - d}{z}\right)
-
-.. math::
-   z = e^{-\beta/\alpha}
-
-where :math:`\nu` is the hub height wind speed, :math:`\alpha` is the best-fit slope
-from a linear regression of wind speeds on vertical heights, :math:`\ln` is the natural logarithm, :math:`H` is the hub height,
-:math:`d` is the zero-plane displacement height, and :math:`\beta` is the intercept
-from the linear regression fit.
-
-Here we estimate :math:`\alpha` and :math:`\beta` fitting a simple linear regression model to the heights and wind speeds in the data.
diff --git a/docs/source/modeling/wind/index.rst b/docs/source/modeling/wind/index.rst
index e6a93019..4ff14d70 100644
--- a/docs/source/modeling/wind/index.rst
+++ b/docs/source/modeling/wind/index.rst
@@ -1,9 +1,9 @@
 Wind Modeling
 =============
 
-Starting from geodata v0.2.0, geodata's capability to model and estimate wind speed have
-been from the cutout module to a separate wind module. This module has the capability to
-estimate wind speed with two modes: interpolation and extrapolation.
+Starting from geodata v0.2.0, geodata's wind modeling capability lives in a separate
+``geodata.model.wind`` module. The **supported ERA5 path** uses vertical spline
+**interpolation** on ``wind_3d_hourly`` data (see :doc:`interpolation`).
 
 How to use the models
 ---------------------
@@ -59,7 +59,33 @@ different dataset, you will need to create a new model.
     model = WindInterpolationModel(ds)
     model.prepare()  # Prepare the model
 
-    print(model.prepared) # Check if the model is prepared. Should return True.
+    print(model.prepared)  # Check if the model is prepared. Should return True.
+
+
+Preparing the model (``prepare``, ``prepared``, ``force``)
+----------------------------------------------------------
+
+Wind models must be **prepared** before ``estimate()``. Preparation reads the
+downloaded ERA5 files, computes month-by-month coefficients (B-spline parameters for
+interpolation), and writes cached results under ``GEODATA_ROOT/models/`` (see
+:doc:`/quick_start/packagesetup`).
+
+- ``model.prepared`` — ``True`` when every month in the model's year/month range has
+  cached outputs on disk.
+- ``model.prepare()`` — run once after ``ds.downloaded`` is ``True``. Safe to skip if
+  already prepared.
+- ``model.prepare(force=True)`` — delete and recompute cached months (use after changing
+  ``years`` / ``months`` / ``bounds`` on the source dataset, or when upgrading geodata).
+
+``estimate()`` raises if the model is not prepared. Pvlib does **not** use this
+prepare step; only wind models do.
+
+.. code:: Python
+
+    if not model.prepared:
+        model.prepare()
+    # After changing the source time range or domain:
+    # model.prepare(force=True)
 
 
 Once the model is prepared, we can use it to estimate wind speed at desired heights.
@@ -76,9 +102,58 @@ Once the model is prepared, we can use it to estimate wind speed at desired heig
 The above demonstrates the typical workflow. More model-specific details can be found
 in each model's respective tutorial as well as in the API reference.
 
+Estimate options
+----------------
+
+Wind models share the same ``estimate()`` subsetting interface (defined on
+``BaseModel`` in ``geodata.model``).
+
+Temporal subsetting
+~~~~~~~~~~~~~~~~~~~
+
+Use ``years`` and ``months`` slices to limit the estimation period. For example,
+``years=slice(2006, 2006), months=slice(1, 1)`` processes January 2006 only.
+
+Spatial subsetting (``xs``, ``ys``)
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+Pass ``xs`` and ``ys`` as ``slice(start, stop)`` to restrict longitude (``x``) and
+latitude (``y``). Omit either argument to use the full horizontal domain of the
+prepared source.
+
+Geodata **normalizes slice bounds** before ``xarray.Dataset.sel()``. You may pass
+``slice(high, low)`` or ``slice(low, high)``; the helper resolves the inclusive
+range and matches the coordinate's ascending or descending order (ERA5 latitude is
+typically descending). Without this, a slice like ``ys=slice(46, 48)`` on a
+descending ``y`` axis can incorrectly return an empty selection.
+
+.. code:: Python
+
+    # Subregion over central Europe — bounds order does not matter
+    wind_speed = model.estimate(
+        height=100.0,
+        years=slice(2006, 2006),
+        months=slice(1, 1),
+        xs=slice(8, 10),
+        ys=slice(48, 46),
+    )
+
+Wind-specific arguments
+~~~~~~~~~~~~~~~~~~~~~~~
+
+Pass **either**:
+
+- ``height=<meters>`` — hub-height or AGL wind speed (``WindInterpolationModel`` on
+  ``wind_3d_hourly``), or
+- ``turbine="<name>"`` — capacity factor (``cf``) from a turbine YAML under
+  ``geodata.resources.windturbine``. The name is the YAML stem (e.g.
+  ``Vestas_V112_3MW``). See :doc:`interpolation` Step 5 for usage and what
+  ``cf`` represents.
+
+List available turbines with ``geodata.resource.get_available_windturbines()``.
+
 .. toctree::
    :maxdepth: 1
    :caption: Tutorials on Specific Models
 
    interpolation
-   extrapolation
diff --git a/docs/source/modeling/wind/interpolation.rst b/docs/source/modeling/wind/interpolation.rst
index 788bba56..04751db6 100644
--- a/docs/source/modeling/wind/interpolation.rst
+++ b/docs/source/modeling/wind/interpolation.rst
@@ -127,16 +127,27 @@ by the original dataset, we can do this as follows:
     )
 
 
-This will return an xarray Dataset containing the estimated wind speed values. Note
-that you can also select a subset area by passing in :code:`xs=slice(start, end)`
-and/or :code:`ys=slice(start, end)` parameters to the `estimate` method.
+This will return an xarray Dataset containing the estimated wind speed values. You can
+also restrict the horizontal domain with ``xs`` and ``ys`` (see
+:doc:`/modeling/wind/index` — **Estimate options** for slice-order behavior on
+ERA5 grids).
+
+.. code:: Python
+
+    estimated_wind_speed = model.estimate(
+        height=60.0,
+        years=slice(2006, 2006),
+        months=slice(1, 1),
+        xs=slice(8, 10),
+        ys=slice(48, 46),
+    )
 
 Step 5: Estimate Wind Turbine Capacity Factor (CF) using the interpolation model
 --------------------------------------------------------------------------------
 
 Geodata also supports a limited set of wind turbine models to estimate the capacity
 factor (CF) of a wind turbine directly. To get a list of available wind turbine models,
-you can use the `get_available_windturbines` function:
+you can use the ``get_available_windturbines`` function:
 
 .. code:: Python
 
@@ -146,20 +157,39 @@ you can use the `get_available_windturbines` function:
     print(turbines)  # List of available wind turbine configurations
 
 
-To estimate the capacity factor of a wind turbine, you can use the `estimate` method
-and passign in the `turbine` parameter with the name of the wind turbine model.
+Pass the YAML **stem** (filename without ``.yaml``) as ``turbine`` — for example
+``Vestas_V112_3MW`` for ``src/geodata/resources/windturbine/Vestas_V112_3MW.yaml``.
 
 .. code:: Python
 
-    # Estimate the capacity factor for a specific wind turbine model
-    estimated_cf: xr.Dataset = model.estimate(
-        turbine="Vestas_V112_3MW",  # Example wind turbine model
+    estimated_cf = model.estimate(
+        turbine="Vestas_V112_3MW",
         years=slice(2006, 2006),
         months=slice(1, 1),
     )
 
-    print(estimated_cf)  # Display the estimated capacity factor
+    print(estimated_cf)
+
+Understanding the output
+~~~~~~~~~~~~~~~~~~~~~~
+
+``estimate(turbine=...)`` returns an ``xarray.DataArray`` named ``cf`` with dimensions
+``(time, x, y)`` when those coordinates are present.
+
+Geodata computes CF in three steps:
+
+1. **Hub-height wind speed** — interpolate to the turbine's ``HUB_HEIGHT`` from the
+   YAML (same vertical spline as Step 4, but at the turbine height rather than a
+   height you pass manually).
+2. **Power from the power curve** — map wind speed to power (MW) by interpolating the
+   tabulated ``V`` / ``POW`` pairs in the turbine YAML.
+3. **Normalize** — ``cf = power / P``, where ``P`` is the rated power (maximum value
+   in ``POW``).
 
+So ``cf`` is a **dimensionless capacity factor** in ``[0, 1]`` (values can exceed 1
+briefly if the curve extrapolates above rated power). Values outside the tabulated
+wind-speed range use SciPy's ``interp1d`` extrapolation — treat edge cases with care
+in sensitivity analysis.
 
-The output will be an xarray Dataset containing the estimated capacity factor values
-for the specified wind turbine model over the given time period and region.
+For implementation details, see ``WindBaseModel._estimate_power`` in the
+:ref:`API reference <modindex>`.
diff --git a/docs/source/quick_start/input_output.md b/docs/source/quick_start/input_output.md
index e08ddea0..03be7cb7 100644
--- a/docs/source/quick_start/input_output.md
+++ b/docs/source/quick_start/input_output.md
@@ -10,6 +10,9 @@
 
 ### MERRA2
 
+MERRA-2 is supported through the **legacy** ``Dataset`` / ``Cutout`` API only. See
+[Legacy workflow → MERRA2](../legacy/merra2/index.md) for download and cutout tutorials.
+
 * [MERRA2 hourly, single-level surface flux diagnostics](https://disc.gsfc.nasa.gov/datasets/M2T1NXFLX_5.12.4/summary)
 * [MERRA2 monthly mean, single-level surface flux diagnostics](https://disc.gsfc.nasa.gov/datasets/M2TMNXFLX_5.12.4/summary)
 * [MERRA2 daily mean, single-level diagnostics](https://disc.gsfc.nasa.gov/datasets/M2SDNXSLV_5.12.4/summary)
@@ -31,15 +34,15 @@ The following outputs are currently supported for climate data:
 
 **Wind**
 
-* Wind generation time-series ([MERRA2](../datasets/merra2/merra2_outputs.md#wind-generation-time-series), [ERA5](../datasets/era5/era5_outputs.md#wind-generation-time-series))
-* Wind speed time-series ([MERRA2](../datasets/merra2/merra2_outputs.md#wind-speed-time-series), [ERA5](../datasets/era5/era5_outputs.md#wind-speed-time-series))
-* Wind power density time-series ([MERRA2 only](../datasets/merra2/merra2_outputs.md#wind-power-density-time-series))
+* Wind generation time-series ([ERA5 model outputs](../modeling/era5_outputs.md#wind-generation-time-series), [wind modeling](../modeling/wind/index.rst), [ERA5 setup](../datasets/era5.rst))
+* Wind speed time-series ([ERA5 model outputs](../modeling/era5_outputs.md#wind-speed-time-series), [wind modeling](../modeling/wind/index.rst), [ERA5 setup](../datasets/era5.rst))
+* Wind power density time-series ([MERRA2 only (legacy)](../legacy/merra2/merra2_outputs.md#wind-power-density-time-series))
 
 
 **Solar**
 
-* Solar photovoltaic generation time-series ([ERA5 only](../datasets/era5/era5_outputs.md#solar-photovoltaic-generation-time-series))
-* PV generation time-series ([MERRA2 only](../datasets/merra2/merra2_outputs.md#pv-generation-time-series))
+* Solar photovoltaic generation time-series ([ERA5 model outputs](../modeling/era5_outputs.md#solar-photovoltaic-generation-time-series), [PVLib modeling](../modeling/pvlib/index.rst))
+* PV generation time-series ([MERRA2 only (legacy)](../legacy/merra2/merra2_outputs.md#pv-generation-time-series))
 
 
 **Temperature**
@@ -49,7 +52,7 @@ The following outputs are currently supported for climate data:
 
 **Aerosols**
 
-* PM2.5 time series ([MERRA2 only](../datasets/merra2/merra2_outputs.md#pm25-time-series))
+* PM2.5 time series ([MERRA2 only (legacy)](../legacy/merra2/merra2_outputs.md#pm25-time-series))
 
 ### Mask Specific