Merge pull request #113 from GeoStat-Framework/latlon_support

Add support for geographical coordinates
GeoStat-Framework · Jan 13, 2021 · c1b56da · c1b56da
2 parents 26c727d + 606c4ff
commit c1b56da
Show file tree

Hide file tree

Showing 39 changed files with 2,000 additions and 231 deletions.
diff --git a/README.md b/README.md
@@ -81,6 +81,7 @@ The documentation also includes some [tutorials][tut_link], showing the most imp
 - [Kriging][tut5_link]
 - [Conditioned random field generation][tut6_link]
 - [Field transformations][tut7_link]
+- [Geographic Coordinates][tut8_link]
 - [Miscellaneous examples][tut0_link]
 
 The associated python scripts are provided in the `examples` folder.
@@ -112,23 +113,47 @@ srf.plot()
 <img src="https://raw.githubusercontent.com/GeoStat-Framework/GSTools/master/docs/source/pics/gau_field.png" alt="Random field" width="600px"/>
 </p>
 
+GSTools also provides support for [geographic coordinates](https://en.wikipedia.org/wiki/Geographic_coordinate_system).
+This works perfectly well with [cartopy](https://scitools.org.uk/cartopy/docs/latest/index.html).
+
+```python
+import matplotlib.pyplot as plt
+import cartopy.crs as ccrs
+import gstools as gs
+# define a structured field by latitude and longitude
+lat = lon = range(-80, 81)
+model = gs.Gaussian(latlon=True, len_scale=777, rescale=gs.EARTH_RADIUS)
+srf = gs.SRF(model, seed=12345)
+field = srf.structured((lat, lon))
+# Orthographic plotting with cartopy
+ax = plt.subplot(projection=ccrs.Orthographic(-45, 45))
+cont = ax.contourf(lon, lat, field, transform=ccrs.PlateCarree())
+ax.coastlines()
+ax.set_global()
+plt.colorbar(cont)
+```
+
+<p align="center">
+<img src="https://github.com/GeoStat-Framework/GeoStat-Framework.github.io/raw/master/img/GS_globe.png" alt="lat-lon random field" width="600px"/>
+</p>
+
 A similar example but for a three dimensional field is exported to a [VTK](https://vtk.org/) file, which can be visualized with [ParaView](https://www.paraview.org/) or [PyVista](https://docs.pyvista.org) in Python:
 
 ```python
 import gstools as gs
 # structured field with a size 100x100x100 and a grid-size of 1x1x1
 x = y = z = range(100)
-model = gs.Gaussian(dim=3, var=0.6, len_scale=20)
+model = gs.Gaussian(dim=3, len_scale=[16, 8, 4], angles=(0.8, 0.4, 0.2))
 srf = gs.SRF(model)
 srf((x, y, z), mesh_type='structured')
 srf.vtk_export('3d_field') # Save to a VTK file for ParaView
 
 mesh = srf.to_pyvista() # Create a PyVista mesh for plotting in Python
-mesh.threshold_percent(0.5).plot()
+mesh.contour(isosurfaces=8).plot()
 ```
 
 <p align="center">
-<img src="https://raw.githubusercontent.com/GeoStat-Framework/GSTools/master/docs/source/pics/3d_gau_field.png" alt="3d Random field" width="600px"/>
+<img src="https://github.com/GeoStat-Framework/GeoStat-Framework.github.io/raw/master/img/GS_pyvista.png" alt="3d Random field" width="600px"/>
 </p>
 
 
@@ -335,6 +360,7 @@ You can contact us via <[email protected]>.
 [tut5_link]: https://geostat-framework.readthedocs.io/projects/gstools/en/stable/examples/05_kriging/index.html
 [tut6_link]: https://geostat-framework.readthedocs.io/projects/gstools/en/stable/examples/06_conditioned_fields/index.html
 [tut7_link]: https://geostat-framework.readthedocs.io/projects/gstools/en/stable/examples/07_transformations/index.html
+[tut8_link]: https://geostat-framework.readthedocs.io/projects/gstools/en/stable/examples/08_geo_coordinates/index.html
 [tut0_link]: https://geostat-framework.readthedocs.io/projects/gstools/en/stable/examples/00_misc/index.html
 [cor_link]: https://en.wikipedia.org/wiki/Autocovariance#Normalization
 [vtk_link]: https://www.vtk.org/
diff --git a/docs/source/conf.py b/docs/source/conf.py
@@ -265,6 +265,7 @@ def setup(app):
 from sphinx_gallery.sorting import FileNameSortKey
 
 sphinx_gallery_conf = {
+    "remove_config_comments": True,
     # only show "print" output as output
     "capture_repr": (),
     # path to your examples scripts
@@ -277,6 +278,7 @@ def setup(app):
         "../../examples/05_kriging/",
         "../../examples/06_conditioned_fields/",
         "../../examples/07_transformations/",
+        "../../examples/08_geo_coordinates/",
     ],
     # path where to save gallery generated examples
     "gallery_dirs": [
@@ -288,6 +290,7 @@ def setup(app):
         "examples/05_kriging/",
         "examples/06_conditioned_fields/",
         "examples/07_transformations/",
+        "examples/08_geo_coordinates/",
     ],
     # Pattern to search for example files
     "filename_pattern": r"\.py",

diff --git a/docs/source/index.rst b/docs/source/index.rst
@@ -97,6 +97,7 @@ showing the most important use cases of GSTools, which are
 - `Kriging <examples/05_kriging/index.html>`__
 - `Conditioned random field generation <examples/06_conditioned_fields/index.html>`__
 - `Field transformations <examples/07_transformations/index.html>`__
+- `Geographic Coordinates <examples/08_geo_coordinates/index.html>`__
 - `Miscellaneous examples <examples/00_misc/index.html>`__
 
 Some more examples are provided in the examples folder.
@@ -133,6 +134,30 @@ with a :any:`Gaussian` covariance model.
    :width: 400px
    :align: center
 
+GSTools also provides support for `geographic coordinates <https://en.wikipedia.org/wiki/Geographic_coordinate_system>`_.
+This works perfectly well with `cartopy <https://scitools.org.uk/cartopy/docs/latest/index.html>`_.
+
+.. code-block:: python
+
+    import matplotlib.pyplot as plt
+    import cartopy.crs as ccrs
+    import gstools as gs
+    # define a structured field by latitude and longitude
+    lat = lon = range(-80, 81)
+    model = gs.Gaussian(latlon=True, len_scale=777, rescale=gs.EARTH_RADIUS)
+    srf = gs.SRF(model, seed=12345)
+    field = srf.structured((lat, lon))
+    # Orthographic plotting with cartopy
+    ax = plt.subplot(projection=ccrs.Orthographic(-45, 45))
+    cont = ax.contourf(lon, lat, field, transform=ccrs.PlateCarree())
+    ax.coastlines()
+    ax.set_global()
+    plt.colorbar(cont)
+
+.. image:: https://github.com/GeoStat-Framework/GeoStat-Framework.github.io/raw/master/img/GS_globe.png
+   :width: 400px
+   :align: center
+
 A similar example but for a three dimensional field is exported to a
 `VTK <https://vtk.org/>`__ file, which can be visualized with
 `ParaView <https://www.paraview.org/>`_ or
@@ -143,15 +168,15 @@ A similar example but for a three dimensional field is exported to a
     import gstools as gs
     # structured field with a size 100x100x100 and a grid-size of 1x1x1
     x = y = z = range(100)
-    model = gs.Gaussian(dim=3, var=0.6, len_scale=20)
+    model = gs.Gaussian(dim=3, len_scale=[16, 8, 4], angles=(0.8, 0.4, 0.2))
     srf = gs.SRF(model)
     srf((x, y, z), mesh_type='structured')
     srf.vtk_export('3d_field') # Save to a VTK file for ParaView
 
     mesh = srf.to_pyvista() # Create a PyVista mesh for plotting in Python
-    mesh.threshold_percent(0.5).plot()
+    mesh.contour(isosurfaces=8).plot()
 
-.. image:: https://raw.githubusercontent.com/GeoStat-Framework/GSTools/master/docs/source/pics/3d_gau_field.png
+.. image:: https://github.com/GeoStat-Framework/GeoStat-Framework.github.io/raw/master/img/GS_pyvista.png
    :width: 400px
    :align: center
 

diff --git a/docs/source/tutorials.rst b/docs/source/tutorials.rst
@@ -19,4 +19,5 @@ explore its whole beauty and power.
    examples/05_kriging/index
    examples/06_conditioned_fields/index
    examples/07_transformations/index
+   examples/08_geo_coordinates/index
    examples/00_misc/index
diff --git a/examples/00_misc/03_precipitation.py b/examples/00_misc/03_precipitation.py
@@ -20,12 +20,12 @@
 # fix the seed for reproducibility
 seed = 20170521
 # half daily timesteps over three months
-t = np.arange(0., 90., 0.5)
+t = np.arange(0.0, 90.0, 0.5)
 # spatial axis of 50km with a resolution of 1km
 x = np.arange(0, 50, 1.0)
 
 # an exponential variogram with a corr. lengths of 2d and 5km
-model = gs.Exponential(dim=2, var=1.0, len_scale=2., anis=2.5)
+model = gs.Exponential(dim=2, var=1.0, len_scale=2.0, anis=2.5)
 # create a spatial random field instance
 srf = gs.SRF(model)
 
@@ -68,46 +68,46 @@
 
 fig, axs = plt.subplots(2, 2, sharex=True, sharey=True)
 
-axs[0,0].set_title('Gaussian')
-axs[0,0].plot(t, P_gau[:, 20])
-axs[0,0].set_ylabel(r'$P$ / mm')
+axs[0, 0].set_title("Gaussian")
+axs[0, 0].plot(t, P_gau[:, 20])
+axs[0, 0].set_ylabel(r"$P$ / mm")
 
-axs[0,1].set_title('Cut Gaussian')
-axs[0,1].plot(t, P_cut[:, 20])
+axs[0, 1].set_title("Cut Gaussian")
+axs[0, 1].plot(t, P_cut[:, 20])
 
-axs[1,0].set_title('Cut Gaussian Anamorphosis')
-axs[1,0].plot(t, P_ana[:, 20])
-axs[1,0].set_xlabel(r'$t$ / d')
-axs[1,0].set_ylabel(r'$P$ / mm')
+axs[1, 0].set_title("Cut Gaussian Anamorphosis")
+axs[1, 0].plot(t, P_ana[:, 20])
+axs[1, 0].set_xlabel(r"$t$ / d")
+axs[1, 0].set_ylabel(r"$P$ / mm")
 
-axs[1,1].set_title('Different Cross Section')
-axs[1,1].plot(t, P_ana[:, 10])
-axs[1,1].set_xlabel(r'$t$ / d')
+axs[1, 1].set_title("Different Cross Section")
+axs[1, 1].plot(t, P_ana[:, 10])
+axs[1, 1].set_xlabel(r"$t$ / d")
 
 plt.tight_layout()
 
 fig, axs = plt.subplots(2, 2, sharex=True, sharey=True)
 
-axs[0,0].set_title('Gaussian')
-cont = axs[0,0].contourf(t, x, P_gau.T, cmap='PuBu')
-cbar = fig.colorbar(cont, ax=axs[0,0])
-cbar.ax.set_ylabel(r'$P$ / mm')
-axs[0,0].set_ylabel(r'$x$ / km')
-
-axs[0,1].set_title('Cut Gaussian')
-cont = axs[0,1].contourf(t, x, P_cut.T, cmap='PuBu')
-cbar = fig.colorbar(cont, ax=axs[0,1])
-cbar.ax.set_ylabel(r'$P$ / mm')
-axs[0,1].set_xlabel(r'$t$ / d')
-
-axs[1,0].set_title('Cut Gaussian Anamorphosis')
-cont = axs[1,0].contourf(t, x, P_ana.T, cmap='PuBu')
-cbar = fig.colorbar(cont, ax=axs[1,0])
-cbar.ax.set_ylabel(r'$P$ / mm')
-axs[1,0].set_xlabel(r'$t$ / d')
-axs[1,0].set_ylabel(r'$x$ / km')
-
-fig.delaxes(axs[1,1])
+axs[0, 0].set_title("Gaussian")
+cont = axs[0, 0].contourf(t, x, P_gau.T, cmap="PuBu")
+cbar = fig.colorbar(cont, ax=axs[0, 0])
+cbar.ax.set_ylabel(r"$P$ / mm")
+axs[0, 0].set_ylabel(r"$x$ / km")
+
+axs[0, 1].set_title("Cut Gaussian")
+cont = axs[0, 1].contourf(t, x, P_cut.T, cmap="PuBu")
+cbar = fig.colorbar(cont, ax=axs[0, 1])
+cbar.ax.set_ylabel(r"$P$ / mm")
+axs[0, 1].set_xlabel(r"$t$ / d")
+
+axs[1, 0].set_title("Cut Gaussian Anamorphosis")
+cont = axs[1, 0].contourf(t, x, P_ana.T, cmap="PuBu")
+cbar = fig.colorbar(cont, ax=axs[1, 0])
+cbar.ax.set_ylabel(r"$P$ / mm")
+axs[1, 0].set_xlabel(r"$t$ / d")
+axs[1, 0].set_ylabel(r"$x$ / km")
+
+fig.delaxes(axs[1, 1])
 plt.tight_layout()
 
 ###############################################################################
@@ -123,14 +123,14 @@
 # fix the seed for reproducibility
 seed = 20170521
 # half daily timesteps over three months
-t = np.arange(0., 90., 0.5)
+t = np.arange(0.0, 90.0, 0.5)
 # 1st spatial axis of 50km with a resolution of 1km
 x = np.arange(0, 50, 1.0)
 # 2nd spatial axis of 40km with a resolution of 1km
 y = np.arange(0, 40, 1.0)
 
 # an exponential variogram with a corr. lengths of 2d, 5km, and 5km
-model = gs.Exponential(dim=3, var=1.0, len_scale=2., anis=(2.5, 2.5))
+model = gs.Exponential(dim=3, var=1.0, len_scale=2.0, anis=(2.5, 2.5))
 # create a spatial random field instance
 srf = gs.SRF(model)
 

diff --git a/examples/00_misc/README.rst b/examples/00_misc/README.rst
@@ -5,9 +5,5 @@ More examples which do not really fit into other categories. Some are not more
 than a code snippet, while others are more complex and more than one part of
 GSTools is involved.
 
-.. only:: html
-
-   Gallery
-   -------
-
-   Below is a gallery of examples
+Examples
+--------
diff --git a/examples/01_random_field/04_srf_merge.py b/examples/01_random_field/04_srf_merge.py
@@ -6,6 +6,7 @@
 to merge two unstructured rectangular fields.
 
 """
+# sphinx_gallery_thumbnail_number = 2
 import numpy as np
 import gstools as gs
 

diff --git a/examples/01_random_field/README.rst b/examples/01_random_field/README.rst
@@ -13,9 +13,5 @@ and its discretised modes are evaluated at random frequencies.
 
 GSTools supports arbitrary and non-isotropic covariance models.
 
-.. only:: html
-
-   Gallery
-   -------
-
-   Below is a gallery of examples
+Examples
+--------
diff --git a/examples/01_random_field/mesh_ensemble.vtk b/examples/01_random_field/mesh_ensemble.vtk
diff --git a/examples/02_cov_model/README.rst b/examples/02_cov_model/README.rst
@@ -71,18 +71,17 @@ The following standard covariance models are provided by GSTools
     HyperSpherical
     SuperSpherical
     JBessel
+    TPLSimple
 
 As a special feature, we also provide truncated power law (TPL) covariance models
 
 .. autosummary::
     TPLGaussian
     TPLExponential
     TPLStable
-    TPLSimple
-
-.. only:: html
 
-   Gallery
-   -------
+These models provide a lower and upper length scale truncation
+for superpositioned models.
 
-   Below is a gallery of examples
+Examples
+--------
diff --git a/examples/03_variogram/README.rst b/examples/03_variogram/README.rst
@@ -10,9 +10,5 @@ The same `(semi-)variogram <https://en.wikipedia.org/wiki/Variogram#Semivariogra
 :ref:`tutorial_02_cov` is being used
 by this subpackage.
 
-.. only:: html
-
-   Gallery
-   -------
-
-   Below is a gallery of examples
+Examples
+--------
diff --git a/examples/04_vector_field/README.rst b/examples/04_vector_field/README.rst
@@ -32,9 +32,6 @@ with the projector
 By calculating :math:`\nabla \cdot \mathbf U = 0`, it can be verified, that
 the resulting field is indeed incompressible.
 
-.. only:: html
 
-   Gallery
-   -------
-
-   Below is a gallery of examples
+Examples
+--------
diff --git a/examples/05_kriging/README.rst b/examples/05_kriging/README.rst
@@ -87,9 +87,5 @@ submodule :any:`gstools.krige`.
     ExtDrift
     Detrended
 
-.. only:: html
-
-   Gallery
-   -------
-
-   Below is a gallery of examples
+Examples
+--------
diff --git a/examples/06_conditioned_fields/README.rst b/examples/06_conditioned_fields/README.rst
@@ -30,9 +30,5 @@ or:
 
     srf.set_condition(cond_pos, cond_val, "ordinary")
 
-.. only:: html
-
-   Gallery
-   -------
-
-   Below is a gallery of examples
+Examples
+--------
diff --git a/examples/07_transformations/README.rst b/examples/07_transformations/README.rst
@@ -33,9 +33,5 @@ Simply import the transform submodule and apply a transformation to the srf clas
     ...
     tf.normal_to_lognormal(srf)
 
-.. only:: html
-
-   Gallery
-   -------
-
-   Below is a gallery of examples
+Examples
+--------