From 307fc9b9e25ad5503704123876c6122ed17410f9 Mon Sep 17 00:00:00 2001 From: KULcoder Date: Tue, 9 Jun 2026 19:04:45 -0700 Subject: [PATCH] docs: remove obsolete Jupyter notebooks from documentation - Deleted several outdated Jupyter notebooks, including `mask_on_cutout.ipynb`, `merra2.ipynb`, `mask_creation_workflow.ipynb`, `xarray_mask_tutorial.ipynb`, and `visualization.ipynb`, to streamline the documentation and remove redundant content. - This cleanup enhances the clarity and focus of the documentation, ensuring users have access to relevant and up-to-date resources. --- .../legacy/mask_on_cutout.ipynb | 444 ------- .../legacy/merra2/merra2.ipynb | 615 --------- .../mask/mask_creation_workflow.ipynb | 1151 ----------------- .../mask/xarray_mask_tutorial.ipynb | 331 ----- .../visualization/visualization.ipynb | 451 ------- 5 files changed, 2992 deletions(-) delete mode 100644 docs/jupyter_execute/legacy/mask_on_cutout.ipynb delete mode 100644 docs/jupyter_execute/legacy/merra2/merra2.ipynb delete mode 100644 docs/jupyter_execute/mask/mask_creation_workflow.ipynb delete mode 100644 docs/jupyter_execute/mask/xarray_mask_tutorial.ipynb delete mode 100644 docs/jupyter_execute/visualization/visualization.ipynb diff --git a/docs/jupyter_execute/legacy/mask_on_cutout.ipynb b/docs/jupyter_execute/legacy/mask_on_cutout.ipynb deleted file mode 100644 index b3d138a..0000000 --- a/docs/jupyter_execute/legacy/mask_on_cutout.ipynb +++ /dev/null @@ -1,444 +0,0 @@ -{ - "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 deleted file mode 100644 index 92b1cf1..0000000 --- a/docs/jupyter_execute/legacy/merra2/merra2.ipynb +++ /dev/null @@ -1,615 +0,0 @@ -{ - "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 deleted file mode 100644 index 9607f7a..0000000 --- a/docs/jupyter_execute/mask/mask_creation_workflow.ipynb +++ /dev/null @@ -1,1151 +0,0 @@ -{ - "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 deleted file mode 100644 index 1722242..0000000 --- a/docs/jupyter_execute/mask/xarray_mask_tutorial.ipynb +++ /dev/null @@ -1,331 +0,0 @@ -{ - "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 deleted file mode 100644 index 2542376..0000000 --- a/docs/jupyter_execute/visualization/visualization.ipynb +++ /dev/null @@ -1,451 +0,0 @@ -{ - "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