GeoStat-Framework / GSTools

GSTools - A geostatistical toolbox: random fields, variogram estimation, covariance models, kriging and much more
https://geostat-framework.org
GNU Lesser General Public License v3.0
544 stars 71 forks source link

Add 'temporal', 'spatial_dim' and 'geo_scale' attributes to Covmodel #308

Closed MuellerSeb closed 1 year ago

MuellerSeb commented 1 year ago

Closes #223, #282

This PR adds a boolean temporal and spatial_dim attribute to the CovModel class as well as geo_scale to CovModel and vario_estimate.

Setting temporal=True indicates that the model should be seen as a metric spatio-temporal model. This was possible before by increasing the model dimension by one but this feature basically enables using spatio-temporal model with latlon coordinates. spatial_dim was added so the spatial dimension could be set explicitly and not implicitly by increasing only dim by 1.

Also a geo_scale attribute was added to be able to get meaningful length-scales and bins for latlon models. Before we always used the rescale attribute to do that, but that was rather a hack. Now both can be used to set the sphere radius and to rescale the length-scale. In addition to that I added KM_SCALE, DEGREE_SCALE and RADIAN_SCALE constant that can be used as radius to get a length-scale in km, degree or radian as some people expect the length-scale to be in degree when working latlon coordinates.

Another dimension indicator attribute was added to the CovModel class:

Works like this:

import numpy as np
import gstools as gs

model = gs.Matern(
    temporal=True,
    latlon=True,
    len_scale=[1000, 100],
    geo_scale=gs.KM_SCALE,
)

lat = lon = np.linspace(-80, 81, 50)
time = np.linspace(0, 777, 50)
srf = gs.SRF(model, seed=1234)
field = srf.structured((lat, lon, time))
srf.plot()

Figure_1 Figure_2

MuellerSeb commented 1 year ago

@LSchueler would it make sense to you to move the time axis position in the pos tuple to the front in case of latlon coordinates? This would then follow CF-Conventions: T-Lat-Lon (and not Lat-Lon-T) which is basically T-Y-X (fully inversed axis order in contrast to normal (x,y,z,...,t)).

For non latlon fields, I would keep (x,y,z,...,t).

MuellerSeb commented 1 year ago

@LSchueler Maybe radius is a bit misleading, since it is a rather technical term with respect to the internal representation of the yadrenko-variogram. May we could call this attribute geo_scale witch is a bit more self-descriptive. What do you think?

MuellerSeb commented 1 year ago

@LSchueler I like the geo_scale attribute. Next question, following #282 is, if we want to add geo_scale also to the vario_estimate routine so the bin_centers have a meaningful unit?

MuellerSeb commented 1 year ago

vario_estimate was update with geo_scale and **std_bins arguments. Now everything feels really straight forward.

MuellerSeb commented 1 year ago

@LSchueler we may disable rotation between spatial dimensions and time if time=True. Otherwise I think the metric spatio-temporal model is violated.

LSchueler commented 1 year ago

Sorry for responding so late. I was sick.

This looks interesting! Just a few very brief comments:

MuellerSeb commented 1 year ago

Thanks for your thoughts @LSchueler. Maybe time_dim could be a compromise? Since we are basically adding a "time-dimension"?

MuellerSeb commented 1 year ago

Or only temporal?

LSchueler commented 1 year ago

You find spatio-temporal too long? I guess I could live with temporal ;-)

MuellerSeb commented 1 year ago

@LSchueler all fixes implemented. I would also vote against the CF conventions and maybe have that with #140

MuellerSeb commented 1 year ago

@LSchueler Only question I still have is, if dim should really be automatically increased by 1 if temporal=True.

For example if you want a model with 3 spatial dimensions and a time dimension, which option looks better to you:

model = gs.Gaussian(dim=3, temporal=True)

or

model = gs.Gaussian(dim=4, temporal=True)

In both cases model.dim would be 4 and model.spatial_dim would be 3.

MuellerSeb commented 1 year ago

I added spatial_dim to solve this. dim now always really refers to the model dimension. spatial_dim could be used instead, when defining a spatio-temporal model:

model = gs.Gaussian(temporal=True, spatial_dim=3)
model.dim == 4
model.spatial_dim == 3

This is equal to:

model = gs.Gaussian(temporal=True, dim=4)