diff --git a/Changelog.rst b/Changelog.rst index a5eefa214e..efef383fb9 100644 --- a/Changelog.rst +++ b/Changelog.rst @@ -1,7 +1,7 @@ -Version NEXTVERSION +Version 3.19.0 -------------- -**2026-01-??** +**2026-01-19** * Write Zarr v3 datasets with `cf.write` (https://github.com/NCAS-CMS/cf-python/issues/895) @@ -9,7 +9,6 @@ Version NEXTVERSION `cf.read` (https://github.com/NCAS-CMS/cf-python/issues/894) * Reduce the time taken to import `cf` (https://github.com/NCAS-CMS/cf-python/issues/902) -* New optional dependency: ``zarr>=3.1.3`` * New function to control the creation of cached elements during data display: `cf.display_data` (https://github.com/NCAS-CMS/cf-python/issues/913) @@ -23,7 +22,8 @@ Version NEXTVERSION (https://github.com/NCAS-CMS/cf-python/issues/902) * Reduce the time taken to import `cf` (https://github.com/NCAS-CMS/cf-python/issues/902) -* Changed dependency: ``cfdm>=1.12.4.0, <1.12.5.0`` +* New optional dependency: ``zarr>=3.1.3`` +* Changed dependency: ``cfdm>=1.13.0.0, <1.13.1.0`` ---- @@ -1214,7 +1214,6 @@ version 2.2.4 version 2.2.3 -------------- ----- **2018-06-21** @@ -1253,7 +1252,6 @@ version 2.2.1 version 2.2.0 ------------- ----- **2018-06-04** @@ -1320,7 +1318,6 @@ version 2.1.3 version 2.1.2 ------------- ----- **2017-11-28** @@ -1353,7 +1350,6 @@ version 2.1 version 2.0.6 ------------- ----- **2017-09-28** diff --git a/README.md b/README.md index 140e6faeb9..44f2c3b51c 100644 --- a/README.md +++ b/README.md @@ -96,7 +96,6 @@ of its array manipulation and can: * test whether two field constructs are the same, * modify field construct metadata and data, * create subspaces of field constructs, -* write field constructs to netCDF datasets on disk, * incorporate, and create, metadata stored in external files, * read, write, and create data that have been compressed by convention (i.e. ragged or gathered arrays, or coordinate arrays compressed by diff --git a/RELEASE.md b/RELEASE.md index 0e6713120e..f76e96fcc4 100644 --- a/RELEASE.md +++ b/RELEASE.md @@ -26,8 +26,7 @@ - [ ] Ensure that the requirements on dependencies & their versions are up-to-date and consistent in both the `requirements.txt` and in `docs/source/installation.rst` (paying particular attention to - `cfdm`); and in the `_requires` list and `Version` checks in - `cf/__init__.py`. + `cfdm`); and in the cfdm-only `Version` checks in `cf/__init__.py`. - [ ] Make sure that `README.md` is up to date. @@ -40,7 +39,7 @@ deprecated methods and keyword arguments that can be completely removed, i.e. those with a ``removed_at`` version that is at or before the version being released. Remove any reference to them in - the method, class, or fucntion (including, if appropriate, the + the method, class, or function (including, if appropriate, the ``@_deprecated_kwarg_check`` decorator), and remove them from the relevant documentation ``.rst`` files. @@ -67,10 +66,6 @@ ./test_tutorial_code ``` -- [ ] **Follow all of the steps outlined externally in [`RECIPES.md`](./RECIPES.md)**, - notably so that the correct Sphinx-related environment is prepared for - documentation building. - - [ ] Ensure that the [PDF for Cheat Sheet](docs/_downloads/cheatsheet.pdf) is updated to include any API changes. The PDF is created using Canva keeping in mind the colours and fonts of the website. The same could diff --git a/cf/__init__.py b/cf/__init__.py index deaac4a813..125aa6286f 100644 --- a/cf/__init__.py +++ b/cf/__init__.py @@ -85,16 +85,16 @@ from packaging.version import Version -__date__ = "2025-10-16" -__version__ = "3.18.2" +__date__ = "2026-01-19" +__version__ = "3.19.0" __cf_version__ = cfdm.__cf_version__ __Conventions__ = f"CF-{__cf_version__}" # Check the version of cfdm (this is worth doing because of the very # tight coupling between cf and cfdm, and the risk of bad things # happening at run time if the versions are mismatched). -_minimum_vn = "1.12.3.1" -_maximum_vn = "1.12.4.0" +_minimum_vn = "1.13.0.0" +_maximum_vn = "1.13.1.0" _cfdm_vn = Version(cfdm.__version__) if _cfdm_vn < Version(_minimum_vn) or _cfdm_vn >= Version(_maximum_vn): raise RuntimeError( diff --git a/cf/functions.py b/cf/functions.py index 6b1850a388..c4f7abe5fb 100644 --- a/cf/functions.py +++ b/cf/functions.py @@ -251,7 +251,7 @@ def configuration( The new display data option. The default is to not change the current behaviour. - .. versionadded:: NEXTVERSION + .. versionadded:: 3.19.0 regrid_logging: `bool` or `Constant`, optional The new value (either True to enable logging or False to @@ -1504,7 +1504,7 @@ def total_memory(): def is_log_level_info(logger): """Return True if and only if log level is at least as verbose as INFO. - Deprecated at version NEXTVERSION and is no longer available. Use + Deprecated at version 3.19.0 and is no longer available. Use `cfdm.is_log_level_info` instead. .. versionadded:: 3.16.3 @@ -1525,7 +1525,7 @@ def is_log_level_info(logger): _DEPRECATION_ERROR_FUNCTION( "is_log_level_info", message="Use cfdm.is_log_level_info instead", - version="NEXTVERSION", + version="3.19.0", removed_at="5.0.0", ) # pragma: no cover diff --git a/cf/read_write/read.py b/cf/read_write/read.py index 52cc132e7e..c920b72b97 100644 --- a/cf/read_write/read.py +++ b/cf/read_write/read.py @@ -318,7 +318,7 @@ class read(cfdm.read): {{read store_dataset_shards: `bool`, optional}} - .. versionadded:: NEXTVERSION + .. versionadded:: 3.19.0 {{read cfa: `dict`, optional}} @@ -334,7 +334,7 @@ class read(cfdm.read): {{read group_dimension_search: `str`, optional}} - .. versionadded:: NEXTVERSION + .. versionadded:: 3.19.0 umversion: deprecated at version 3.0.0 Use the *um* parameter instead. diff --git a/docs/_templates/sponsors.html b/docs/_templates/sponsors.html index e5e4d4f6be..0ca7485220 100644 --- a/docs/_templates/sponsors.html +++ b/docs/_templates/sponsors.html @@ -16,8 +16,8 @@ and by NCAS.

- - - + + +
- + diff --git a/docs/_templates/logo_EC.png b/docs/source/_static/logo_EC.png similarity index 100% rename from docs/_templates/logo_EC.png rename to docs/source/_static/logo_EC.png diff --git a/docs/_templates/logo_ERC.png b/docs/source/_static/logo_ERC.png similarity index 100% rename from docs/_templates/logo_ERC.png rename to docs/source/_static/logo_ERC.png diff --git a/docs/_templates/logo_NCAS.png b/docs/source/_static/logo_NCAS.png similarity index 100% rename from docs/_templates/logo_NCAS.png rename to docs/source/_static/logo_NCAS.png diff --git a/docs/_templates/logo_NERC.png b/docs/source/_static/logo_NERC.png similarity index 100% rename from docs/_templates/logo_NERC.png rename to docs/source/_static/logo_NERC.png diff --git a/docs/source/aggregation_rules.rst b/docs/source/aggregation_rules.rst index c9732ee35e..6015cc4952 100644 --- a/docs/source/aggregation_rules.rst +++ b/docs/source/aggregation_rules.rst @@ -6,8 +6,6 @@ **Aggregation rules** ===================== ----- - Version |release| for version |version| of the CF conventions. *David Hassell and Jonathan Gregory (2012)* diff --git a/docs/source/api_reference.rst b/docs/source/api_reference.rst index 24a60e03de..af73b1a377 100644 --- a/docs/source/api_reference.rst +++ b/docs/source/api_reference.rst @@ -4,8 +4,6 @@ **API reference** ================= ----- - Version |release| for version |version| of the CF conventions. * **Construct classes** diff --git a/docs/source/cheat_sheet.rst b/docs/source/cheat_sheet.rst index 8a9b8dead7..dce63c93e0 100644 --- a/docs/source/cheat_sheet.rst +++ b/docs/source/cheat_sheet.rst @@ -7,8 +7,6 @@ **Cheat Sheet** =============== ----- - Version |release| for version |version| of the CF conventions. This cheat sheet provides a summary of some key functions and methods in diff --git a/docs/source/class/cf.Data.rst b/docs/source/class/cf.Data.rst index b89f262de9..41b4c787b7 100644 --- a/docs/source/class/cf.Data.rst +++ b/docs/source/class/cf.Data.rst @@ -659,6 +659,11 @@ Performance ~cf.Data.fits_in_memory ~cf.Data.section ~cf.Data.persist + ~cf.Data.cache_elements + ~cf.Data.get_cached_elements + ~cf.Data.nc_clear_dataset_shards + ~cf.Data.nc_dataset_shards + ~cf.Data.nc_set_dataset_shards .. rubric:: Attributes diff --git a/docs/source/function.rst b/docs/source/function.rst index 8892f325c1..31de3bce8e 100644 --- a/docs/source/function.rst +++ b/docs/source/function.rst @@ -20,7 +20,6 @@ Reading and writing cf.read cf.write - cf.netcdf_lock Aggregation ----------- @@ -157,7 +156,6 @@ Active storage reductions cf.active_storage cf.active_storage_url cf.active_storage_max_requests - cf.netcdf_lock Miscellaneous ------------- diff --git a/docs/source/installation.rst b/docs/source/installation.rst index 4ec36e0d67..13c9f0ef88 100644 --- a/docs/source/installation.rst +++ b/docs/source/installation.rst @@ -225,8 +225,8 @@ Required * `scipy `_, version 1.10.0 or newer. -* `cfdm `_, version 1.12.3.1 or up to, - but not including, 1.12.4.0. +* `cfdm `_, version 1.13.0.0 or up to, + but not including, 1.13.1.0. * `cfunits `_, version 3.3.7 or newer. diff --git a/docs/source/introduction.rst b/docs/source/introduction.rst index 9e9b3208a9..537ff8a088 100644 --- a/docs/source/introduction.rst +++ b/docs/source/introduction.rst @@ -93,8 +93,6 @@ The `cf` package can: * create subspaces of field constructs, -* write field constructs to netCDF datasets on disk, - * incorporate, and create, metadata stored in external files, * read, write, and create data that have been compressed by convention diff --git a/docs/source/releases.rst b/docs/source/releases.rst index ff5fe0923b..74986da574 100644 --- a/docs/source/releases.rst +++ b/docs/source/releases.rst @@ -13,10 +13,15 @@ Documentation for all versions of cf. :local: :backlinks: entry +**CF-1.13** +----------- + +* `Version 3.19.0 `_ (2026-01-19) + **CF-1.12** ----------- -* `Version 3.18.2 `_ (2025-10-16) +* `Version 3.18.2 `_ (2025-10-16) * `Version 3.18.1 `_ (2025-08-20) * `Version 3.18.0 `_ (2025-06-05) * `Version 3.17.0 `_ (2025-04-02) diff --git a/docs/source/spelling_false_positives.txt b/docs/source/spelling_false_positives.txt index 656661c165..36d5ada565 100644 --- a/docs/source/spelling_false_positives.txt +++ b/docs/source/spelling_false_positives.txt @@ -449,8 +449,9 @@ seterr setitem setmask setprop -Signell shapefiles +sharding +Signell sinh spacings spatiotemporal diff --git a/recipes-docs/source/conf.py b/recipes-docs/source/conf.py index e5bfa04418..cc6782ef0b 100644 --- a/recipes-docs/source/conf.py +++ b/recipes-docs/source/conf.py @@ -1,8 +1,11 @@ +import faulthandler import os import sys from sphinx_gallery.sorting import FileNameSortKey +faulthandler.enable() # seg fault detection - as recipes prone to seg faulting + # Make main 'docs' conf.py importable sys.path.insert(0, os.path.abspath("../../docs/source")) @@ -19,8 +22,9 @@ # sphinx-gallery configuration sphinx_gallery_conf = { - "examples_dirs": "recipes", # path to recipe files + "examples_dirs": "recipes-source", # path to recipe files "gallery_dirs": "recipes", # path to save gallery generated output + "ignore_pattern": "/exclusions/", "run_stale_examples": False, # Below setting can be buggy: see: # https://github.com/sphinx-gallery/sphinx-gallery/issues/967 @@ -29,15 +33,22 @@ "doc_module": ("cf",), "inspect_global_variables": True, "within_subsection_order": FileNameSortKey, - "default_thumb_file": "_static/cf-recipe-placeholder-squarecrop.png", + "default_thumb_file": "../../docs/source/_static/cf-recipe-placeholder-squarecrop.png", "image_scrapers": ( "matplotlib", ), # Ensures Matplotlib images are captured "plot_gallery": True, # Enables plot rendering "reset_modules": ("matplotlib",), # Helps with memory management "capture_repr": (), + # "filename_pattern": r"plot", } +exclude_patterns = [ + "exclusions/**", +] + html_static_path = ["../../docs/source/_static"] html_logo = "../../docs/source/images/logo.svg" html_favicon = "../../docs/source/_static/favicon.ico" + +templates_path = ["../../docs/_templates"] diff --git a/recipes-docs/source/index.rst b/recipes-docs/source/index.rst index ee2bddf946..3bbd15b42f 100644 --- a/recipes-docs/source/index.rst +++ b/recipes-docs/source/index.rst @@ -1,95 +1,9 @@ -:orphan: - **Recipes using cf** ==================== ---- -Version |release| for version |version| of the CF conventions. - -Click on the keywords below to filter the recipes according to their function: - -.. raw:: html - - - -.. raw:: html - - - -
- - - - - - - - - - - - - -
- -.. raw:: html - - -.. raw:: html - -
- - -.. raw:: html - -
- - -.. only:: html - - .. container:: sphx-glr-footer sphx-glr-footer-gallery - - .. container:: sphx-glr-download sphx-glr-download-python - - :download:`Download all examples in Python source code: recipes_python.zip ` - - .. container:: sphx-glr-download sphx-glr-download-jupyter - - :download:`Download all examples in Jupyter notebooks: recipes_jupyter.zip ` - - -.. only:: html - - .. rst-class:: sphx-glr-signature +.. toctree:: + :maxdepth: 2 - `Gallery generated by Sphinx-Gallery `_ + recipes-source/gallery diff --git a/recipes-docs/source/recipes-source/README.rst b/recipes-docs/source/recipes-source/README.rst new file mode 100644 index 0000000000..69ef485e26 --- /dev/null +++ b/recipes-docs/source/recipes-source/README.rst @@ -0,0 +1,64 @@ +**Recipes using cf** +==================== + +---- + +Version |release| for version |version| of the CF conventions. + +Click on the keywords below to filter the recipes according to their function: + +.. raw:: html + + + +.. raw:: html + + + +
+ + + + + + + + + + + + + +
+ +.. raw:: html + diff --git a/recipes-docs/source/recipes-source/cfplot.png b/recipes-docs/source/recipes-source/cfplot.png new file mode 100644 index 0000000000..9e4718c22d Binary files /dev/null and b/recipes-docs/source/recipes-source/cfplot.png differ diff --git a/recipes-docs/source/recipes-source/gallery.rst b/recipes-docs/source/recipes-source/gallery.rst new file mode 100644 index 0000000000..511e268f52 --- /dev/null +++ b/recipes-docs/source/recipes-source/gallery.rst @@ -0,0 +1,121 @@ +**Recipes using cf** +==================== + +---- + +Version |release| for version |version| of the CF conventions. + +Click on the keywords below to filter the recipes according to their function: + +.. raw:: html + + + +.. raw:: html + + + +
+ + + + + + + + + + + + + +
+ +.. raw:: html + + + + +.. raw:: html + +
+ +.. thumbnail-parent-div-open + +.. raw:: html + +
+ +.. only:: html + + .. image:: /recipes/images/thumb/sphx_glr_plot_22_recipe_thumb.png + :alt: + + :doc:`/recipes/plot_22_recipe` + +.. raw:: html + +
Plotting a wind rose as a scatter plot
+
+ + +.. thumbnail-parent-div-close + +.. raw:: html + +
+ + +.. toctree:: + :hidden: + + /recipes/plot_22_recipe + + +.. only:: html + + .. container:: sphx-glr-footer sphx-glr-footer-gallery + + .. container:: sphx-glr-download sphx-glr-download-python + + :download:`Download all examples in Python source code: recipes_python.zip ` + + .. container:: sphx-glr-download sphx-glr-download-jupyter + + :download:`Download all examples in Jupyter notebooks: recipes_jupyter.zip ` + + +.. only:: html + + .. rst-class:: sphx-glr-signature + + `Gallery generated by Sphinx-Gallery `_ diff --git a/recipes-docs/source/recipes-source/index.rst b/recipes-docs/source/recipes-source/index.rst new file mode 100644 index 0000000000..646c716435 --- /dev/null +++ b/recipes-docs/source/recipes-source/index.rst @@ -0,0 +1,7 @@ +cf-python recipes +================= + +.. toctree:: + :maxdepth: 2 + + recipes/gallery diff --git a/recipes-docs/source/recipes-source/plot_01_recipe.py b/recipes-docs/source/recipes-source/plot_01_recipe.py new file mode 100644 index 0000000000..93b9dde417 --- /dev/null +++ b/recipes-docs/source/recipes-source/plot_01_recipe.py @@ -0,0 +1,77 @@ +""" +Calculating global mean temperature timeseries +============================================== + +In this recipe we will calculate and plot monthly and annual global mean temperature timeseries. +""" + +# %% +# 1. Import cf-python and cf-plot: + +import cfplot as cfp + +import cf + +# %% +# 2. Read the field constructs: + +f = cf.read("~/recipes/cru_ts4.06.1901.2021.tmp.dat.nc") +print(f) + +# %% +# 3. Select near surface temperature by index and look at its contents: + +temp = f[1] +print(temp) + +# %% +# 4. Select latitude and longitude dimensions by identities, with two different techniques: + +lon = temp.coordinate("long_name=longitude") +lat = temp.coordinate("Y") + +# %% +# 5. Print the description of near surface temperature using the dump method to show properties of all constructs: + +temp.dump() + +# %% +# 6. Latitude and longitude dimension coordinate cell bounds are absent, which are created and set: + +a = lat.create_bounds() +lat.set_bounds(a) +lat.dump() + +# %% + +b = lon.create_bounds() +lon.set_bounds(b) +lon.dump() + +# %% + +print(b.array) + +# %% +# 7. Time dimension coordinate cell bounds are similarly created and set for cell sizes of one calendar month: + +time = temp.coordinate("long_name=time") +c = time.create_bounds(cellsize=cf.M()) +time.set_bounds(c) +time.dump() + +# %% +# 8. Calculate and plot the area weighted mean surface temperature for each time: + +global_avg = temp.collapse("area: mean", weights=True) +cfp.lineplot(global_avg, color="red", title="Global mean surface temperature") + +# %% +# 9. Calculate and plot the annual global mean surface temperature: + +annual_global_avg = global_avg.collapse("T: mean", group=cf.Y()) +cfp.lineplot( + annual_global_avg, + color="red", + title="Annual global mean surface temperature", +) diff --git a/recipes-docs/source/recipes-source/plot_02_recipe.py b/recipes-docs/source/recipes-source/plot_02_recipe.py new file mode 100644 index 0000000000..902c796788 --- /dev/null +++ b/recipes-docs/source/recipes-source/plot_02_recipe.py @@ -0,0 +1,96 @@ +""" +Calculating and plotting the global average temperature anomalies +================================================================= + +In this recipe we will calculate and plot the global average temperature anomalies. +""" + +# %% +# 1. Import cf-python and cf-plot: + +import cfplot as cfp + +import cf + +# %% +# 2. Read the field constructs: + +f = cf.read("~/recipes/cru_ts4.06.1901.2021.tmp.dat.nc") +print(f) + +# %% +# 3. Select near surface temperature by index and look at its contents: + +temp = f[1] +print(temp) + +# %% +# 4. Select latitude and longitude dimensions by identities, with two different techniques: + +lon = temp.coordinate("long_name=longitude") +lat = temp.coordinate("Y") + +# %% +# 5. Print the description of near surface temperature to show properties of all constructs: + +temp.dump() + +# %% +# 6. Latitude and longitude dimension coordinate cell bounds are absent, which are created and set: + +a = lat.create_bounds() +lat.set_bounds(a) +lat.dump() + +# %% + +b = lon.create_bounds() +lon.set_bounds(b) +lon.dump() + +# %% + +print(b.array) + +# %% +# 7. Time dimension coordinate cell bounds are similarly created and set for cell sizes of one calendar month: + +time = temp.coordinate("long_name=time") +c = time.create_bounds(cellsize=cf.M()) +time.set_bounds(c) +time.dump() + +# %% +# 8. Calculate the area weighted mean surface temperature for each time using the collapse method: + +global_avg = temp.collapse("area: mean", weights=True) + +# %% +# 9. Calculate the annual global mean surface temperature: + +annual_global_avg = global_avg.collapse("T: mean", group=cf.Y()) + +# %% +# 10. The temperature values are averaged for the climatological period of 1961-1990 by defining a subspace within these years using `cf.wi` query instance over subspace and doing a statistical collapse with the collapse method: + +annual_global_avg_61_90 = annual_global_avg.subspace( + T=cf.year(cf.wi(1961, 1990)) +) +print(annual_global_avg_61_90) + +# %% + +temp_clim = annual_global_avg_61_90.collapse("T: mean") +print(temp_clim) + +# %% +# 11. The temperature anomaly is then calculated by subtracting these climatological temperature values from the annual global average temperatures and plotted: + +temp_anomaly = annual_global_avg - temp_clim +cfp.lineplot( + temp_anomaly, + color="red", + title="Global Average Temperature Anomaly (1901-2021)", + ylabel="1961-1990 climatology difference ", + yunits="degree Celcius", +) diff --git a/recipes-docs/source/recipes-source/plot_03_recipe.py b/recipes-docs/source/recipes-source/plot_03_recipe.py new file mode 100644 index 0000000000..25633b1c54 --- /dev/null +++ b/recipes-docs/source/recipes-source/plot_03_recipe.py @@ -0,0 +1,35 @@ +""" +Plotting global mean temperatures spatially +=========================================== + +In this recipe, we will plot the global mean temperature spatially. +""" + +# %% +# 1. Import cf-python and cf-plot: + +import cfplot as cfp + +import cf + +# %% +# 2. Read the field constructs: + +f = cf.read("~/recipes/cru_ts4.06.1901.2021.tmp.dat.nc") +print(f) + +# %% +# 3. Select near surface temperature by index and look at its contents: + +temp = f[1] +print(temp) + +# %% +# 4. Average the monthly mean surface temperature values by the time axis using the collapse method: + +global_avg = temp.collapse("mean", axes="long_name=time") + +# %% +# 5. Plot the global mean surface temperatures: + +cfp.con(global_avg, lines=False, title="Global mean surface temperature") diff --git a/recipes-docs/source/recipes-source/plot_04_recipe.py b/recipes-docs/source/recipes-source/plot_04_recipe.py new file mode 100644 index 0000000000..a50296bb9e --- /dev/null +++ b/recipes-docs/source/recipes-source/plot_04_recipe.py @@ -0,0 +1,40 @@ +""" +Comparing two datasets with different resolutions using regridding +================================================================== + +In this recipe, we will regrid two different datasets with different resolutions. An example use case could be one where the observational dataset with a higher resolution needs to be regridded to that of the model dataset so that they can be compared with each other. +""" + +# %% +# 1. Import cf-python: + +import cf + +# %% +# 2. Read the field constructs: + +obs = cf.read("~/recipes/cru_ts4.06.1901.2021.tmp.dat.nc", dask_chunks=None) +print(obs) + +# %% + +model = cf.read( + "~/recipes/tas_Amon_HadGEM3-GC3-1_hist-1p0_r3i1p1f2_gn_185001-201412.nc" +) +print(model) + +# %% +# 3. Select observation and model temperature fields by identity and index respectively, and look at their contents: + +obs_temp = obs.select_field("long_name=near-surface temperature") +print(obs_temp) + +# %% + +model_temp = model[0] +print(model_temp) + +# %% +# 4. Regrid observational data to that of the model data and create a new low resolution observational data using bilinear interpolation: +obs_temp_regrid = obs_temp.regrids(model_temp, method="linear") +print(obs_temp_regrid) diff --git a/recipes-docs/source/recipes-source/plot_05_recipe.py b/recipes-docs/source/recipes-source/plot_05_recipe.py new file mode 100644 index 0000000000..e40f1c23ad --- /dev/null +++ b/recipes-docs/source/recipes-source/plot_05_recipe.py @@ -0,0 +1,64 @@ +""" +Plotting wind vectors overlaid on precipitation data +==================================================== + +In this recipe we will plot wind vectors, derived from northward and eastward wind components, over precipitation data. +""" + +# %% +# 1. Import cf-python and cf-plot: + +import cfplot as cfp + +import cf + +# %% +# 2. Read the field constructs: + +f1 = cf.read("~/recipes/northward.nc") +print(f1) + +# %% + +f2 = cf.read("~/recipes/eastward.nc") +print(f2) + +# %% + +f3 = cf.read("~/recipes/monthly_precipitation.nc") +print(f3) + +# %% +# 3. Select wind vectors and precipitation data by index and look at their contents: +v = f1[0] +print(v) + +# %% + +u = f2[0] +print(u) + +# %% + +pre = f3[0] +print(pre) + +# %% +# 4. Plot the wind vectors on top of precipitation data for June 1995 by creating a subspace with a date-time object and using `cfplot.con `_. Here `cfplot.gopen `_ is used to define the parts of the plot area, which is closed by `cfplot.gclose `_; `cfplot.cscale `_ is used to choose one of the colour maps amongst many available; `cfplot.levs `_ is used to set the contour levels for precipitation data; and `cfplot.vect `_ is used to plot the wind vectors for June 1995: +june_95 = cf.year(1995) & cf.month(6) +cfp.gopen() +cfp.cscale("precip4_11lev") +cfp.levs(step=100) +cfp.con( + pre.subspace(T=june_95), + lines=False, + title="June 1995 monthly global precipitation", +) +cfp.vect( + u=u.subspace(T=june_95), + v=v.subspace(T=june_95), + key_length=10, + scale=35, + stride=5, +) +cfp.gclose() diff --git a/recipes-docs/source/recipes-source/plot_06_recipe.py b/recipes-docs/source/recipes-source/plot_06_recipe.py new file mode 100644 index 0000000000..91ca4ab100 --- /dev/null +++ b/recipes-docs/source/recipes-source/plot_06_recipe.py @@ -0,0 +1,65 @@ +""" +Converting from rotated latitude-longitude to regular latitude-longitude +======================================================================== + +In this recipe, we will be regridding from a rotated latitude-longitude source domain to a regular latitude-longitude destination domain. +""" + +# %% +# 1. Import cf-python, cf-plot and numpy: + +import cfplot as cfp + +import cf + +# %% +# 2. Read the field constructs using read function: + +f = cf.read("~/recipes/au952a.pd20510414.pp") +print(f) + +# %% +# 3. Select the field by index and print its description to show properties of all constructs: + +gust = f[0] +gust.dump() + +# %% +# 4. Access the time coordinate of the gust field and retrieve the datetime values of the time coordinate: + +print(gust.coordinate("time").datetime_array) + +# %% +# 5. Create a new instance of the `cf.dt` class with a specified year, month, day, hour, minute, second and microsecond. Then store the result in the variable ``test``: +test = cf.dt(2051, 4, 14, 1, 30, 0, 0) +print(test) + +# %% +# 6. Plot the wind gust by creating a subspace for the specified variable ``test`` using `cfplot.con `_. Here `cfplot.mapset `_ is used to set the mapping parameters like setting the map resolution to 50m: +cfp.mapset(resolution="50m") +cfp.con(gust.subspace(T=test), lines=False) + +# %% +# 7. To see the rotated pole data on the native grid, the above steps are repeated and projection is set to rotated in `cfplot.mapset `_: +cfp.mapset(resolution="50m", proj="rotated") +cfp.con(gust.subspace(T=test), lines=False) + +# %% +# 8. Create dimension coordinates for the destination grid with the latitude and +# longitude values for Europe. `cf.Domain.create_regular +# `_ +# method is used to +# create a regular grid with longitudes and latitudes. Spherical regridding is +# then performed on the gust variable by passing the target domain as argument. +# The method also takes an argument ``'linear'`` which specifies the type of +# regridding method to use. The description of the ``regridded_data`` is finally +# printed to show properties of all its constructs: + +target_domain = cf.Domain.create_regular((-25, 45, 10), (32, 72, 10)) +regridded_data = gust.regrids(target_domain, "linear") +regridded_data.dump() + +# %% +# 9. Step 6 is similarly repeated for the ``regridded_data`` to plot the wind gust on a regular latitude-longitude domain: +cfp.mapset(resolution="50m") +cfp.con(regridded_data.subspace(T=test), lines=False) diff --git a/recipes-docs/source/recipes-source/plot_07_recipe.py b/recipes-docs/source/recipes-source/plot_07_recipe.py new file mode 100644 index 0000000000..02908e83ee --- /dev/null +++ b/recipes-docs/source/recipes-source/plot_07_recipe.py @@ -0,0 +1,65 @@ +""" +Plotting members of a model ensemble +==================================== + +In this recipe, we will plot the members of a model ensemble. +""" + +# %% +# 1. Import cf-python and cf-plot: + +import cfplot as cfp + +import cf + +# %% +# 2. Read the field constructs using read function and store it in the variable ``f``. The * in the filename is a wildcard character which means the function reads all files in the directory that match the specified pattern. [0:5] selects the first five elements of the resulting list: + +f = cf.read("~/recipes/realization/PRMSL.1941_mem*.nc")[0:5] +print(f) + +# %% +# 3. The description of one of the fields from the list shows ``'realization'`` as a property by which the members of the model ensemble are labelled: + +f[1].dump() + +# %% +# 4. An ensemble of the members is then created by aggregating the data in ``f`` along a new ``'realization'`` axis using the cf.aggregate function, and storing the result in the variable ``ensemble``. ``'relaxed_identities=True'`` allows for missing coordinate identities to be inferred. [0] selects the first element of the resulting list. ``id%realization`` now shows as an auxiliary coordinate for the ensemble: + +ensemble = cf.aggregate( + f, dimension=("realization",), relaxed_identities=True +)[0] +print(ensemble) + + +# %% +# 5. To see the constructs for the ensemble, print the *constructs* attribute: + +print(ensemble.constructs) + +# %% +# 6. Loop over the realizations in the ensemble using the *range* function and the *domain_axis* to determine the size of the realization dimension. For each realization, extract a subspace of the ensemble using the *subspace* method and the ``'id%realization'`` keyword argument along a specific latitude and longitude and plot the realizations from the 4D field using `cfplot.lineplot `_. +# A moving average of the ensemble along the time axis, with a window size of 90 (i.e. an approximately 3-month moving average) is calculated using the *moving_window* method. The ``mode='nearest'`` parameter is used to specify how to pad the data outside of the time range. The *squeeze* method removes any dimensions of size 1 from the field to produce a 2D field: + +cfp.gopen() + +for realization in range(1, ensemble.domain_axis("id%realization").size + 1): + cfp.lineplot( + ensemble.subspace( + **{"id%realization": realization}, latitude=[0], longitude=[0] + ).squeeze(), + label=f"Member {realization}", + linewidth=1.0, + ) + +cfp.lineplot( + ensemble.moving_window( + method="mean", window_size=90, axis="T", mode="nearest" + )[0, :, 0, 0].squeeze(), + label="Ensemble mean", + linewidth=2.0, + color="black", + title="Model Ensemble Pressure", +) + +cfp.gclose() diff --git a/recipes-docs/source/recipes-source/plot_08_recipe.py b/recipes-docs/source/recipes-source/plot_08_recipe.py new file mode 100644 index 0000000000..63427f62a7 --- /dev/null +++ b/recipes-docs/source/recipes-source/plot_08_recipe.py @@ -0,0 +1,143 @@ +""" +Plotting statistically significant temperature trends with stippling +==================================================================== + +In this recipe, we will analyse and plot temperature trends from the HadCRUT.5.0.1.0 dataset for two different time periods. The plotted maps also include stippling, which is used to highlight areas where the temperature trends are statistically significant. +""" + +# %% +# 1. Import cf-python, cf-plot, numpy and scipy.stats: + +import cfplot as cfp +import cf + +import numpy as np +import scipy.stats as stats + + +# %% +# 2. Three functions are defined: + +# %% +# * ``linear_trend(data, time_axis)``: This function calculates the linear regression slope and p-value for the input data along the time axis. It takes two arguments: ``'data'``, which represents the temperature anomalies or any other data you want to analyse, and ``'time_axis'``, which represents the corresponding time points for the data. The function uses the `stats.linregress `_ method from the `scipy.stats `_ library to calculate the slope and p-value of the linear regression. It returns these two values as a tuple: + + +def linear_trend(data, time_axis): + slope, _, _, p_value, _ = stats.linregress(time_axis, data) + return slope, p_value + + +# %% +# * ``create_trend_stipple_obj(temp_data, input_data)``: This function creates a new object with the input data provided and *collapses* the time dimension by taking the mean. It takes two arguments: ``'temp_data'``, which represents the temperature data object, and ``'input_data'``, which is the data to be set in the new object. The function creates a copy of the ``'temp_data'`` object by selecting the first element using index and *squeezes* it to remove the size 1 axis. It then sets the input data with the ``'Y'`` (latitude) and ``'X'`` (longitude) axes, and then *collapses* the time dimension using the ``"T: mean"`` operation: + + +def create_trend_stipple_obj(temp_data, input_data): + trend_stipple_obj = temp_data[0].squeeze() + trend_stipple_obj.set_data(input_data, axes=["Y", "X"]) + return trend_stipple_obj + + +# %% +# * ``process_subsets(subset_mask)``: This function processes the subsets of data by applying the ``linear_trend`` function along a specified axis. It takes one argument, ``'subset_mask'``, which is a boolean mask representing the time points to be considered in the analysis. The function first extracts the masked subset of data and then applies the ``linear_trend`` function along the time axis (axis 0) using the `numpy.ma.apply_along_axis `_ function. The result is an array containing the slope and p-value for each grid point in the dataset: + + +def process_subsets(subset_mask): + subset_data = masked_data[subset_mask, :, :] + return np.ma.apply_along_axis( + linear_trend, 0, subset_data, time_axis[subset_mask] + ) + + +# %% +# 3. Read the field constructs: + +temperature_data = cf.read( + "~/recipes/HadCRUT.5.0.1.0.analysis.anomalies.ensemble_mean.nc" +)[0] +print(temperature_data) + +# %% +# 4. Calculate the annual mean temperature anomalies. The ``'weights=True'`` argument is used take the varying lengths of months into account which ensures that the calculated mean is more accurate. A masked array is created for the annual mean temperature anomalies, masking any invalid values: + +annual_temperature = temperature_data.collapse( + "T: mean", weights=True, group=cf.Y() +) +time_axis = annual_temperature.coordinate("T").year.array +masked_data = np.ma.masked_invalid(annual_temperature.array) + +# %% +# 5. Define two time periods for analysis: 1850-2020 and 1980-2020, along with a significance level (alpha) of 0.05: + +time_periods = [(1850, 2020, "sub_1850_2020"), (1980, 2020, "sub_1980_2020")] +alpha = 0.05 +results = {} + +# %% +# 6. Loop through the time periods, processing the subsets, calculating trend p-values, and creating stipple objects. For each time period, the script calculates the trends and p-values using the ``process_subsets`` function. If the p-value is less than the significance level (alpha = 0.05), a stippling mask is created. The script then creates a new object for the trend and stippling mask using the ``create_trend_stipple_obj`` function: + +for start, end, prefix in time_periods: + subset_mask = (time_axis >= start) & (time_axis <= end) + subset_trend_pvalue = process_subsets(subset_mask) + results[prefix + "_trend_pvalue"] = subset_trend_pvalue + results[prefix + "_stipple"] = subset_trend_pvalue[1] < alpha + results[prefix + "_trend"] = create_trend_stipple_obj( + temperature_data, subset_trend_pvalue[0] + ) + results[prefix + "_stipple_obj"] = create_trend_stipple_obj( + temperature_data, results[prefix + "_stipple"] + ) + +# %% +# 7. Create two plots - one for the 1850-2020 time period and another for the 1980-2020 time period using `cfplot.con `_. +# The results are multiplied by 10 so that each plot displays the temperature trend in K/decade with stippling to indicate areas where the trend is statistically significant (p-value < 0.05). +# Here `cfplot.gopen `_ is used to define the parts of the plot area with two rows and one column, and setting the bottom margin to 0.2. +# It is closed by `cfplot.gclose `_; +# `cfplot.gpos `_ is used to set the plotting position of both the plots; +# `cfplot.mapset `_ is used to set the map projection to Robinson; +# `cfplot.cscale `_ is used to choose one of the colour maps amongst many available; +# `cfplot.levs `_ is used to set the contour levels; +# and `cfplot.stipple `_ is used to add stippling to show statistically significant areas: + +cfp.gopen(rows=2, columns=1, bottom=0.2) + +cfp.gpos(1) +cfp.mapset(proj="robin") +cfp.cscale("temp_19lev") +cfp.levs(min=-1, max=1, step=0.1) +cfp.con( + results["sub_1850_2020_trend"] * 10, + lines=False, + colorbar=None, + title="Temperature Trend 1850-2020", +) +cfp.stipple( + results["sub_1850_2020_stipple_obj"], + min=1, + max=1, + size=5, + color="k", + marker=".", +) + +cfp.gpos(2) +cfp.mapset(proj="robin") +cfp.cscale("temp_19lev") +cfp.levs(min=-1, max=1, step=0.1) +cfp.con( + results["sub_1980_2020_trend"] * 10, + lines=False, + title="Temperature Trend 1980-2020", + colorbar_position=[0.1, 0.1, 0.8, 0.02], + colorbar_orientation="horizontal", + colorbar_title="K/decade", +) +cfp.stipple( + results["sub_1980_2020_stipple_obj"], + min=1, + max=1, + size=5, + color="k", + marker=".", +) + +cfp.gclose() diff --git a/recipes-docs/source/recipes-source/plot_09_recipe.py b/recipes-docs/source/recipes-source/plot_09_recipe.py new file mode 100644 index 0000000000..42fe49cbd9 --- /dev/null +++ b/recipes-docs/source/recipes-source/plot_09_recipe.py @@ -0,0 +1,82 @@ +""" +Plotting a joint histogram +========================== + +In this recipe, we will be creating a joint histogram of PM2.5 mass +concentration and 2-metre temperature. +""" + +# %% +# 1. Import cf-python, numpy and matplotlib.pyplot: + +import matplotlib.pyplot as plt +import numpy as np + +import cf + +# %% +# 2. Read the field constructs using read function: +f = cf.read("~/recipes/levtype_sfc.nc") +print(f) + +# %% +# 3. Select the PM2.5 mass concentration and 2-metre temperature fields by +# index and print the description to show properties of all constructs: +pm25_field = f[2] +pm25_field.dump() + +# %% + +temp_field = f[3] +temp_field.dump() + +# %% +# 4. Convert the units to degree celsius for the temperature field: +temp_field.units = "degC" +temp_field.get_property("units") + +# %% +# 5. Digitize the PM2.5 mass concentration and 2-metre temperature fields. +# This step counts the number of values in each of the 10 equally sized bins +# spanning the range of the values. The ``'return_bins=True'`` argument makes +# sure that the calculated bins are also returned: +pm25_indices, pm25_bins = pm25_field.digitize(10, return_bins=True) +temp_indices, temp_bins = temp_field.digitize(10, return_bins=True) + +# %% +# 6. Create a joint histogram of the digitized fields: +joint_histogram = cf.histogram(pm25_indices, temp_indices) + +# %% +# 7. Get histogram data from the ``joint_histogram``: +histogram_data = joint_histogram.array + +# %% +# 8. Calculate bin centres for PM2.5 and temperature bins: +pm25_bin_centers = [(a + b) / 2 for a, b in pm25_bins] +temp_bin_centers = [(a + b) / 2 for a, b in temp_bins] + + +# %% +# 9. Create grids for PM2.5 and temperature bins using `numpy.meshgrid +# `_: +temp_grid, pm25_grid = np.meshgrid(temp_bin_centers, pm25_bin_centers) + +# %% +# 10. Plot the joint histogram using `matplotlib.pyplot.pcolormesh +# `_. Use `cf.Field.unique +# `_ to get +# the unique data array values and show the bin boundaries as ticks on the +# plot. ``'style='plain'`` in `matplotlib.axes.Axes.ticklabel_format +# `_ +# disables the scientific notation on the y-axis: +plt.pcolormesh(temp_grid, pm25_grid, histogram_data, cmap="viridis") +plt.xlabel("2-metre Temperature (degC)") +plt.ylabel("PM2.5 Mass Concentration (kg m**-3)") +plt.xticks(temp_bins.unique().array, rotation=45) +plt.yticks(pm25_bins.unique().array) +plt.colorbar(label="Frequency") +plt.title("Joint Histogram of PM2.5 and 2-metre Temperature", y=1.05) +plt.ticklabel_format(style="plain") +plt.show() diff --git a/recipes-docs/source/recipes-source/plot_10_recipe.py b/recipes-docs/source/recipes-source/plot_10_recipe.py new file mode 100644 index 0000000000..8c1a3f65dd --- /dev/null +++ b/recipes-docs/source/recipes-source/plot_10_recipe.py @@ -0,0 +1,96 @@ +""" +Calculating and plotting the relative vorticity +=============================================== + +Vorticity, the microscopic measure of rotation in a fluid, is a vector field +defined as the curl of velocity +`(James R. Holton, Gregory J. Hakim, An Introduction to Dynamic Meteorology, +2013, Elsevier : Academic Press p95-125) +`_. +In this recipe, we will be calculating and plotting the relative vorticity +from the wind components. +""" + +# %% +# 1. Import cf-python and cf-plot: + +import cfplot as cfp + +import cf + +# %% +# 2. Read the field constructs: +f = cf.read("~/recipes/ERA5_monthly_averaged_pressure_levels.nc") +print(f) + +# %% +# 3. Select wind components and look at their contents: +u = f.select_field("eastward_wind") +print(u) + +# %% + +v = f.select_field("northward_wind") +print(v) + +# %% +# 4. Create a date-time object for the required time period: +jan_2023 = cf.year(2023) & cf.month(1) + +# %% +# 5. The relative vorticity is calculated using `cf.curl_xy +# `_ and +# plotted using `cfplot.con `_. +# The ``with cf.relaxed_identities(True)`` context manager statement prevents +# the curl operation broadcasting across the two ``expver`` dimensions because +# it can't be certain that they are the same as they lack the standardised +# metadata. Setting +# ``cf.relaxed_identities(True)`` allows the ``long_name`` to be treated +# as standardised metadata. Since the horizontal coordinates are latitude and +# longitude, the +# `cf.curl_xy `_ +# function automatically accounts for the Earth's spherical geometry when +# calculating the spatial derivatives in the horizontal directions, and for this +# it requires the Earth's radius. In this case the radius is not stored in the +# wind fields, so must be provided by setting ``radius="earth"`` keyword +# parameter. While plotting, the relative vorticity is subspaced for January +# 2023 and one of the `experiment versions` using the dictionary unpacking +# operator (``**``) as there is an equal to sign in the identifier +# (``"long_name=expver"``): + +with cf.relaxed_identities(True): + rv = cf.curl_xy(u, v, radius="earth") + +cfp.con( + rv.subspace(T=jan_2023, **{"long_name=expver": 1}), + lines=False, + title="Relative Vorticity", +) + +# %% +# 6. Although the X axis is cyclic, it is not recognised as such, owing to the +# fact that the longitude coordinate bounds are missing. This results in +# discontinuities in the calculated vorticity field on the plot at the +# wrap-around location of 0 degrees east. The cyclicity could either be set on +# the field itself or just in the curl command by setting ``'x_wrap=True'`` +# while calculating the relative vorticity. Setting ``rv.units = "s-1"``, +# ensures that the units of the relative vorticity field are consistent with +# the calculation and the physical interpretation of the quantity: + +print(v.coordinate("X").has_bounds()) + +# %% + +with cf.relaxed_identities(True): + rv = cf.curl_xy(u, v, x_wrap=True, radius="earth") + +rv.units = "s-1" +print(rv) + +# %% + +cfp.con( + rv.subspace(T=jan_2023, **{"long_name=expver": 1}), + lines=False, + title="Relative Vorticity", +) diff --git a/recipes-docs/source/recipes-source/plot_11_recipe.py b/recipes-docs/source/recipes-source/plot_11_recipe.py new file mode 100644 index 0000000000..b147e51809 --- /dev/null +++ b/recipes-docs/source/recipes-source/plot_11_recipe.py @@ -0,0 +1,92 @@ +""" +Plotting the Warming Stripes +============================ + +In this recipe, we will plot the `Warming Stripes (Climate Stripes) +`_ created by +Professor Ed Hawkins at NCAS, University of Reading. Here we will use the +ensemble mean of the +`HadCRUT.5.0.1.0 analysis gridded data +`_ for +the same. + +""" + +# %% +# 1. Import cf-python and matplotlib.pyplot: + +import matplotlib.pyplot as plt + +import cf + +# %% +# 2. Read the field constructs: +temperature_data = cf.read( + "~/recipes/HadCRUT.5.0.1.0.analysis.anomalies.ensemble_mean.nc" +)[0] +print(temperature_data) + +# %% +# 3. Calculate the annual mean temperature anomalies. The ``'weights=True'`` +# argument is used to take the varying lengths of months into account which +# ensures that the calculated mean is more accurate: +annual_temperature = temperature_data.collapse( + "T: mean", weights=True, group=cf.Y() +) + +# %% +# 4. Select the data from 1850 to 2022: +period = annual_temperature.subspace(T=cf.year(cf.wi(1850, 2022))) + +# %% +# 5. Calculate the global average temperature for each year: +global_temperature = period.collapse("X: Y: mean") + +# %% +# 6. Get the global average temperature and squeeze it to remove the size 1 axis: +global_avg_temp = global_temperature.array.squeeze() + +# %% +# 7. Create a normalisation function that maps the interval from the minimum to +# the maximum temperature to the interval [0, 1] for colouring: +norm_global = plt.Normalize(global_avg_temp.min(), global_avg_temp.max()) + +# %% +# 8. Set the colormap instance: +cmap = plt.get_cmap("RdBu_r") + +# %% +# 9. Create the figure and the axes for the global plot. Loop over the selected +# years, plot a colored vertical stripe for each and remove the axes: +fig_global, ax_global = plt.subplots(figsize=(10, 2)) + +for i in range(global_avg_temp.shape[0]): + ax_global.axvspan( + xmin=i - 0.5, xmax=i + 0.5, color=cmap(norm_global(global_avg_temp[i])) + ) + +ax_global.axis("off") + +plt.show() + +# %% +# 10. For the regional warming stripes, steps 5 to 9 are repeated for the +# specific region. Here, we define the bounding box for UK by subspacing over +# a domain spanning 49.9 to 59.4 degrees north and -10.5 to 1.8 degrees east: +uk_temperature = period.subspace(X=cf.wi(-10.5, 1.8), Y=cf.wi(49.9, 59.4)) +uk_avg_temperature = uk_temperature.collapse("X: Y: mean") +uk_avg_temp = uk_avg_temperature.array.squeeze() +norm_uk = plt.Normalize(uk_avg_temp.min(), uk_avg_temp.max()) + +# %% + +fig_uk, ax_uk = plt.subplots(figsize=(10, 2)) + +for i in range(uk_avg_temp.shape[0]): + ax_uk.axvspan( + xmin=i - 0.5, xmax=i + 0.5, color=cmap(norm_uk(uk_avg_temp[i])) + ) + +ax_uk.axis("off") + +plt.show() diff --git a/recipes-docs/source/recipes-source/plot_12_recipe.py b/recipes-docs/source/recipes-source/plot_12_recipe.py new file mode 100644 index 0000000000..b09db0b29f --- /dev/null +++ b/recipes-docs/source/recipes-source/plot_12_recipe.py @@ -0,0 +1,117 @@ +""" +Using mask to plot Aerosol Optical Depth +======================================== + +In this recipe, we will make use of a +`masked array +`_ +to plot the `high-quality` retrieval of Aerosol Optical Depth (AOD) from all other +retrievals. + +""" + +# %% +# 1. Import cf-python, cf-plot and matplotlib.pyplot: + +import matplotlib.pyplot as plt +import cfplot as cfp + +import cf + +# %% +# 2. Read the field constructs: +fl = cf.read( + "~/recipes/JRR-AOD_v3r0_npp_s202012310752331_e202012310753573_c202100000000000.nc" +) +print(fl) + +# %% +# 3. Select AOD from the field list by identity and look at the contents: +aod = fl.select_field("long_name=AOT at 0.55 micron for both ocean and land") +print(aod) + +# %% +# 4. Select AOD retrieval quality by index and look at the quality flags: +quality = fl[13] +print(quality) + +# %% +# 5. Select latitude and longitude dimensions by identities, with two different +# techniques: +lon = aod.coordinate("long_name=Longitude") +lat = aod.coordinate("Y") + +# %% +# 6. Plot the AOD for all the retrievals using +# `cfplot.con `_. Here the argument +# ``'ptype'`` specifies the type of plot to use (latitude-longitude here) and +# the argument ``'lines=False'`` does not draw contour lines: +cfp.con(f=aod.array, x=lon.array, y=lat.array, ptype=1, lines=False) + +# %% +# 7. Create a mask for AOD based on the quality of the retrieval. The +# ``'__ne__'`` method is an implementation of the ``!=`` operator. It is used to +# create a mask where all the `high-quality` AOD points (with the flag 0) are +# marked as ``False``, and all the other data points (medium quality, low +# quality, or no retrieval) are marked as ``True``: +mask = quality.array.__ne__(0) + +# %% +# 8. Apply the mask to the AOD dataset. The ``'where'`` function takes the +# mask as an input and replaces all the values in the AOD dataset that +# correspond to ``True`` in the mask with a masked value using `cf.masked +# `_. +# In this case, all AOD values that are not of `high-quality` (since they were +# marked as ``True`` in the mask) are masked. This means that the ``high`` +# variable contains only the AOD data that was retrieved with `high-quality`: +high = aod.where(mask, cf.masked) + +# %% +# 9. Now plot both the AOD from `high-quality` retrieval and all other retrievals +# using `cfplot.con `_. Here: +# +# - `cfplot.gopen `_ is used to +# define the parts of the plot area, specifying that the figure should have +# 1 row and 2 columns, which is closed by +# `cfplot.gclose `_; +# - `plt.suptitle `_ +# is used to add a title for the whole figure; +# - the subplots for plotting are selected using +# `cfplot.gpos `_ after which +# `cfplot.mapset `_ is used to +# set the map limits and resolution for the subplots; +# - and as cf-plot stores the plot in a plot object with the name +# ``cfp.plotvars.plot``, country borders are added using normal +# `Cartopy operations `_ +# on the ``cfp.plotvars.mymap`` object: +import cartopy.feature as cfeature + +cfp.gopen(rows=1, columns=2, bottom=0.2) +plt.suptitle("AOD for both ocean and land", fontsize=20) +cfp.gpos(1) +cfp.mapset(resolution="50m", lonmin=68, lonmax=98, latmin=7, latmax=36) +cfp.con( + f=aod.array, + x=lon.array, + y=lat.array, + ptype=1, + lines=False, + title="All retrievals", + colorbar=None, +) +cfp.plotvars.mymap.add_feature(cfeature.BORDERS) +cfp.gpos(2) +cfp.mapset(resolution="50m", lonmin=68, lonmax=98, latmin=7, latmax=36) +cfp.con( + f=high.array, + x=lon.array, + y=lat.array, + ptype=1, + lines=False, + title="High quality retrieval", + colorbar_position=[0.1, 0.20, 0.8, 0.02], + colorbar_orientation="horizontal", + colorbar_title="AOD at 0.55 $\mu$", +) +cfp.plotvars.mymap.add_feature(cfeature.BORDERS) +cfp.gclose() diff --git a/recipes-docs/source/recipes-source/plot_13_recipe.py b/recipes-docs/source/recipes-source/plot_13_recipe.py new file mode 100644 index 0000000000..bf0398713e --- /dev/null +++ b/recipes-docs/source/recipes-source/plot_13_recipe.py @@ -0,0 +1,268 @@ +""" +Calculate and plot the Niño 3.4 Index +===================================== + +In this recipe, we will calculate and plot the sea surface temperature (SST) +anomaly in the Niño 3.4 region. According to `NCAR Climate Data Guide +`_, +the Niño 3.4 anomalies may be thought of as representing the average equatorial +SSTs across the Pacific from about the dateline to the South American coast. +The Niño 3.4 index typically uses a 5-month running mean, and El Niño or La +Niña events are defined when the Niño 3.4 SSTs exceed +/- 0.4 degrees Celsius for a +period of six months or more. + +""" + +# %% +# 1. Import cf-python and cf-plot, as well as some other libraries for use +# in next steps. + +import cartopy.crs as ccrs +import matplotlib.patches as mpatches + +import cfplot as cfp + +import cf + + +# %% +# 2. Read and select the SST by index and look at its contents: +sst = cf.read("~/recipes/ERA5_monthly_averaged_SST.nc")[0] +print(sst) + +# %% +# 3. Set the units from Kelvin to degrees Celsius: +sst.Units = cf.Units("degreesC") + +# %% +# 4. SST is subspaced for the Niño 3.4 region (5N-5S, 170W-120W) and as the +# dataset is using longitudes in 0-360 degrees East format, they are subtracted +# from 360 to convert them: +region = sst.subspace(X=cf.wi(360 - 170, 360 - 120), Y=cf.wi(-5, 5)) + +# %% +# 5. Plot the various Niño regions using cf-plot. Here: +# +# - `cfplot.gopen `_ is used to +# define the parts of the plot area, which is closed by +# `cfplot.gclose `_; +# - `cfplot.mapset `_ is used to +# set the map limits and projection; +# - `cfplot.setvars `_ is used to +# set various attributes of the plot, like setting the land colour to grey; +# - `cfplot.cscale `_ is used to +# choose one of the colour maps amongst many available; +# - `cfplot.con `_ plots contour data +# from the ``region`` subspace at a specific time with no contour lines and a +# title; +# - next, four Niño regions and labels are defined using +# `Matplotlib's Rectangle `_ +# and +# `Text `_ +# function with cf-plot plot object (``cfp.plotvars.plot``): + +cfp.gopen() +cfp.mapset(proj="cyl", lonmin=0, lonmax=360, latmin=-90, latmax=90) +cfp.setvars(land_color="grey") +cfp.cscale(scale="scale1") +cfp.con( + region.subspace(T=cf.dt(2022, 12, 1, 0, 0, 0, 0)), + lines=False, + title="Niño Index Regions", +) + +# Niño 3.4 region(5N-5S, 170W-120W): +rectangle = mpatches.Rectangle( + (-170, -5), + 50, + 10, + fill=False, + linewidth=1, + edgecolor="black", + transform=ccrs.PlateCarree(), +) +cfp.plotvars.mymap.add_patch(rectangle) +cfp.plotvars.mymap.text( + -145, + 7, + "3.4", + horizontalalignment="center", + fontsize=14, + weight="bold", + transform=ccrs.PlateCarree(), +) + +# Niño 1+2 region (0-10S, 90W-80W): +rectangle = mpatches.Rectangle( + (-90, 0), + 10, + 10, + hatch="**", + fill=False, + linewidth=1, + edgecolor="black", + alpha=0.3, + transform=ccrs.PlateCarree(), +) +cfp.plotvars.mymap.add_patch(rectangle) +cfp.plotvars.mymap.text( + -85, + 3, + "1+2", + horizontalalignment="center", + fontsize=8, + weight="bold", + transform=ccrs.PlateCarree(), +) + +# Niño 3 region (5N-5S, 150W-90W): +rectangle = mpatches.Rectangle( + (-150, -5), + 60, + 10, + hatch="xxx", + fill=False, + linewidth=1, + edgecolor="black", + alpha=0.3, + transform=ccrs.PlateCarree(), +) +cfp.plotvars.mymap.add_patch(rectangle) +cfp.plotvars.mymap.text( + -120, + -3, + "3", + horizontalalignment="center", + fontsize=14, + weight="bold", + transform=ccrs.PlateCarree(), +) + +# Niño 4 region (5N-5S, 160E-150W): +rectangle = mpatches.Rectangle( + (-200, -5), + 50, + 10, + hatch="oo", + fill=False, + linewidth=1, + edgecolor="black", + alpha=0.3, + transform=ccrs.PlateCarree(), +) +cfp.plotvars.mymap.add_patch(rectangle) +cfp.plotvars.mymap.text( + -175, + -3, + "4", + horizontalalignment="center", + fontsize=14, + weight="bold", + transform=ccrs.PlateCarree(), +) +cfp.gclose() + +# %% +# 6. Calculate the Niño 3.4 index and standardise it to create an anomaly index. +# The `collapse `_ +# method is used to calculate the mean over the longitude (X) and latitude (Y) +# dimensions: +nino34_index = region.collapse("X: Y: mean") + +# %% +# 7. The result, ``nino34_index``, represents the average SST in the defined +# Niño 3.4 region for each time step. In the variable ``base_period``, +# ``nino34_index`` is subset to only include data from the years 1961 to 1990. +# This period is often used as a reference period for calculating anomalies. +# The variables ``climatology`` and ``std_dev`` include the mean and the +# standard deviation over the time (T) dimension of the ``base_period`` data +# respectively: +base_period = nino34_index.subspace(T=cf.year(cf.wi(1961, 1990))) +climatology = base_period.collapse("T: mean") +std_dev = base_period.collapse("T: sd") + +# %% +# 8. The line for variable ``nino34_anomaly`` calculates the standardised +# anomaly for each time step in the ``nino34_index`` data. It subtracts the +# ``climatology`` from the ``nino34_index`` and then divides by the ``std_dev``. +# The resulting ``nino34_anomaly`` data represents how much the SST in the Niño +# 3.4 region deviates from the 1961-1990 average, in units of standard +# deviations. This is a common way to quantify climate anomalies like El Niño +# and La Niña events: +nino34_anomaly = (nino34_index - climatology) / std_dev + +# %% +# 9. A moving average of the ``nino34_anomaly`` along the time axis, with a +# window size of 5 (i.e. an approximately 5-month moving average) is calculated +# using the +# `moving_window `_ +# method. The ``mode='nearest'`` parameter is used to specify how to pad the +# data outside of the time range. The resulting ``nino34_rolling`` variable +# represents a smoothed version of the ``nino34_anomaly`` data. It removes +# short-term fluctuations and highlights longer-term trends or cycles: +nino34_rolling = nino34_anomaly.moving_window( + method="mean", window_size=5, axis="T", mode="nearest" +) + +# %% +# 10. Define El Niño and La Niña events by creating Boolean masks to identify +# El Niño and La Niña events. Now plot SST anomalies in the Niño 3.4 region over +# time using cf-plot. Here: +# +# - `cfplot.gset `_ sets the limits +# of the x-axis (years from 1940 to 2022) and y-axis (anomalies from -3 +# degrees C to 3 degrees C) for the plot; +# - `cfplot.gopen `_ is used to +# define the parts of the plot area, which is closed by +# `cfplot.gclose `_; +# - `cfplot.lineplot `_ plots +# the rolling Niño 3.4 index over time; +# - a zero line and also horizontal dashed lines are drawn for El Niño and +# La Niña thresholds using +# `Matplotlib's axhline `_ +# with cf-plot plot object (``cfp.plotvars.plot``); +# - `fill_between `_ +# from Matplotlib is used with cf-plot plot object (``cfp.plotvars.plot``) +# to fill the area between the Niño 3.4 index and the El Niño/La Niña +# thresholds; +# - similarly, +# `cfplot.plotvars.plot.legend `_ +# is used to add a legend in the end: +elnino = nino34_rolling >= 0.4 +lanina = nino34_rolling <= -0.4 + +cfp.gset(xmin="1940-1-1", xmax="2022-12-31", ymin=-3, ymax=3) + +cfp.gopen(figsize=(10, 6)) +cfp.lineplot( + nino34_rolling, + color="black", + title="SST Anomaly in Niño 3.4 Region (5N-5S, 120-170W)", + ylabel="Temperature anomaly ($\degree C$)", + xlabel="Year", +) +cfp.plotvars.plot.axhline( + 0.4, color="red", linestyle="--", label="El Niño Threshold" +) +cfp.plotvars.plot.axhline( + -0.4, color="blue", linestyle="--", label="La Niña Threshold" +) +cfp.plotvars.plot.axhline(0, color="black", linestyle="-", linewidth=1) +cfp.plotvars.plot.fill_between( + nino34_rolling.coordinate("T").array, + 0.4, + nino34_rolling.array.squeeze(), + where=elnino.squeeze(), + color="red", + alpha=0.3, +) +cfp.plotvars.plot.fill_between( + nino34_rolling.coordinate("T").array, + -0.4, + nino34_rolling.array.squeeze(), + where=lanina.squeeze(), + color="blue", + alpha=0.3, +) +cfp.plotvars.plot.legend(frameon=False, loc="lower center", ncol=2) +cfp.gclose() diff --git a/recipes-docs/source/recipes-source/plot_14_recipe.py b/recipes-docs/source/recipes-source/plot_14_recipe.py new file mode 100644 index 0000000000..957a3f38d7 --- /dev/null +++ b/recipes-docs/source/recipes-source/plot_14_recipe.py @@ -0,0 +1,181 @@ +""" +Overlay Geopotential height contours over Temperature anomalies +=============================================================== + +In this recipe, we will overlay Geopotential height contours over Temperature +anomalies to help analyse meteorological conditions during July 2018, +specifically focusing on the significant concurrent extreme events that occurred +during the 2018 boreal spring/summer season in the Northern Hemisphere. + +""" + +# %% +# 1. Import cf-python and cf-plot: +import cfplot as cfp + +import cf + +# %% +# 2. Read and select the 200 hPa geopotential by index and look at its contents: +gp = cf.read("~/recipes/ERA5_monthly_averaged_z200.nc")[0] +print(gp) + +# %% +# 3. Convert the geopotential data to geopotential height by dividing it by the +# acceleration due to gravity (approximated as 9.81 :math:`m \cdot {s}^{-2}`): +gph = gp / 9.81 + +# %% +# 4. Subset the geopotential height to extract data specifically for July 2018, +# a significant month due to heat extremes and heavy rainfall: +gph_july = gph.subspace(T=cf.month(7) & cf.year(2018)).squeeze() + +# %% +# 5. Plot contour lines of this geopotential height for July 2018. Here: +# +# - `cfplot.gopen `_ is used to +# define the parts of the plot area, which is closed by +# `cfplot.gclose `_; +# - `cfplot.mapset `_ is used to +# set the map projection to North Polar Stereographic; +# - `cfplot.setvars `_ is used to +# set various attributes of the plot, like setting the thickness of the lines +# that represent continents; +# - `cfplot.con `_ plots the contour +# lines representing the 200 hPa geopotential height values without filling +# between the contour lines (``fill=False``) and no colour bar +# (``colorbar=False``); +# - `cfplot.levs `_ is used to +# specify two contour levels, 12000 and 12300 m, corresponding to the +# approximate polar-front jet and subtropical jet respectively; +# - `cfplot.con `_ is again used to +# plot the contour lines for polar-front jet and subtropical jet with a +# thicker line width; +# - `cfp.plotvars.mymap.stock_img() `_ +# then finally visualises the Earth's surface in cf-plot's +# ``cfp.plotvars.mymap`` plot object: +cfp.gopen() +cfp.mapset(proj="npstere") +cfp.setvars(continent_thickness=0.5) + +cfp.con( + f=gph_july, + fill=False, + lines=True, + line_labels=False, + colors="black", + linewidths=1, + colorbar=False, +) + +cfp.levs(manual=[12000, 12300]) +cfp.con( + f=gph_july, + fill=False, + lines=True, + colors="black", + linewidths=3.0, + colorbar=False, +) + +cfp.plotvars.mymap.stock_img() +cfp.gclose() + +# %% +# 6. Read and select the 2-metre temperature by index and look at its contents: +t2m = cf.read("~/recipes/ERA5_monthly_averaged_t2m.nc")[0] +print(t2m) + +# %% +# 7. Set the units from Kelvin to degrees Celsius: +t2m.Units = cf.Units("degreesC") + +# %% +# 8. Extract a subset for July across the years for ``t2m``: +t2m_july = t2m.subspace(T=cf.month(7)) + +# %% +# 9. The 2-meter temperature climatology is then calculated for the month of +# July over the period from 1981 to 2010, which provides a baseline against +# which anomalies in later years are compared: +t2m_july_climatology = t2m_july.subspace( + T=cf.year(cf.wi(1981, 2010)) +).collapse("T: mean") + +# %% +# 10. Calculate the temperature anomaly for the month of July in the year 2018 +# relative to the climatological baseline (``t2m_july_climatology``). This +# indicates how much the temperatures for that month in that year deviated from +# the long-term average for July across the 1981-2010 period: +t2m_july_anomaly_2018 = ( + t2m_july.subspace(T=cf.year(2018)).squeeze() - t2m_july_climatology +) + +# %% +# 11. +# The July 2018 season experienced extreme heat in many parts of the Northern +# Hemisphere. This period's extreme events were related to unusual +# meteorological conditions, particularly abnormalities in the jet stream. To +# provide an insight into the atmospheric conditions, the temperature anomalies +# and the geopotential height contours are plotted using cf-plot. Here: +# +# - `cfplot.gopen `_ is used to +# define the parts of the plot area, which is closed by +# `cfplot.gclose `_; +# - `cfplot.mapset `_ is used to +# set the map projection to Robinson; +# - `cfplot.setvars `_ is used to +# set various attributes of the plot, like setting the thickness of the lines +# that represent continents and master title properties; +# - `cfplot.levs `_ is used to +# specify the contour levels for temperature anomalies, starting from -2 to 2 +# with an interval of 0.5; +# - `cfplot.cscale `_ is used to +# choose one of the colour maps amongst many available; +# - `cfplot.con `_ plots contour fill +# of temperature anomalies without contour lines (``lines=False``); +# - `cfplot.levs() `_ is used to +# reset contour levels to default after which the steps to plot the contour +# lines representing the 200 hPa geopotential height values, the approximate +# polar-front jet and subtropical jet from Step 5 are repeated: +cfp.gopen() +cfp.mapset(proj="robin") +cfp.setvars( + continent_thickness=0.5, + master_title="July 2018", + master_title_fontsize=22, + master_title_location=[0.53, 0.83], +) + +cfp.levs(min=-2, max=2, step=0.5) +cfp.cscale("temp_19lev") +cfp.con( + f=t2m_july_anomaly_2018, + lines=False, + colorbar_title="Temperature anomaly relative to 1981-2010 ($\degree C$)", + colorbar_fontsize=13, + colorbar_thick=0.04, +) + +cfp.levs() +cfp.con( + f=gph_july, + fill=False, + lines=True, + line_labels=False, + colors="black", + linewidths=1, + colorbar=False, +) + +cfp.levs(manual=[12000, 12300]) +cfp.con( + f=gph_july, + fill=False, + lines=True, + colors="black", + linewidths=3.0, + colorbar=False, +) + +cfp.gclose() diff --git a/recipes-docs/source/recipes-source/plot_15_recipe.py b/recipes-docs/source/recipes-source/plot_15_recipe.py new file mode 100644 index 0000000000..d20f622826 --- /dev/null +++ b/recipes-docs/source/recipes-source/plot_15_recipe.py @@ -0,0 +1,163 @@ +""" +Resampling Land Use Flags to a Coarser Grid +=========================================== + +In this recipe, we will compare the land use distribution in different countries +using a land use data file and visualize the data as a histogram. This will help +to understand the proportion of different land use categories in each country. + +The land use data is initially available at a high spatial resolution of +approximately 100 m, with several flags defined with numbers representing the +type of land use. Regridding the data to a coarser resolution of approximately +25 km would incorrectly represent the flags on the new grids. + +To avoid this, we will resample the data to the coarser resolution by +aggregating the data within predefined spatial regions or bins. This approach +will give a dataset where each 25 km grid cell contains a histogram of land use +flags, as determined by the original 100 m resolution data. It retains the +original spatial extent of the data while reducing its spatial complexity. +Regridding, on the other hand, involves interpolating the data onto a new grid, +which can introduce artefacts and distortions in the data. + +""" + +import cartopy.io.shapereader as shpreader +import matplotlib.pyplot as plt +import numpy as np + +# %% +# 1. Import the required libraries. We will use Cartopy's ``shapereader`` to +# work with shapefiles that define country boundaries: +import cf + +# %% +# 2. Read and select land use data by index and see properties of +# all constructs: +f = cf.read("~/recipes/output.tif.nc")[0] +f.dump() + + +# %% +# 3. Define a function to extract data for a specific country: +# +# - The ``extract_data`` function is defined to extract land use data for a +# specific country, specified by the ``country_name`` parameter. +# - It uses the `Natural Earth `_ +# shapefile to get the bounding coordinates of the selected country. +# - The `shpreader.natural_earth `_ +# function is called to access the Natural +# Earth shapefile of country boundaries with a resolution of 10 m. +# - The `shpreader.Reader `_ +# function reads the shapefile, and the selected country's record is retrieved +# by filtering the records based on the ``NAME_LONG`` attribute. +# - The bounding coordinates are extracted using the ``bounds`` attribute of the +# selected country record. +# - The land use data file is then read and subset using these bounding +# coordinates with the help of the ``subspace`` function. The subset data is +# stored in the ``f`` variable. + + +def extract_data(country_name): + shpfilename = shpreader.natural_earth( + resolution="10m", category="cultural", name="admin_0_countries" + ) + reader = shpreader.Reader(shpfilename) + country = [ + country + for country in reader.records() + if country.attributes["NAME_LONG"] == country_name + ][0] + lon_min, lat_min, lon_max, lat_max = country.bounds + + f = cf.read("~/recipes/output.tif.nc")[0] + f = f.subspace(X=cf.wi(lon_min, lon_max), Y=cf.wi(lat_min, lat_max)) + + return f + + +# %% +# 4. Define a function to plot a histogram of land use distribution for a +# specific country: +# +# - The `digitize `_ +# function of the ``cf.Field`` object is called to convert the land use data +# into indices of bins. It takes an array of bins (defined by +# the `np.linspace `_ function) +# and the ``return_bins=True`` parameter, which returns the actual bin values +# along with the digitized data. +# - The `np.linspace `_ +# function is used to create an array of evenly spaced bin edges from 0 to 50, +# with 51 total values. This creates bins of width 1. +# - The ``digitized`` variable contains the bin indices for each data point, +# while the bins variable contains the actual bin values. +# - The `cf.histogram `_ +# function is called on the digitized data to create a histogram. This +# function returns a field object with the histogram data. +# - The `squeeze `_ +# function applied to the histogram ``array`` extracts the histogram data as a NumPy +# array and removes any single dimensions. +# - The ``total_valid_sub_cells`` variable calculates the total number of valid +# subcells (non-missing data points) by summing the histogram data. +# - The last element of the bin_counts array is removed with slicing +# (``bin_counts[:-1]``) to match the length of the ``bin_indices`` array. +# - The ``percentages`` variable calculates the percentage of each bin by +# dividing the ``bin_counts`` by the ``total_valid_sub_cells`` and multiplying +# by 100. +# - The ``bin_indices`` variable calculates the centre of each bin by averaging +# the bin edges. This is done by adding the ``bins.array[:-1, 0]`` and +# ``bins.array[1:, 0]`` arrays and dividing by 2. +# - The ``ax.bar`` function is called to plot the histogram as a bar chart on +# the provided axis. The x-axis values are given by the ``bin_indices`` array, +# and the y-axis values are given by the ``percentages`` array. +# - The title, x-axis label, y-axis label, and axis limits are set based on the +# input parameters. + + +def plot_histogram(field, ax, label, ylim, xlim): + digitized, bins = field.digitize(np.linspace(0, 50, 51), return_bins=True) + + h = cf.histogram(digitized) + bin_counts = h.array.squeeze() + + total_valid_sub_cells = bin_counts.sum() + + bin_counts = bin_counts[:-1] + + percentages = bin_counts / total_valid_sub_cells * 100 + + bin_indices = (bins.array[:-1, 0] + bins.array[1:, 0]) / 2 + + ax.bar(bin_indices, percentages, label=label) + ax.set_title(label) + ax.set_xlabel("Land Use Flag") + ax.set_ylabel("Percentage") + ax.set_ylim(ylim) + ax.set_xlim(xlim) + + +# %% +# 5. Define the countries of interest: +countries = ["Ireland", "Belgium", "Switzerland"] + +# %% +# 6. Set up the figure and axes for plotting the histograms: +# +# - The ``plt.subplots`` function is called to set up a figure with three +# subplots, with a figure size of 8 inches by 10 inches. +# - A loop iterates over each country in the countries list and for each +# country, the ``extract_data`` function is called to extract its land use +# data. +# - The ``plot_histogram`` function is then called to plot the histogram of land +# use distribution on the corresponding subplot. +# - The ``plt.tight_layout`` function is called to ensure that the subplots are +# properly spaced within the figure and finally, the ``plt.show`` function +# displays the figure with the histograms. +fig, axs = plt.subplots(3, 1, figsize=(8, 10)) + +for i, country in enumerate(countries): + ax = axs[i] + data = extract_data(country) + plot_histogram(data, ax, label=country, ylim=(0, 50), xlim=(0, 50)) + +plt.tight_layout() +plt.show() diff --git a/recipes-docs/source/recipes-source/plot_16_recipe.py b/recipes-docs/source/recipes-source/plot_16_recipe.py new file mode 100644 index 0000000000..00a7f2e1d0 --- /dev/null +++ b/recipes-docs/source/recipes-source/plot_16_recipe.py @@ -0,0 +1,57 @@ +""" +Plotting contour subplots with different projections +==================================================== + +In this recipe, we will plot the same data using different projections +as subplots to illustrate visually some available possibilities. + +""" + +# %% +# 1. Import cf-python and cf-plot: + +import cfplot as cfp + +import cf + +# %% +# 2. Read the field in: +f = cf.read("~/recipes/ggap.nc")[0] + +# %% +# 3. List the projection types to use. Here we are using +# Cylindrical/Default, North Pole Stereographic, South Pole Stereographic, +# Mollweide, Mercator and Robinson. However there are several other choices +# possible, see: +# https://ncas-cms.github.io/cf-plot/build/user_guide.html#appendixc. Our +# chosen list is: +projtypes = ["cyl", "npstere", "spstere", "moll", "merc", "robin"] + +# %% +# 4. Create the file with subplots. If changing the number of subplots, +# ensure the number of rows * number of columns = the number of projections. +# Here we are doing 6 projections so 2 x 3 is fine. Then loop through the +# list of projection types and plot each as a sub-plot: +cfp.gopen(rows=2, columns=3, bottom=0.2) +for i, proj in enumerate(projtypes): + # gpos has 1 added to the index because it takes 1 as its first value + cfp.gpos(i + 1) + cfp.mapset(proj=proj) + + # For the final plot only, add a colour bar to cover all the sub-plots + if i == len(projtypes) - 1: + cfp.con( + f.subspace(pressure=850), + lines=False, + title=proj, + colorbar_position=[0.1, 0.1, 0.8, 0.02], + colorbar_orientation="horizontal", + ) + else: + cfp.con( + f.subspace(pressure=850), + lines=False, + title=proj, + colorbar=False, + ) +cfp.gclose() diff --git a/recipes-docs/source/recipes-source/plot_17_recipe.py b/recipes-docs/source/recipes-source/plot_17_recipe.py new file mode 100644 index 0000000000..c94769e2ba --- /dev/null +++ b/recipes-docs/source/recipes-source/plot_17_recipe.py @@ -0,0 +1,109 @@ +""" +Plotting contour subplots with different colour maps/scales +=========================================================== + +In this recipe, we will plot the same data with different colour maps from +three categories in separate subplots to illustrate the importance of +choosing a suitable one for given data. To avoid unintended bias and +misrepresentation, or lack of accessibility, a careful choice must be made. +""" + +# %% +# 1. Import cf-python and cf-plot: + +import matplotlib.pyplot as plt +import cfplot as cfp + +import cf + +# %% +# 2. Read the field in: +f = cf.read("~/recipes/ggap.nc")[0] + +# %% +# 3. Choose a set of predefined colour scales to view. These can be chosen +# from the selection at: +# https://ncas-cms.github.io/cf-plot/build/colour_scales.html or you +# can define your own, see: +# https://ncas-cms.github.io/cf-plot/build/colour_scales.html#user-defined-colour-scales. +# Here we take three colour scales each from three different general +# categories, to showcase some differences in representation. +# Note colour scale levels can be adjusted using 'cscale' and keywords such as: +# cfp.cscale(, ncols=16, below=2, above=14) + +# %% +# a. Perceptually uniform colour scales, with no zero value, see: +# https://ncas-cms.github.io/cf-plot/build/colour_scales.html#perceptually-uniform-colour-scales. +colour_scales_pu = ["viridis", "magma", "plasma"] + +# %% +# b. NCAR Command Language colour scales enhanced to help with colour +# blindness, see: +# https://ncas-cms.github.io/cf-plot/build/colour_scales.html#ncar-command-language-enhanced-to-help-with-colour-blindness. +# These colour maps are better for accessibility. +colour_scales_ncl = ["posneg_1", "GreenMagenta16", "StepSeq25"] + +# %% +# c. Orography/bathymetry colour scales, see: +# https://ncas-cms.github.io/cf-plot/build/colour_scales.html#orography-bathymetry-colour-scales. +# These are used to show the shape/contour of land masses, but bear in mind the +# data we show here is air temperature so doesn't represent this and +# therefore it is not a good choice in this case: +colour_scales_ob = ["wiki_1_0_2", "wiki_2_0", "wiki_2_0_reduced"] + + +# %% +# 4. We plot each category of colourmap in a given columns of the subplot, +# but given the 'gpos' function positions subplots from left to right, row by +# row from the top, we need to interleave the values in a list. We can use +# zip to do this: +colour_scales_columns = [ + cscale + for category in zip(colour_scales_pu, colour_scales_ncl, colour_scales_ob) + for cscale in category +] + + +# %% +# 5. Create the figure and give it an overall title. Ensure the +# number of rows * number of columns = number of colour scales. +# Then we loop through all the different colour maps defined and plot +# as subplots, with each category in the same column, labelling each column +# with the colour scale category: +cfp.gopen(rows=3, columns=3, bottom=0.1, top=0.85) +plt.suptitle( + ( + "Air temperature (K) at 850 mbar pressure shown in different " + "categories of colour scale" + ), + fontsize=18, +) +for i, colour_scale in enumerate(colour_scales_columns): + cfp.gpos(i + 1) + cfp.cscale(colour_scale, ncols=15) + + # For the topmost plots, label the column with the colour scale category + # using the 'title' argument, otherwise don't add a title. + # Ensure the order the titles are written in corresponds to the + # order unzipped in step 4, so the columns match up correctly. + if i == 0: + set_title = "Perceptually uniform\ncolour maps" + elif i == 1: + set_title = "NCL colour maps enhanced to \nhelp with colour blindness" + elif i == 2: + set_title = "Orography/bathymetry\ncolour maps" + else: + set_title = "" + + cfp.con( + f.subspace(pressure=850), + title=set_title, + lines=False, + axes=False, + colorbar_drawedges=False, + colorbar_title=f"Shown in '{colour_scale}'", + colorbar_fraction=0.04, + colorbar_thick=0.02, + colorbar_fontsize=11, + ) +cfp.gclose() diff --git a/recipes-docs/source/recipes-source/plot_18_recipe.py b/recipes-docs/source/recipes-source/plot_18_recipe.py new file mode 100644 index 0000000000..f0eae36e35 --- /dev/null +++ b/recipes-docs/source/recipes-source/plot_18_recipe.py @@ -0,0 +1,143 @@ +""" +Calculating the Pearson correlation coefficient between datasets +================================================================ + +In this recipe, we will take two datasets, one for an independent variable +(in this example elevation) and one for a dependent variable (snow +cover over a particular day), regrid them to the same resolution then +calculate the correlation coefficient, to get a measure of the relationship +between them. + +""" + +# %% +# 1. Import cf-python, cf-plot and other required packages: +import matplotlib.pyplot as plt +import scipy.stats.mstats as mstats +import cfplot as cfp + +import cf + + +# %% +# 2. Read the data in and unpack the Fields from FieldLists using indexing. +# In our example We are investigating the influence of the land height on +# the snow cover extent, so snow cover is the dependent variable. The snow +# cover data is the +# 'Snow Cover Extent 2017-present (raster 500 m), Europe, daily – version 1' +# sourced from the Copernicus Land Monitoring Service which is described at: +# https://land.copernicus.eu/en/products/snow/snow-cover-extent-europe-v1-0-500m +# and the elevation data is the 'NOAA NGDC GLOBE topo: elevation data' dataset +# which can be sourced from the IRI Data Library, or details found, at: +# http://iridl.ldeo.columbia.edu/SOURCES/.NOAA/.NGDC/.GLOBE/.topo/index.html. +orog = cf.read("~/recipes/1km_elevation.nc")[0] +snow = cf.read("~/recipes/snowcover")[0] + +# %% +# 3. Choose the day of pre-aggregated snow cover to investigate. We will +# take the first datetime element corresponding to the first day from the +# datasets, 1st January 2024, but by changing the indexing you can explore +# other days by changing the index. We also get the string corresponding to +# the date, to reference later: +snow_day = snow[0] +snow_day_dt = snow_day.coordinate("time")[0].data +snow_day_daystring = f"{snow_day_dt.datetime_as_string[0].split(' ')[0]}" + +# %% +# 4. Choose the region to consider to compare the relationship across, +# which must be defined across both datasets, though not necessarily on the +# same grid since we regrid to the same grid next and subspace to the same +# area for both datasets ready for comparison in the next steps. By changing +# the latitude and longitude points in the tuple below, you can change the +# area that is used: +region_in_mid_uk = ((-3.0, -1.0), (52.0, 55.0)) +sub_orog = orog.subspace( + longitude=cf.wi(*region_in_mid_uk[0]), latitude=cf.wi(*region_in_mid_uk[1]) +) +sub_snow = snow_day.subspace( + longitude=cf.wi(*region_in_mid_uk[0]), latitude=cf.wi(*region_in_mid_uk[1]) +) + +# %% +# 5. Ensure data quality, since the standard name here corresponds to a +# unitless fraction, but the values are in the tens, so we need to +# normalise these to all lie between 0 and 1 and change the units +# appropriately: +sub_snow = (sub_snow - sub_snow.minimum()) / (sub_snow.range()) +sub_snow.override_units("1", inplace=True) + +# %% +# 6. Regrid the data so that they lie on the same grid and therefore each +# array structure has values with corresponding geospatial points that +# can be statistically compared. Here the elevation field is regridded to the +# snow field since the snow is higher-resolution, but the other way round is +# possible by switching the field order: +regridded_orog = sub_orog.regrids(sub_snow, method="linear") + +# %% +# 7. Squeeze the snow data to remove the size 1 axes so we have arrays of +# the same dimensions for each of the two fields to compare: +sub_snow = sub_snow.squeeze() + +# %% +# 8. Finally, perform the statistical calculation by using the SciPy method +# to find the Pearson correlation coefficient for the two arrays now they are +# in comparable form. Note we need to use 'scipy.stats.mstats' and not +# 'scipy.stats' for the 'pearsonr' method, to account for masked +# data in the array(s) properly: +coefficient = mstats.pearsonr(regridded_orog.array, sub_snow.array) +print(f"The Pearson correlation coefficient is: {coefficient}") + +# %% +# 9. Make a final plot showing the two arrays side-by-side and quoting the +# determined Pearson correlation coefficient to illustrate the relationship +# and its strength visually. We use 'gpos' to position the plots in two +# columns and apply some specific axes ticks and labels for clarity. +cfp.gopen( + rows=1, + columns=2, + top=0.85, + user_position=True, +) + +# Joint configuration of the plots, including adding an overall title +plt.suptitle( + ( + "Snow cover compared to elevation for the same area of the UK " + f"aggregated across\n day {snow_day_daystring} with correlation " + "coefficient (on the same grid) of " + f"{coefficient.statistic:.4g} (4 s.f.)" + ), + fontsize=17, +) +cfp.mapset(resolution="10m") +cfp.setvars(ocean_color="white", lake_color="white") +label_info = { + "xticklabels": ("3W", "2W", "1W"), + "yticklabels": ("52N", "53N", "54N", "55N"), + "xticks": (-3, -2, -1), + "yticks": (52, 53, 54, 55), +} + +# Plot the two contour plots as columns +cfp.gpos(1) +cfp.cscale("wiki_2_0_reduced", ncols=11) +cfp.con( + regridded_orog, + lines=False, + title="Elevation (from 1km-resolution orography)", + colorbar_drawedges=False, + **label_info, +) +cfp.gpos(2) +# Don't add extentions on the colourbar since it can only be 0 to 1 inclusive +cfp.levs(min=0, max=1, step=0.1, extend="neither") +cfp.cscale("precip_11lev", ncols=11, reverse=1) +cfp.con( + sub_snow, + lines=False, + title="Snow cover extent (from satellite imagery)", + colorbar_drawedges=False, + **label_info, +) +cfp.gclose() diff --git a/recipes-docs/source/recipes-source/plot_19_recipe.py b/recipes-docs/source/recipes-source/plot_19_recipe.py new file mode 100644 index 0000000000..02d493dc21 --- /dev/null +++ b/recipes-docs/source/recipes-source/plot_19_recipe.py @@ -0,0 +1,108 @@ +""" +Plotting per-season trends in global sea surface temperature extrema +==================================================================== + +In this recipe we find the area-based extrema of global sea surface +temperature per month and, because it is very difficult to +interpret for trends when in a monthly form, we calculate and plot +on top of this the mean across each season for both the minima and the +maxima. +""" + +# %% +# 1. Import cf-python, cf-plot and other required packages: +import matplotlib.pyplot as plt +import cfplot as cfp + +import cf + +# %% +# 2. Read the dataset in extract the SST Field from the FieldList: +f = cf.read("~/recipes/ERA5_monthly_averaged_SST.nc") +sst = f[0] # this gives the sea surface temperature (SST) + +# %% +# 3. Collapse the SST data by area extrema (extrema over spatial dimensions): +am_max = sst.collapse("area: maximum") # equivalent to "X Y: maximum" +am_min = sst.collapse("area: minimum") # equivalent to "X Y: minimum" + +# %% +# 4. Reduce all timeseries down to just 1980+ since there are some data +# quality issues before 1970 and also this window is about perfect size +# for viewing the trends without the line plot becoming too cluttered: +am_max = am_max.subspace(T=cf.ge(cf.dt("1980-01-01"))) +am_min = am_min.subspace(T=cf.ge(cf.dt("1980-01-01"))) + +# %% +# 5. Create a mapping which provides the queries we need to collapse on +# the four seasons, along with our description of them, as a value, with +# the key of the string encoding the colour we want to plot these +# trend lines in. This structure will be iterated over to make our plot: +colours_seasons_mapping = { + "red": (cf.mam(), "Mean across MAM: March, April and May"), + "blue": (cf.jja(), "Mean across JJA: June, July and August"), + "green": (cf.son(), "Mean across SON: September, October and November"), + "purple": (cf.djf(), "Mean across DJF: December, January and February"), +} + +# %% +# 6. Create and open the plot file. Put maxima subplot at top since these +# values are higher, given increasing x axis. +# Note we set limits manually with 'gset' only to +# allow space so the legend doesn't overlap the data, which isn't +# possible purely from positioning it anywhere within the default plot. +# Otherwise cf-plot handles this for us. To plot the per-season means +# of the maxima, we loop through the season query mapping and do a +# "T: mean" collapse setting the season as the grouping: +cfp.gopen( + rows=2, + columns=1, + bottom=0.1, + top=0.85, +) +cfp.gpos(1) +cfp.gset(xmin="1980-01-01", xmax="2022-12-01", ymin=304, ymax=312) +for colour, season_query in colours_seasons_mapping.items(): + query_on_season, season_description = season_query + am_max_collapse = am_max.collapse("T: mean", group=query_on_season) + cfp.lineplot( + am_max_collapse, + color=colour, + markeredgecolor=colour, + marker="o", + label=season_description, + title="Maxima per month or season", + ) +cfp.lineplot( + am_max, + color="grey", + xlabel="", + label="All months", +) +# Create and add minima subplot below the maxima one. Just like for the +# maxima case, we plot per-season means by looping through the season query +# mapping and doing a "T: mean" collapse setting the season as the grouping +cfp.gpos(2) +cfp.gset(xmin="1980-01-01", xmax="2022-12-01", ymin=269, ymax=272) +for colour, season_query in colours_seasons_mapping.items(): + query_on_season, season_description = season_query + am_min_collapse = am_min.collapse("T: mean", group=query_on_season) + cfp.lineplot( + am_min_collapse, + color=colour, + markeredgecolor=colour, + marker="o", + xlabel="", + title="Minima per month or season", + ) +cfp.lineplot( + am_min, + color="grey", +) +# Add an overall title to the plot and close the file to save it +plt.suptitle( + "Global mean sea surface temperature (SST) monthly\nminima and maxima " + "showing seasonal means of these extrema", + fontsize=18, +) +cfp.gclose() diff --git a/recipes-docs/source/recipes-source/plot_20_recipe.py b/recipes-docs/source/recipes-source/plot_20_recipe.py new file mode 100644 index 0000000000..95e02b01d7 --- /dev/null +++ b/recipes-docs/source/recipes-source/plot_20_recipe.py @@ -0,0 +1,97 @@ +""" +Calculating and plotting the divergence of sea currents +======================================================= + +In this recipe, we will calculate the divergence of depth-averaged +currents in the Irish Sea, then plot the divergence as a contour +fill plot underneath the vectors themselves in the form of a vector plot. +""" + +# %% +# 1. Import cf-python and cf-plot: +import cfplot as cfp + +import cf + +# %% +# 2. Read the fields in. This dataset consists of depth-averaged eastward and +# northward current components plus the sea surface height above sea level and +# is a gridded dataset, with grid resolution of 1.85 km, covering the entire +# Irish Sea area. It was found via the CEDA Archive at the location of: +# https://catalogue.ceda.ac.uk/uuid/1b89e025eedd49e8976ee0721ec6e9b5, with +# DOI of https://dx.doi.org/10.5285/031e7ca1-9710-280d-e063-6c86abc014a0: +f = cf.read("~/recipes/POLCOMS_WAM_ZUV_01_16012006.nc") + +# %% +# 3. Get the separate vector components, which are stored as separate fields. +# The first, 'u', corresponds to the eastward component and the second, 'v', +# the northward component: +u = f[0] +v = f[1] + +# %% +# 4. Squeeze the fields to remove the size 1 axes in each case: +u = u.squeeze() +v = v.squeeze() + +# %% +# 5. Consider the currents at a set point in time. To do this we +# select one of the 720 datetime sample points in the fields to +# investigate, in this case by subspacing to pick out a particular +# datetime value we saw within the time coordinate data of the field (but +# you could also use indexing or filtering to select a specific value). +# Once we subspace to one datetime, we squeeze out the size 1 time axis +# in each case: +chosen_time = "2006-01-15 23:30:00" # 720 choices to pick from, try this one! +u_1 = u.subspace(T=cf.dt(chosen_time)) +v_1 = v.subspace(T=cf.dt(chosen_time)) +u_1 = u_1.squeeze() +v_1 = v_1.squeeze() + +# %% +# 6. +# When inspecting the u and v fields using cf inspection methods such as +# from print(u_1.data) and u_1.data.dump(), for example, we can see that there are +# lots of -9999 values in their data array, apparently used as a +# fill/placeholder value, including to indicate undefined data over the land. +# In order for these to not skew the data and dominate the plot, we need +# to mask values matching this, so that only meaningful values remain. +u_2 = u_1.where(cf.lt(-9000), cf.masked) +v_2 = v_1.where(cf.lt(-9000), cf.masked) + +# %% +# 7. Calculate the divergence using the 'div_xy' function operating on the +# vector eastward and northward components as the first and second argument +# respectively. We need to calculate this for the latitude-longitude plane +# of the Earth, defined in spherical polar coordinates, so we must specify +# the Earth's radius for the appropriate calculation: +div = cf.div_xy(u_2, v_2, radius="earth") + +# %% +# 8. First we configure the overall plot by +# making the map higher resolution, to show the coastlines of the UK and +# Ireland in greater detail, and changing the colourmap to better reflect +# the data which can be positive or negative, i.e. has 0 as the 'middle' +# value of significance, so should use a diverging colour map. +cfp.mapset(resolution="10m") +cfp.cscale("ncl_default", ncols=21) + +# %% +# 9. Now generate the final plot. Plot the current vectors, noting we had +# to play around with the 'stride' and 'scale' parameter values to adjust +# the vector spacing and size so that the vector field is best represented +# and visible without over-cluttering the plot. Finally we plot the +# divergence as a contour plot without any lines showing. This compound +# plot is saved on one canvas using 'gopen' and 'gclose' to wrap the two +# plotting calls: +cfp.gopen() +cfp.vect(u=u_2, v=v_2, stride=6, scale=3, key_length=1) +cfp.con( + div, + lines=False, + title=( + f"Depth-averaged Irish Sea currents at {chosen_time} with " + "their divergence" + ), +) +cfp.gclose() diff --git a/recipes-docs/source/recipes-source/plot_22_recipe.py b/recipes-docs/source/recipes-source/plot_22_recipe.py new file mode 100644 index 0000000000..20a77eeb86 --- /dev/null +++ b/recipes-docs/source/recipes-source/plot_22_recipe.py @@ -0,0 +1,121 @@ +""" +Plotting a wind rose as a scatter plot +====================================== + +Given a file containing northerly and easterly wind components, we can +calculate the magnitude and bearing of the resultant wind at each point +in the region and plot them using a scatter plot on a polar grid to +create a wind rose representing wind vectors in the given area. +""" + +# %% +# 1. Import cf-python, Dask.array, NumPy, and Matplotlib: +import dask.array as da +import matplotlib.pyplot as plt +import numpy as np + +import cf + +# %% +# 2. Read the field constructs and load the wind speed component fields: + +f = cf.read("~/recipes/data1.nc") +print(f) + +U = f[2].squeeze() # Easterly wind speed component +V = f[3].squeeze() # Northerly wind speed component + + +# %% +# 3. Set a bounding region for the data and discard readings outside of it: + +tl = (41, 72) # (long, lat) of top left of bounding box. +br = (65, 46) # (long, lat) of bottom right of bounding box. + +U_region = U.subspace(X=cf.wi(tl[0], br[0]), Y=cf.wi(br[1], tl[1])) +V_region = V.subspace(X=cf.wi(tl[0], br[0]), Y=cf.wi(br[1], tl[1])) + +# %% +# 4. Select measurements for a specific pressure using the subspace method, +# then use squeeze to remove the size 1 axis: + +U_sub = U_region.subspace(pressure=500.0) +V_sub = V_region.subspace(pressure=500.0) + +U_sub.squeeze(inplace=True) +V_sub.squeeze(inplace=True) + +# %% +# 5. Calculate the magnitude of each resultant vector using Dask's hypot +# function: + +magnitudes = da.hypot(U_sub.data, V_sub.data) + +# %% +# 6. Calculate the angle of the resultant vector (relative to an Easterly ray) +# using Dask's arctan2 function, then convert to a clockwise bearing: + +azimuths = da.arctan2(V_sub.data, U_sub.data) + +bearings = ((np.pi / 2) - azimuths) % (np.pi * 2) + +# %% +# 7. Flatten the two dimensions of each array for plotting with Matplotlib: + +bearings_flattened = da.ravel(bearings) + +magnitudes_flattened = da.ravel(magnitudes) + +# %% +# 8. Draw the scatter plot using Matplotlib: + +plt.figure(figsize=(5, 6)) +# sphinx_gallery_start_ignore +fig = plt.gcf() +# sphinx_gallery_end_ignore +ax = plt.subplot(polar=True) +ax.set_theta_zero_location("N") # Place 0 degrees at the top. +ax.set_theta_direction(-1) # Arrange bearings clockwise around the plot. +# sphinx_gallery_start_ignore +plt.close() +# sphinx_gallery_end_ignore + +# %% +# 9. Draw a scatter plot on the polar plot, using the wind direction bearing as +# the angle and the magnitude of the resultant wind speed as the distance +# from the pole. +ax.scatter(bearings_flattened.compute(), magnitudes_flattened.compute(), s=1.2) +# sphinx_gallery_start_ignore +plt.close() +# sphinx_gallery_end_ignore + +# %% +# 10. Label the axes and add a title. + +# sphinx_gallery_start_ignore +plt.figure(fig) +# sphinx_gallery_end_ignore +plt.title( + f"Wind Rose Scatter Plot\nLat: {br[1]}°-{tl[1]}°, Long: {tl[0]}°-{br[0]}°" +) + +ax.set_xlabel("Bearing [°]") + +ax.set_ylabel("Speed [m/s]", rotation=45, labelpad=30, size=8) + +ax.yaxis.set_label_coords(0.45, 0.45) + +ax.yaxis.set_tick_params(which="both", labelrotation=45, labelsize=8) + +ax.set_rlabel_position(45) +# sphinx_gallery_start_ignore +plt.close() +# sphinx_gallery_end_ignore + +# %% +# 11. Display the plot. + +# sphinx_gallery_start_ignore +plt.figure(fig) +# sphinx_gallery_end_ignore +plt.show() diff --git a/recipes-docs/source/recipes-source/plot_23_recipe.py b/recipes-docs/source/recipes-source/plot_23_recipe.py new file mode 100644 index 0000000000..4ae11c3863 --- /dev/null +++ b/recipes-docs/source/recipes-source/plot_23_recipe.py @@ -0,0 +1,174 @@ +""" +Combining cf and Matplotlib plots in one figure +=============================================== + +Presently, cf-plot has very few exposed interfaces to its Matplotlib and +Cartopy backend. This makes it difficult to combine plots from the three +in one figure, but not impossible. + +A combined cf and Matplotlib plot can be achieved by amending the figure +stored at ``cfp.plotvars.master_plot``, and then redrawing it with the +new subplots. +""" + +# %% +# 1. Import cf-python, cf-plot, Matplotlib, NumPy, and Dask.array: + +# sphinx_gallery_start_ignore +# sphinx_gallery_thumbnail_number = 2 +# sphinx_gallery_end_ignore + +import matplotlib.pyplot as plt +import numpy as np +import dask.array as da + +import cfplot as cfp +import cf + + +# %% +# 2. Read example data field constructs, and set region for our plots: + +f = cf.read(f"~/recipes/data1.nc") + +u = f.select_by_identity("eastward_wind")[0] +v = f.select_by_identity("northward_wind")[0] +t = f.select_by_identity("air_temperature")[0] + +# Subspace to get values for a specified pressure, here 500 mbar +u = u.subspace(pressure=500) +v = v.subspace(pressure=500) +t = t.subspace(pressure=500) + +lonmin, lonmax, latmin, latmax = 10, 120, -30, 30 + +# %% [markdown] +# +# Outlining the figure with cf-plot +# --------------------------------- +# + +# %% +# 1. Set desired dimensions for our final figure: + +rows, cols = 2, 2 + +# %% +# 2. Create a figure of set dimensions with ``cfp.gopen()``, then set the +# position of the cf plot: + + +cfp.gopen(rows, cols) + +pos = 2 # Second position in the figure + +cfp.gpos(pos) +# sphinx_gallery_start_ignore +plt.close() +# sphinx_gallery_end_ignore + +# %% +# 3. Create a simple vector plot: + +cfp.mapset(lonmin=lonmin, lonmax=lonmax, latmin=latmin, latmax=latmax) +cfp.vect(u=u, v=v, key_length=10, scale=120, stride=4) + +# %% [markdown] +# +# Creating our Matplotlib plots +# ----------------------------- +# + +# %% +# 1. Access the newly-created figure: + +fig = cfp.plotvars.master_plot + +# %% +# 2. Reduce fields down to our test data for a wind rose scatter plot: + +# Limit to specific geographic region +u_region = u.subspace(X=cf.wi(lonmin, lonmax), Y=cf.wi(latmin, latmax)) +v_region = v.subspace(X=cf.wi(lonmin, lonmax), Y=cf.wi(latmin, latmax)) +t_region = t.subspace(X=cf.wi(lonmin, lonmax), Y=cf.wi(latmin, latmax)) + +# Remove size 1 axes +u_squeeze = u_region.squeeze() +v_squeeze = v_region.squeeze() +t_squeeze = t_region.squeeze() + +# Flatten to one dimension for plot +u_f = da.ravel(u_squeeze.data) +v_f = da.ravel(v_squeeze.data) +t_f = da.ravel(t_squeeze.data) + +# %% +# 3. Perform calculations to create appropriate plot data: + +mag_f = da.hypot(u_f, v_f) # Wind speed magnitude + +azimuths_f = da.arctan2(v_f, u_f) +rad_f = ((np.pi / 2) - azimuths_f) % (np.pi * 2) # Wind speed bearing + +# Normalise temperature data into a range appropriate for setting point sizes (1-10pt). +temp_scaled = 1 + (t_f - t_f.min()) / (t_f.max() - t_f.min()) * (10 - 1) + +# %% +# 4. Add Matplotlib subplot to our existing cf figure: + +pos = 1 # First position in the figure + +ax = fig.add_subplot(rows, cols, pos, polar=True) +ax.set_theta_zero_location("N") +ax.set_theta_direction(-1) + +ax.scatter( + rad_f.compute(), + mag_f.compute(), + s=temp_scaled.compute(), + c=temp_scaled.compute(), + alpha=0.5, +) + +ax.set_xlabel("Bearing [°]") +ax.set_ylabel("Speed [m/s]", rotation=45, labelpad=30, size=8) +ax.yaxis.set_label_coords(0.45, 0.45) +ax.yaxis.set_tick_params(which="both", labelrotation=45, labelsize=8) +ax.set_rlabel_position(45) + +# %% +# 5. Create and add a third plot, for example: + +x = np.linspace(0, 10, 100) +y = np.sin(x) + +pos = 3 # Third position in the figure + +ax1 = fig.add_subplot(rows, cols, pos) + +ax1.plot(x, y, label="sin(x)") +ax1.legend() + +# %% [markdown] +# +# Drawing the new figure +# ---------------------- +# + +# %% +# 1. Draw final figure: + +fig = plt.figure(fig) +fig.tight_layout() +fig.show() + +# %% [markdown] +# +# Summary +# ------- +# +# In summary, to use other plotting libraries with cf-plot, you must first +# create your figure with cf-plot with placeholders for your other plots, +# then add subplots by accessing the ``cfp.plotvars.master_plot`` object, +# and finally redraw the figure containing the new plots. + diff --git a/requirements.txt b/requirements.txt index 13ddf3cf40..66cd2d0e36 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,7 +1,7 @@ netCDF4==1.7.2 cftime>=1.6.4 numpy>=2.0.0 -cfdm>=1.12.3.1, <1.12.4.0 +cfdm>=1.13.0.0, <1.13.1.0 psutil>=0.6.0 cfunits>=3.3.7 dask>=2025.5.1 diff --git a/setup.py b/setup.py index 2ed416976d..4620e64a4a 100755 --- a/setup.py +++ b/setup.py @@ -198,8 +198,6 @@ def compile(): * create subspaces of field constructs, -* write field constructs to netCDF datasets on disk, - * incorporate, and create, metadata stored in external files, * read, write, and create data that have been compressed by convention