Lat-Lon support (geographical coordinates)

[MuellerSeb]

To really earn the name "Geo"Stat-tools, it would be nice to support lat-lon coordinates in the kriging and field-generation routines.

The RandMeth routines are not able to work with latlon coordinates since it's an approximation of the fourier-transformation in euclidean space which differs from the FT on the manifold of a sphere.

A trick to prevent this problem could be to use 3D fields and evaluate the field on a sphere. For rather small correlation lengths this shouldn't be a problem but with increasing length-scales the error between the spatial distance and the great-circle distance gets bigger.

Within the kriging routines it should be quite simple to implement this, since only the pairwise distances of the points are calculated. At this point we could install a switch to use the great-circle distance (or sth more fancy from geopy)

Also the variogram estimation should be easy adoptable, since only the distance calculation has to be adjusted.

Since no anisotropy and consequently no rotation is allowed for geographical coordinates, it would be a good idea to implement this as a switch in the CovModel. Then the following should be set:

• dim=2 (here we could get a problem, since the above mentioned fix in RandMeth needs dim=3)
• anis=[1]
• angles=[0]

Maybe it could be a good idea to freeze them.

We could add a keyword-argument latlon or geo_coord to the CovModel.__init__, that is False by default.

What do you think? @LSchueler, @fhesze, @AlrauneZ, @banesullivan

[LSchueler]

Before putting too much work into this, we should actually check the differences between a naive approach like the one pasted below and the correct computation on a sphere in dependence of the area on which the fields are computed.

Does anyone have an idea what the largest areas are, on which random fields are generated? At least in oceanography, nobody cares about spherical coordinates for areas smaller than 50km x 50km.

Are you suggesting to neither support anisotropic nor rotated anisotropic fields? - That would reduce the potential applications dramatically.
For global fields (which are the size of a significant chunk of the globe), I get that anisotropy cannot be defined. But as an approximation for much smaller areas, which are still defined on lat-lon-coords, anisotropic fields should be supported.

This is just for visualisation, I would not suggest to generate a field for a quarter of the globe.

import numpy as np
import plotly.graph_objects as go
from plotly.subplots import make_subplots

from gstools import SRF, Gaussian

def spherical_linspace(start, stop, steps):
    """Takes degrees, spits out rads."""
    return np.linspace(start * np.pi/180.0, stop * np.pi/180.0, steps)

def spherical_coords(theta, phi, R):
    x = R * np.cos(theta) * np.cos(phi)
    y = R * np.cos(theta) * np.sin(phi)
    z = R * np.sin(theta) * np.ones_like(phi)
    return x, y, z


theta_deg = np.array((-20.0, 50.0))
phi_deg = np.array((0.0, 180.0))
N_theta = 40
N_phi = 80

theta = spherical_linspace(theta_deg[0], theta_deg[1], N_theta)
phi = spherical_linspace(phi_deg[0], phi_deg[1], N_phi)

theta = np.reshape(theta, (len(theta), 1))
phi = np.reshape(phi, (1, len(phi)))
R = 1.0

x, y, z = spherical_coords(theta, phi, R)

model = Gaussian(dim=3, var=1, len_scale=0.01)
srf = SRF(model, seed=474636)

# at the moment, it is a bit ugly to use an unstructured grid
# and then reshape it to a structured one
field = srf((x, y, z), mesh_type='unstructured')
field_reshape = field.reshape((N_theta, N_phi))

mesh = srf.to_pyvista()
mesh.plot()

fig = make_subplots(rows=1, cols=1,
                    specs=[[{'is_3d': True}]],
                    )

fig.add_trace(go.Surface(x=x, y=y, z=z, surfacecolor=field_reshape), 1, 1)
fig.update_layout(title_text="Spherical Coords")
fig.show()

srf.plot()

[MuellerSeb]

Approximation of small areas

We could provide a procedure to convert lat-lon coordinates with an euclidean geodesic approximation like described here: link, when the error introduced by earths curvature is relatively small.

We could also use the tangent-plane at the center of the given latlon coordinates that is directing to the north-pole and use orthogonal projection onto this plane (like done in some alternative world maps):

Milestones for the general procedure:

• add latlon switch to CovModel (setting dim=2, anis=1, angels=0)
• suppress anisotropy in the latlon case
• apply a hack to the Randmeth class:
  

  GSTools/gstools/field/generator.py

  Line 449 in 32aa702

  def _reshape_pos(x, y=None, z=None, dtype=np.double):

  • this procedure could convert latlon data to 3D point sets on a sphere (to apply the workflow from @LSchueler)
  • there we can also apply a correction of the length-scale similar to this procedure: https://github.com/bsmurphy/PyKrige/blob/a31dbbb9ee22d89da2378f28dc11c0d882ae564d/pykrige/core.py#L93
  • To correctly sample in the randmeth class, we have to modify
    

    GSTools/gstools/field/generator.py

    Line 203 in 32aa702

    def reset_seed(self, seed=np.nan):

    • switch of latlon for sampling and set dim=3
    • since we have a deepcopy of the model in generator, we don't harm the model in the SRF class
    • after generating the rad-samples, switch back to latlon
• since incompressible random vector fields can not have a unique mean flow direction on a sphere, we should prohibit them for the latlon case
• to provide kriging on latlon coordinates we only have to provide the Harversine equivalent to scipy's cdist (should be quite easy as a small cython module)

[MuellerSeb]

To resolve the "latlon" vs "lonlat" discussion, look here:
https://en.wikipedia.org/wiki/ISO_6709

We should stick to latlon

[MuellerSeb]

This could be interesting to generate isotropic random fields on the sphere:
https://www.degruyter.com/view/journals/mcma/24/1/article-p1.xml

[MuellerSeb]

This one comes even closer to our workflow:

"Spectral Simulation of Isotropic Gaussian Random Fields on a Sphere"

https://link.springer.com/article/10.1007%2Fs11004-019-09799-4

[MuellerSeb]

As pointed out by @mjziebarth here: GeoStat-Framework/PyKrige#149, the described Yadrenko models produce the same fields like the "cut-out" fields described by @LSchueler. Therefore we could provide these fields as the standard geographic fields. We just need to reference the Yadrenko models as the standard way of interpreting geo-coordinates.

We would need to provide some convenient functionality to take lat-lon input and generate 3D fields.

Regarding the underlying Mesh class from #65, we could store the 3D field and save the lat-lon information as additional point data. What do you think @LSchueler?

[MuellerSeb]

We could also produce fields with additional elevation above the "earth-radius" (e.g. mean sea level) elev. Then we can provide anisotropy in vertical direction. There we need to provide the rescale factor mentioned in #90. Then the isotropified positions in space are:

x = (1 + elev / (rescale * anis)) * cos(lat) * cos(lon)
y = (1 + elev / (rescale * anis)) * cos(lat) * sin(lon)
z = (1 + elev / (rescale * anis)) * sin(lat)

I would suggest to add a geo_dim parameter to the CovModel class, that is None by default and can be set to 2 (planar lat-lon fields) or 3 (lat-lon fields with elevation). If geo_dim is 2 or 3 the spatial dim is set to 3 (since the sphere is living in the 3-dim space and so the RandMeth class can handle it).

The CovModel class should provide routines to transform a given pos array to the isotropified positions. There we can also take care of the special geo_dim variants.

[MuellerSeb]

Closed by #113