From 01819b59167ef2eb649e4214022099f22fa841bf Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Wed, 25 Nov 2020 00:09:24 +0100 Subject: [PATCH 01/39] Variogram: add support for geographical variogram estimation --- gstools/variogram/estimator.pyx | 45 ++++++++++++++++++++++++++++++--- gstools/variogram/variogram.py | 15 +++++++++++ 2 files changed, 56 insertions(+), 4 deletions(-) diff --git a/gstools/variogram/estimator.pyx b/gstools/variogram/estimator.pyx index 4b435216..bd814cd7 100644 --- a/gstools/variogram/estimator.pyx +++ b/gstools/variogram/estimator.pyx @@ -10,7 +10,7 @@ import numpy as np cimport cython from cython.parallel import prange, parallel from libcpp.vector cimport vector -from libc.math cimport fabs, sqrt, isnan, acos, M_PI +from libc.math cimport fabs, sqrt, isnan, acos, pow, sin, cos, atan2, M_PI cimport numpy as np @@ -18,7 +18,7 @@ DTYPE = np.double ctypedef np.double_t DTYPE_t -cdef inline double distance( +cdef inline double dist_euclid( const int dim, const double[:,:] pos, const int i, @@ -31,6 +31,33 @@ cdef inline double distance( return sqrt(dist_squared) +cdef inline double dist_haversine( + const int dim, + const double[:,:] pos, + const int i, + const int j +) nogil: + # pos holds lat-lon in deg + cdef double deg_2_rad = M_PI / 180.0 + cdef double diff_lat = (pos[0, j] - pos[0, i]) * deg_2_rad + cdef double diff_lon = (pos[1, j] - pos[1, i]) * deg_2_rad + cdef double arg = ( + pow(sin(diff_lat/2.0), 2) + + cos(pos[0, i]*deg_2_rad) * + cos(pos[0, j]*deg_2_rad) * + pow(sin(diff_lon/2.0), 2) + ) + return 2.0 * atan2(sqrt(arg), sqrt(1.0-arg)) + + +ctypedef double (*_dist_func)( + const int, + const double[:,:], + const int, + const int +) nogil + + cdef inline bint dir_test( const int dim, const double[:,:] pos, @@ -203,7 +230,7 @@ def directional( for i in prange(i_max, nogil=True): for j in range(j_max): for k in range(j+1, k_max): - dist = distance(dim, pos, j, k) + dist = dist_euclid(dim, pos, j, k) if dist < bin_edges[i] or dist >= bin_edges[i+1]: continue # skip if not in current bin for d in range(d_max): @@ -227,8 +254,18 @@ def unstructured( const double[:,:] f, const double[:] bin_edges, const double[:,:] pos, - str estimator_type='m' + str estimator_type='m', + str distance_type='e' ): + cdef _dist_func distance + + if distance_type == 'e': + distance = dist_euclid + else: + distance = dist_haversine + if dim != 2: + raise ValueError('Haversine: dim = {0} != 2'.format(dim)) + if pos.shape[1] != f.shape[1]: raise ValueError('len(pos) = {0} != len(f) = {1} '. format(pos.shape[1], f.shape[1])) diff --git a/gstools/variogram/variogram.py b/gstools/variogram/variogram.py index c7bdb4f0..c6a2249d 100644 --- a/gstools/variogram/variogram.py +++ b/gstools/variogram/variogram.py @@ -71,6 +71,7 @@ def vario_estimate( sampling_size=None, sampling_seed=None, estimator="matheron", + latlon=False, direction=None, angles=None, angles_tol=np.pi / 8, @@ -141,6 +142,14 @@ def vario_estimate( * "cressie": an estimator more robust to outliers Default: "matheron" + latlon : :class:`bool`, optional + Whether the data is representing 2D fields on earths surface described + by latitude and longitude. When using this, the estimator will + use great-circle distance for variogram estimation. + Note, that only an isotropic variogram can be estimated and an + ValueError will be raised, if direction were specified. + Bin edges need to be given in radians in this case. + Default: False direction : :class:`list` of :class:`numpy.ndarray`, optional directions to evaluate a directional variogram. Anglular tolerance is given by `angles_tol`. @@ -223,6 +232,8 @@ def vario_estimate( pos, __, dim = format_unstruct_pos_shape( pos, field.shape, check_stacked_shape=True ) + if latlon and dim != 2: + raise ValueError("Variogram: given field to be 2D for lat-lon.") # prepare the field pnt_cnt = len(pos[0]) field = field.reshape((-1, pnt_cnt)) @@ -261,6 +272,8 @@ def vario_estimate( dir_no = direction.shape[0] # prepare directional variogram if dir_no > 0: + if latlon: + raise ValueError("Directional variogram not allowed for lat-lon.") norms = np.linalg.norm(direction, axis=1) if np.any(np.isclose(norms, 0)): raise ValueError("Zero length direction {}".format(direction)) @@ -280,12 +293,14 @@ def vario_estimate( cython_estimator = _set_estimator(estimator) # run if dir_no == 0: + distance_type = "s" if latlon else "e" estimates, counts = unstructured( dim, field, bin_edges, pos, estimator_type=cython_estimator, + distance_type=distance_type, ) else: estimates, counts = directional( From a216787b4c87cac2ba74991afa05264ff3e35c3b Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Wed, 25 Nov 2020 00:10:20 +0100 Subject: [PATCH 02/39] Tools: add converter from spatial 3D pos to latlon and vice versa; add earth rad constant --- gstools/tools/__init__.py | 12 ++++++++ gstools/tools/export.py | 4 +-- gstools/tools/geometric.py | 61 ++++++++++++++++++++++++++++++++++++++ 3 files changed, 75 insertions(+), 2 deletions(-) diff --git a/gstools/tools/__init__.py b/gstools/tools/__init__.py index 4d8963d9..9c00c993 100644 --- a/gstools/tools/__init__.py +++ b/gstools/tools/__init__.py @@ -45,6 +45,12 @@ matrix_anisometrize ang2dir +Misc +^^^^ + +.. autosummary:: + EARTH_RADIUS + ---- """ @@ -83,6 +89,11 @@ ang2dir, ) + +EARTH_RADIUS = 6371.0 +"""float: earth radius for WGS84 ellipsoid in km""" + + __all__ = [ "vtk_export", "vtk_export_structured", @@ -110,4 +121,5 @@ "matrix_anisometrize", "rotated_main_axes", "ang2dir", + "EARTH_RADIUS", ] diff --git a/gstools/tools/export.py b/gstools/tools/export.py index 92d8dd03..16a50e18 100644 --- a/gstools/tools/export.py +++ b/gstools/tools/export.py @@ -206,7 +206,7 @@ def to_vtk(pos, fields, mesh_type="unstructured"): # pragma: no cover :class:`pyvista.RectilinearGrid` and unstructured meshes will return an :class:`pyvista.UnstructuredGrid` object. """ - if mesh_type == "structured": + if mesh_type != "unstructured": grid = to_vtk_structured(pos=pos, fields=fields) else: grid = to_vtk_unstructured(pos=pos, fields=fields) @@ -233,6 +233,6 @@ def vtk_export( mesh_type : :class:`str`, optional 'structured' / 'unstructured'. Default: structured """ - if mesh_type == "structured": + if mesh_type != "unstructured": return vtk_export_structured(filename=filename, pos=pos, fields=fields) return vtk_export_unstructured(filename=filename, pos=pos, fields=fields) diff --git a/gstools/tools/geometric.py b/gstools/tools/geometric.py index a02e13b0..d34565c5 100644 --- a/gstools/tools/geometric.py +++ b/gstools/tools/geometric.py @@ -24,6 +24,8 @@ format_struct_pos_shape format_unstruct_pos_shape ang2dir + latlon2pos + pos2latlon """ # pylint: disable=C0103 @@ -47,6 +49,8 @@ "format_struct_pos_shape", "format_unstruct_pos_shape", "ang2dir", + "latlon2pos", + "pos2latlon", ] @@ -580,3 +584,60 @@ def ang2dir(angles, dtype=np.double, dim=None): if dim in [2, 3]: vec[:, [0, 1]] = vec[:, [1, 0]] # to match convention in 2D and 3D return vec + + +def latlon2pos(latlon, radius=1.0, dtype=np.double): + """Convert lat-lon geo coordinates to 3D position tuple. + + Parameters + ---------- + latlon : :class:`list` of :class:`numpy.ndarray` + latitude and longitude given in degrees. + radius : :class:`float`, optional + Earth radius. Default: `1.0` + dtype : data-type, optional + The desired data-type for the array. + If not given, then the type will be determined as the minimum type + required to hold the objects in the sequence. Default: None + + Returns + ------- + :class:`numpy.ndarray` + the 3D position array + """ + latlon = np.array(latlon, dtype=dtype).reshape((2, -1)) + lat, lon = np.deg2rad(latlon) + return np.array( + ( + radius * np.cos(lat) * np.cos(lon), + radius * np.cos(lat) * np.sin(lon), + radius * np.sin(lat) * np.ones_like(lon), + ), + dtype=dtype, + ) + + +def pos2latlon(pos, radius=1.0, dtype=np.double): + """Convert 3D position tuple from sphere to lat-lon geo coordinates. + + Parameters + ---------- + pos : :class:`list` of :class:`numpy.ndarray` + The position tuple containing points on a unit-sphere. + radius : :class:`float`, optional + Earth radius. Default: `1.0` + dtype : data-type, optional + The desired data-type for the array. + If not given, then the type will be determined as the minimum type + required to hold the objects in the sequence. Default: None + + Returns + ------- + :class:`numpy.ndarray` + the 3D position array + """ + pos = np.array(pos, dtype=dtype).reshape((3, -1)) + # prevent numerical errors in arcsin + lat = np.arcsin(np.maximum(np.minimum(pos[2] / radius, 1.0), -1.0)) + lon = np.arctan2(pos[1], pos[0]) + return np.rad2deg((lat, lon), dtype=dtype) From a8170064a8418166abfacbb90a3aa0e39e294669 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Wed, 25 Nov 2020 00:11:09 +0100 Subject: [PATCH 03/39] Init: add EARTH_RADIUS; fix minor sphinx bug --- gstools/__init__.py | 16 ++++++++++++---- 1 file changed, 12 insertions(+), 4 deletions(-) diff --git a/gstools/__init__.py b/gstools/__init__.py index c29faee1..d369555d 100644 --- a/gstools/__init__.py +++ b/gstools/__init__.py @@ -38,7 +38,7 @@ ^^^^^^^^^^^^^^^^^^^^^ Class to construct user defined covariance models -.. currentmodule:: gstools.covmodel.base +.. currentmodule:: gstools.covmodel .. autosummary:: CovModel @@ -49,8 +49,6 @@ Standard Covariance Models ~~~~~~~~~~~~~~~~~~~~~~~~~~ -.. currentmodule:: gstools.covmodel.models - .. autosummary:: Gaussian Exponential @@ -67,7 +65,6 @@ Truncated Power Law Covariance Models ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -.. currentmodule:: gstools.covmodel.tpl_models .. autosummary:: TPLGaussian TPLExponential @@ -103,12 +100,22 @@ .. autosummary:: vario_estimate vario_estimate_axis + +Misc +==== + +.. currentmodule:: gstools.tools + +.. autosummary:: + EARTH_RADIUS + """ # Hooray! from gstools import field, variogram, random, covmodel, tools, krige, transform from gstools.field import SRF from gstools.tools import ( rotated_main_axes, + EARTH_RADIUS, vtk_export, vtk_export_structured, vtk_export_unstructured, @@ -179,6 +186,7 @@ __all__ += [ "SRF", "rotated_main_axes", + "EARTH_RADIUS", "vtk_export", "vtk_export_structured", "vtk_export_unstructured", From 1d0c107e41e7229cb89725c4c3a9d76d1646d998 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Wed, 25 Nov 2020 00:12:21 +0100 Subject: [PATCH 04/39] Krige: minor bugfix for internal drift --- gstools/krige/base.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/gstools/krige/base.py b/gstools/krige/base.py index c468bcb6..4fab9d23 100755 --- a/gstools/krige/base.py +++ b/gstools/krige/base.py @@ -267,8 +267,7 @@ def _get_krige_vecs( res[self.cond_no, :] = 1 # drift function need the anisotropic and rotated positions if self.int_drift_no > 0: - pos = self.model.anisometrize(pos) - chunk_pos = pos[:, slice(*chunk_slice)] + chunk_pos = self.model.anisometrize(pos)[:, slice(*chunk_slice)] # apply functional drift for i, f in enumerate(self.drift_functions): res[-self.drift_no + i, :] = f(*chunk_pos) From da7e772b47e6fc26924d96d7aa5c53b507b392ec Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Wed, 25 Nov 2020 00:13:53 +0100 Subject: [PATCH 05/39] CovModel: add latlon argument and adopt internal handling; add plot_dim --- gstools/covmodel/base.py | 181 ++++++++++++++++++-------------------- gstools/covmodel/tools.py | 96 +++++++++++++++++++- 2 files changed, 181 insertions(+), 96 deletions(-) diff --git a/gstools/covmodel/base.py b/gstools/covmodel/base.py index 56b5b924..e909636a 100644 --- a/gstools/covmodel/base.py +++ b/gstools/covmodel/base.py @@ -20,8 +20,9 @@ matrix_anisometrize, matrix_isometrize, rotated_main_axes, + latlon2pos, + pos2latlon, ) -from gstools.tools.misc import list_format from gstools.covmodel.tools import ( InitSubclassMeta, _init_subclass, @@ -33,6 +34,8 @@ set_arg_bounds, check_arg_bounds, set_dim, + compare, + model_repr, ) from gstools.covmodel import plot from gstools.covmodel.fit import fit_variogram @@ -89,6 +92,17 @@ class CovModel(metaclass=InitSubclassMeta): to coincide with e.g. the integral scale. Will be set by each model individually. Default: :any:`None` + latlon : :class:`bool`, optional + Whether the model is describing 2D fields on earths surface described + by latitude and longitude. When using this, the model will internally + use the associated 'Yadrenko' model to represent a valid model. + This means, the spatial distance :math:`r` will be replace with + :math:`2\sin(\alpha/2)`, where :math:`\alpha` is the great-circle + distance, which is equal to the spatial distance of two points in 3D. + As a consequence, `dim` will be set to `3` and anisotropy will be + disabled. `rescale` can be set to earths radius, + to have a meaningful `len_scale` parameter. + Default: False var_raw : :class:`float` or :any:`None`, optional raw variance of the model which will be multiplied with :any:`CovModel.var_factor` to result in the actual variance. @@ -116,6 +130,7 @@ def __init__( angles=0.0, integral_scale=None, rescale=None, + latlon=False, var_raw=None, hankel_kw=None, **opt_arg @@ -140,6 +155,8 @@ def __init__( self._nugget_bounds = None self._anis_bounds = None self._opt_arg_bounds = {} + # Set latlon first + self._latlon = bool(latlon) # SFT class will be created within dim.setter but needs hankel_kw self.hankel_kw = hankel_kw self.dim = dim @@ -157,8 +174,14 @@ def __init__( self.default_rescale() if rescale is None else abs(float(rescale)) ) self._nugget = nugget - self._angles = set_angles(self.dim, angles) - self._len_scale, self._anis = set_len_anis(self.dim, len_scale, anis) + # set anisotropy and len_scale, disable anisotropy for latlon models + self._len_scale, anis = set_len_anis(self.dim, len_scale, anis) + if self.latlon: + self._anis = np.array((self.dim - 1) * [1], dtype=np.double) + self._angles = np.array(self.dim * [0], dtype=np.double) + else: + self._anis = anis + self._angles = set_angles(self.dim, angles) # set var at last, because of the var_factor (to be right initialized) if var_raw is None: self._var = None @@ -225,6 +248,18 @@ def cor_axis(self, r, axis=0): return self.correlation(r) return self.correlation(np.abs(r) / self.anis[axis - 1]) + def vario_yadrenko(self, zeta): + r"""Yadrenko variogram for great-circle distance from geo-coords.""" + return self.variogram(2 * np.sin(zeta / 2)) + + def cov_yadrenko(self, zeta): + r"""Yadrenko covariance for great-circle distance from geo-coords.""" + return self.covariance(2 * np.sin(zeta / 2)) + + def cor_yadrenko(self, zeta): + r"""Yadrenko correlation for great-circle distance from geo-coords.""" + return self.correlation(2 * np.sin(zeta / 2)) + def vario_spatial(self, pos): r"""Spatial variogram respecting anisotropy and rotation.""" return self.variogram(self._get_iso_rad(pos)) @@ -270,6 +305,9 @@ def plot(self, func="variogram", **kwargs): # pragma: no cover * "vario_spatial" * "cov_spatial" * "cor_spatial" + * "vario_yadrenko" + * "cov_yadrenko" + * "cor_yadrenko" * "vario_axis" * "cov_axis" * "cor_axis" @@ -493,12 +531,17 @@ def _has_ppf(self): def isometrize(self, pos): """Make a position tuple ready for isotropic operations.""" - pos = np.array(pos, dtype=np.double).reshape((self.dim, -1)) + dim = 2 if self.latlon else self.dim + pos = np.array(pos, dtype=np.double).reshape((dim, -1)) + if self.latlon: + return latlon2pos(pos) return np.dot(matrix_isometrize(self.dim, self.angles, self.anis), pos) def anisometrize(self, pos): """Bring a position tuple into the anisotropic coordinate-system.""" pos = np.array(pos, dtype=np.double).reshape((self.dim, -1)) + if self.latlon: + return pos2latlon(pos) return np.dot( matrix_anisometrize(self.dim, self.angles, self.anis), pos ) @@ -690,11 +733,7 @@ def set_arg_bounds(self, check_args=True, **kwargs): Default: True **kwargs Parameter name as keyword ("var", "len_scale", "nugget", ) - and a list of 2 or 3 values as value: - - * ``[a, b]`` or - * ``[a, b, ]`` - + and a list of 2 or 3 values: ``[a, b]`` or ``[a, b, ]`` where is one of ``"oo"``, ``"cc"``, ``"oc"`` or ``"co"`` to define if the bounds are open ("o") or closed ("c"). """ @@ -712,11 +751,7 @@ def var_bounds(self): Notes ----- - Is a list of 2 or 3 values: - - * ``[a, b]`` or - * ``[a, b, ]`` - + Is a list of 2 or 3 values: ``[a, b]`` or ``[a, b, ]`` where is one of ``"oo"``, ``"cc"``, ``"oc"`` or ``"co"`` to define if the bounds are open ("o") or closed ("c"). """ @@ -726,9 +761,7 @@ def var_bounds(self): def var_bounds(self, bounds): if not check_bounds(bounds): raise ValueError( - "Given bounds for 'var' are not " - + "valid, got: " - + str(bounds) + "Given bounds for 'var' are not valid, got: " + str(bounds) ) self._var_bounds = bounds @@ -738,11 +771,7 @@ def len_scale_bounds(self): Notes ----- - Is a list of 2 or 3 values: - - * ``[a, b]`` or - * ``[a, b, ]`` - + Is a list of 2 or 3 values: ``[a, b]`` or ``[a, b, ]`` where is one of ``"oo"``, ``"cc"``, ``"oc"`` or ``"co"`` to define if the bounds are open ("o") or closed ("c"). """ @@ -752,8 +781,7 @@ def len_scale_bounds(self): def len_scale_bounds(self, bounds): if not check_bounds(bounds): raise ValueError( - "Given bounds for 'len_scale' are not " - + "valid, got: " + "Given bounds for 'len_scale' are not valid, got: " + str(bounds) ) self._len_scale_bounds = bounds @@ -764,11 +792,7 @@ def nugget_bounds(self): Notes ----- - Is a list of 2 or 3 values: - - * ``[a, b]`` or - * ``[a, b, ]`` - + Is a list of 2 or 3 values: ``[a, b]`` or ``[a, b, ]`` where is one of ``"oo"``, ``"cc"``, ``"oc"`` or ``"co"`` to define if the bounds are open ("o") or closed ("c"). """ @@ -778,9 +802,7 @@ def nugget_bounds(self): def nugget_bounds(self, bounds): if not check_bounds(bounds): raise ValueError( - "Given bounds for 'nugget' are not " - + "valid, got: " - + str(bounds) + "Given bounds for 'nugget' are not valid, got: " + str(bounds) ) self._nugget_bounds = bounds @@ -790,11 +812,7 @@ def anis_bounds(self): Notes ----- - Is a list of 2 or 3 values: - - * ``[a, b]`` or - * ``[a, b, ]`` - + Is a list of 2 or 3 values: ``[a, b]`` or ``[a, b, ]`` where is one of ``"oo"``, ``"cc"``, ``"oc"`` or ``"co"`` to define if the bounds are open ("o") or closed ("c"). """ @@ -804,9 +822,7 @@ def anis_bounds(self): def anis_bounds(self, bounds): if not check_bounds(bounds): raise ValueError( - "Given bounds for 'anis' are not " - + "valid, got: " - + str(bounds) + "Given bounds for 'anis' are not valid, got: " + str(bounds) ) self._anis_bounds = bounds @@ -817,10 +833,7 @@ def opt_arg_bounds(self): Notes ----- Keys are the opt-arg names and values are lists of 2 or 3 values: - - * ``[a, b]`` or - * ``[a, b, ]`` - + ``[a, b]`` or ``[a, b, ]`` where is one of ``"oo"``, ``"cc"``, ``"oc"`` or ``"co"`` to define if the bounds are open ("o") or closed ("c"). """ @@ -832,11 +845,8 @@ def arg_bounds(self): Notes ----- - Keys are the opt-arg names and values are lists of 2 or 3 values: - - * ``[a, b]`` or - * ``[a, b, ]`` - + Keys are the arg names and values are lists of 2 or 3 values: + ``[a, b]`` or ``[a, b, ]`` where is one of ``"oo"``, ``"cc"``, ``"oc"`` or ``"co"`` to define if the bounds are open ("o") or closed ("c"). """ @@ -849,6 +859,18 @@ def arg_bounds(self): res.update(self.opt_arg_bounds) return res + # geographical coordinates related + + @property + def latlon(self): + """:class:`bool`: Whether the model depends on geographical coords.""" + return self._latlon + + @property + def plot_dim(self): + """:class:`int`: The plotting dimension of the model.""" + return 2 if self.latlon else self.dim + # standard parameters @property @@ -900,9 +922,11 @@ def len_scale(self): @len_scale.setter def len_scale(self, len_scale): - self._len_scale, self._anis = set_len_anis( - self.dim, len_scale, self.anis - ) + self._len_scale, anis = set_len_anis(self.dim, len_scale, self.anis) + if self.latlon: + self._anis = np.array((self.dim - 1) * [1], dtype=np.double) + else: + self._anis = anis self.check_arg_bounds() @property @@ -922,9 +946,12 @@ def anis(self): @anis.setter def anis(self, anis): - self._len_scale, self._anis = set_len_anis( - self.dim, self.len_scale, anis - ) + if self.latlon: + self._anis = np.array((self.dim - 1) * [1], dtype=np.double) + else: + self._len_scale, self._anis = set_len_anis( + self.dim, self.len_scale, anis + ) self.check_arg_bounds() @property @@ -934,7 +961,10 @@ def angles(self): @angles.setter def angles(self, angles): - self._angles = set_angles(self.dim, angles) + if self.latlon: + self._angles = np.array(self.dim * [0], dtype=np.double) + else: + self._angles = set_angles(self.dim, angles) self.check_arg_bounds() @property @@ -1100,23 +1130,7 @@ def __eq__(self, other): """Compare CovModels.""" if not isinstance(other, CovModel): return False - # prevent attribute error in opt_arg if the are not equal - if set(self.opt_arg) != set(other.opt_arg): - return False - # prevent dim error in anis and angles - if self.dim != other.dim: - return False - equal = True - equal &= self.name == other.name - equal &= np.isclose(self.var, other.var) - equal &= np.isclose(self.var_raw, other.var_raw) # ?! needless? - equal &= np.isclose(self.nugget, other.nugget) - equal &= np.isclose(self.len_scale, other.len_scale) - equal &= np.all(np.isclose(self.anis, other.anis)) - equal &= np.all(np.isclose(self.angles, other.angles)) - for opt in self.opt_arg: - equal &= np.isclose(getattr(self, opt), getattr(other, opt)) - return equal + return compare(self, other) def __ne__(self, other): """Compare CovModels.""" @@ -1128,23 +1142,4 @@ def __str__(self): def __repr__(self): """Return String representation.""" - opt_str = "" - for opt in self.opt_arg: - opt_str += ( - ", " + opt + "={0:.{1}}".format(getattr(self, opt), self._prec) - ) - return ( - "{0}(dim={1}, var={2:.{p}}, len_scale={3:.{p}}, " - "nugget={4:.{p}}, anis={5}, angles={6}".format( - self.name, - self.dim, - self.var, - self.len_scale, - self.nugget, - list_format(self.anis, 3), - list_format(self.angles, 3), - p=self._prec, - ) - + opt_str - + ")" - ) + return model_repr(self) diff --git a/gstools/covmodel/tools.py b/gstools/covmodel/tools.py index 96a5b0ae..70b1ec6b 100644 --- a/gstools/covmodel/tools.py +++ b/gstools/covmodel/tools.py @@ -20,6 +20,8 @@ set_arg_bounds check_arg_bounds set_dim + compare + model_repr """ # pylint: disable=C0103, W0212 @@ -28,6 +30,7 @@ from scipy.optimize import root from scipy import special as sps from hankel import SymmetricFourierTransform as SFT +from gstools.tools.misc import list_format from gstools.tools.geometric import set_anis, set_angles __all__ = [ @@ -44,6 +47,8 @@ "set_arg_bounds", "check_arg_bounds", "set_dim", + "compare", + "model_repr", ] @@ -539,17 +544,25 @@ def set_dim(model, dim): When dimension is < 1. """ # check if a fixed dimension should be used - if model.fix_dim() is not None: + if model.fix_dim() is not None and model.fix_dim() != dim: warnings.warn( model.name + ": using fixed dimension " + str(model.fix_dim()), AttributeWarning, ) dim = model.fix_dim() + if model.latlon and dim != 3: + raise ValueError( + model.name + + ": using fixed dimension " + + str(model.fix_dim()) + + ", which is not compatible with a latlon model." + ) + # force dim=3 for latlon models + dim = 3 if model.latlon else dim # set the dimension if dim < 1: raise ValueError("Only dimensions of d >= 1 are supported.") - check = model.check_dim(dim) - if not check: + if not model.check_dim(dim): warnings.warn( "Dimension {} is not appropriate for this model.".format(dim), AttributeWarning, @@ -565,3 +578,80 @@ def set_dim(model, dim): if model._angles is not None: model._angles = set_angles(model.dim, model._angles) model.check_arg_bounds() + + +def compare(this, that): + """ + Compare CovModels. + + Parameters + ---------- + this / that : :any:`CovModel` + The covariance models to compare. + """ + # prevent attribute error in opt_arg if the are not equal + if set(this.opt_arg) != set(that.opt_arg): + return False + # prevent dim error in anis and angles + if this.dim != that.dim: + return False + equal = True + equal &= this.name == that.name + equal &= np.isclose(this.var, that.var) + equal &= np.isclose(this.var_raw, that.var_raw) # ?! needless? + equal &= np.isclose(this.nugget, that.nugget) + equal &= np.isclose(this.len_scale, that.len_scale) + equal &= np.all(np.isclose(this.anis, that.anis)) + equal &= np.all(np.isclose(this.angles, that.angles)) + equal &= np.isclose(this.rescale, that.rescale) + equal &= this.latlon == that.latlon + for opt in this.opt_arg: + equal &= np.isclose(getattr(this, opt), getattr(that, opt)) + return equal + + +def model_repr(model): + """ + Generate the model string representation. + + Parameters + ---------- + model : :any:`CovModel` + The covariance model in use. + """ + opt_str = "" + for opt in model.opt_arg: + opt_str += ( + ", " + opt + "={0:.{1}}".format(getattr(model, opt), model._prec) + ) + if model.latlon: + std_str = ( + "{0}(latlon={1}, var={2:.{p}}, len_scale={3:.{p}}, " + "nugget={4:.{p}}".format( + model.name, + model.latlon, + model.var, + model.len_scale, + model.nugget, + p=model._prec, + ) + + opt_str + + ")" + ) + else: + std_str = ( + "{0}(dim={1}, var={2:.{p}}, len_scale={3:.{p}}, " + "nugget={4:.{p}}, anis={5}, angles={6}".format( + model.name, + model.dim, + model.var, + model.len_scale, + model.nugget, + list_format(model.anis, 3), + list_format(model.angles, 3), + p=model._prec, + ) + + opt_str + + ")" + ) + return std_str From a1170a2eba5a059a47c65ff8ce6c060c893fbf92 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Wed, 25 Nov 2020 00:14:48 +0100 Subject: [PATCH 06/39] CovModel.fit: allow geographical estimated variograms for fitting --- gstools/covmodel/fit.py | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/gstools/covmodel/fit.py b/gstools/covmodel/fit.py index 32056063..4f1121a6 100755 --- a/gstools/covmodel/fit.py +++ b/gstools/covmodel/fit.py @@ -286,6 +286,13 @@ def _check_vario(model, x_data, y_data): + "Either provide only one variogram to fit an isotropic model, " + "or directional ones for all main axes to fit anisotropy." ) + if is_dir_vario and model.latlon: + raise ValueError( + "CovModel.fit_variogram: lat-lon models don't support anisotropy." + ) + elif model.latlon: + # convert to yadrenko model + x_data = 2 * np.sin(x_data / 2) return x_data, y_data, is_dir_vario From 7b4359d492d9e0aa9f73936d557e91b4ba8494d8 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Wed, 25 Nov 2020 00:15:17 +0100 Subject: [PATCH 07/39] CovModel.plot: add ploting routines for yadrenko methods --- gstools/covmodel/plot.py | 51 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 51 insertions(+) diff --git a/gstools/covmodel/plot.py b/gstools/covmodel/plot.py index 4cd4aca7..e17b8643 100644 --- a/gstools/covmodel/plot.py +++ b/gstools/covmodel/plot.py @@ -10,6 +10,9 @@ plot_variogram plot_covariance plot_correlation + plot_vario_yadrenko + plot_cov_yadrenko + plot_cor_yadrenko plot_vario_axis plot_cov_axis plot_cor_axis @@ -30,6 +33,9 @@ "plot_variogram", "plot_covariance", "plot_correlation", + "plot_vario_yadrenko", + "plot_cov_yadrenko", + "plot_cor_yadrenko", "plot_vario_axis", "plot_cov_axis", "plot_cor_axis", @@ -151,6 +157,51 @@ def plot_correlation( return ax +def plot_vario_yadrenko( + model, x_min=0.0, x_max=None, fig=None, ax=None, **kwargs +): # pragma: no cover + """Plot Yadrenko variogram of a given CovModel.""" + fig, ax = _get_fig_ax(fig, ax) + if x_max is None: + x_max = min(3 * model.len_rescaled, np.pi) + x_s = np.linspace(x_min, x_max) + kwargs.setdefault("label", model.name + " Yadrenko variogram") + ax.plot(x_s, model.vario_yadrenko(x_s), **kwargs) + ax.legend() + fig.show() + return ax + + +def plot_cov_yadrenko( + model, x_min=0.0, x_max=None, fig=None, ax=None, **kwargs +): # pragma: no cover + """Plot Yadrenko covariance of a given CovModel.""" + fig, ax = _get_fig_ax(fig, ax) + if x_max is None: + x_max = min(3 * model.len_rescaled, np.pi) + x_s = np.linspace(x_min, x_max) + kwargs.setdefault("label", model.name + " Yadrenko covariance") + ax.plot(x_s, model.cov_yadrenko(x_s), **kwargs) + ax.legend() + fig.show() + return ax + + +def plot_cor_yadrenko( + model, x_min=0.0, x_max=None, fig=None, ax=None, **kwargs +): # pragma: no cover + """Plot Yadrenko correlation function of a given CovModel.""" + fig, ax = _get_fig_ax(fig, ax) + if x_max is None: + x_max = min(3 * model.len_rescaled, np.pi) + x_s = np.linspace(x_min, x_max) + kwargs.setdefault("label", model.name + " Yadrenko correlation") + ax.plot(x_s, model.cor_yadrenko(x_s), **kwargs) + ax.legend() + fig.show() + return ax + + def plot_vario_axis( model, axis=0, x_min=0.0, x_max=None, fig=None, ax=None, **kwargs ): # pragma: no cover From 2a3d14fe96c62be1f6686ad3ddfd8dce0f5c7d21 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Wed, 25 Nov 2020 00:16:37 +0100 Subject: [PATCH 08/39] Field: correct dimension handling for latlon fields --- gstools/field/base.py | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/gstools/field/base.py b/gstools/field/base.py index cf74c73f..937e7160 100755 --- a/gstools/field/base.py +++ b/gstools/field/base.py @@ -181,18 +181,21 @@ def pre_pos(self, pos, mesh_type="unstructured"): """ # save mesh-type self.mesh_type = mesh_type + dim = 2 if self.model.latlon else self.model.dim # save pos tuple if mesh_type != "unstructured": - pos, shape = format_struct_pos_dim(pos, self.model.dim) + pos, shape = format_struct_pos_dim(pos, dim) self.pos = pos pos = gen_mesh(pos) else: - pos = np.array(pos, dtype=np.double).reshape(self.model.dim, -1) + pos = np.array(pos, dtype=np.double).reshape(dim, -1) self.pos = pos shape = np.shape(pos[0]) # prepend dimension if we have a vector field if self.value_type == "vector": shape = (self.model.dim,) + shape + if self.model.latlon: + raise ValueError("Field: Vector fields not allowed for latlon") # return isometrized pos tuple and resulting field shape return self.model.isometrize(pos), shape From da723138a07c22ecbe88cb9dda3796f38e99ec47 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Wed, 25 Nov 2020 00:17:10 +0100 Subject: [PATCH 09/39] Field.plot: add plotting routines for latlon fields --- gstools/field/plot.py | 35 ++++++++++++++++++++++------------- 1 file changed, 22 insertions(+), 13 deletions(-) diff --git a/gstools/field/plot.py b/gstools/field/plot.py index 29a77f81..87c7efc3 100644 --- a/gstools/field/plot.py +++ b/gstools/field/plot.py @@ -42,11 +42,13 @@ def plot_field(fld, field="field", fig=None, ax=None): # pragma: no cover """ plot_field = getattr(fld, field) assert not (fld.pos is None or plot_field is None) - if fld.model.dim == 1: + if fld.model.plot_dim == 1: ax = _plot_1d(fld.pos, plot_field, fig, ax) - elif fld.model.dim == 2: - ax = _plot_2d(fld.pos, plot_field, fld.mesh_type, fig, ax) - elif fld.model.dim == 3: + elif fld.model.plot_dim == 2: + ax = _plot_2d( + fld.pos, plot_field, fld.mesh_type, fig, ax, fld.model.latlon + ) + elif fld.model.plot_dim == 3: ax = _plot_3d(fld.pos, plot_field, fld.mesh_type, fig, ax) else: raise ValueError("Field.plot: only possible for dim=1,2,3!") @@ -68,21 +70,28 @@ def _plot_1d(pos, field, fig=None, ax=None): # pragma: no cover return ax -def _plot_2d(pos, field, mesh_type, fig=None, ax=None): # pragma: no cover +def _plot_2d( + pos, field, mesh_type, fig=None, ax=None, latlon=False +): # pragma: no cover """Plot a 2d field.""" fig, ax = _get_fig_ax(fig, ax) title = "Field 2D " + mesh_type + ": " + str(field.shape) - x = pos[0] - y = pos[1] + y = pos[0] if latlon else pos[1] + x = pos[1] if latlon else pos[0] if mesh_type == "unstructured": cont = ax.tricontourf(x, y, field.ravel(), levels=256) else: + plot_field = field if latlon else field.T try: - cont = ax.contourf(x, y, field.T, levels=256) + cont = ax.contourf(x, y, plot_field, levels=256) except TypeError: - cont = ax.contourf(x, y, field.T, 256) - ax.set_xlabel("X") - ax.set_ylabel("Y") + cont = ax.contourf(x, y, plot_field, 256) + if latlon: + ax.set_ylabel("Lat in deg") + ax.set_xlabel("Lon in deg") + else: + ax.set_xlabel("X") + ax.set_ylabel("Y") ax.set_title(title) fig.colorbar(cont) fig.show() @@ -147,7 +156,7 @@ def get_plane(z_val_in, z_dir): z_io = dir2 * z_range + z_min x_io = np.full_like(y_io, z_val_in) - if mesh_type == "structured": + if mesh_type != "unstructured": # contourf plots image like for griddata, therefore transpose plane = inter.interpn( pos, field, np.array((x_io, y_io, z_io)).T, bounds_error=False @@ -231,7 +240,7 @@ def plot_vec_field(fld, field="field", fig=None, ax=None): # pragma: no cover Axes to plot on. If `None`, a new one will be added to the figure. Default: `None` """ - if fld.mesh_type != "structured": + if fld.mesh_type == "unstructured": raise RuntimeError( "Only structured vector fields are supported" + " for plotting. Please create one on a structured grid." From b04314082f7495d4f212f962da4ef3948b2b95ae Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Wed, 25 Nov 2020 00:18:43 +0100 Subject: [PATCH 10/39] Examples: add geo coordinates tutorial --- .../08_geo_coordinates/00_field_generation.py | 63 ++++++++++++++++ examples/08_geo_coordinates/README.rst | 71 +++++++++++++++++++ 2 files changed, 134 insertions(+) create mode 100755 examples/08_geo_coordinates/00_field_generation.py create mode 100644 examples/08_geo_coordinates/README.rst diff --git a/examples/08_geo_coordinates/00_field_generation.py b/examples/08_geo_coordinates/00_field_generation.py new file mode 100755 index 00000000..bb3e787e --- /dev/null +++ b/examples/08_geo_coordinates/00_field_generation.py @@ -0,0 +1,63 @@ +""" +Working with lat-lon random fields +---------------------------------- + +In this example, we demonstrate how to generate a random field on +geographical coordinates. + +First we setup a model, with ``latlon=True``, to get the associated +Yadrenko model. + +In addition we will use the earth radius provided by :any:`EARTH_RADIUS`, +to have a meaningful length scale in km. + +To generate the field, we simply pass ``(lat, lon)`` as position tuple +to the :any:`SRF` class. +""" +# sphinx_gallery_thumbnail_number = 3 +import gstools as gs + +model = gs.Gaussian(latlon=True, var=1, len_scale=777, rescale=gs.EARTH_RADIUS) + +lat = lon = range(-80, 81) +srf = gs.SRF(model, seed=12345) +field = srf.structured((lat, lon)) +srf.plot() + +############################################################################### +# This was easy as always! Now we can use this field to estimate the empirical +# variogram in order to prove, that the generated field has the correct +# geo-statistical properties. +# The :any:`vario_estimate` routine also provides a ``latlon`` switch to +# indicate, that the given field is defined on geographical variables. +# +# As we will see, everthing went well... Phew! + +bin_edges = [0.01 * i for i in range(30)] +bin_center, emp_vario = gs.vario_estimate( + *((lat, lon), field, bin_edges), + latlon=True, + mesh_type="structured", + sampling_size=2000, + sampling_seed=12345, +) + +ax = model.plot("vario_yadrenko", x_max=0.3) +model.fit_variogram(bin_center, emp_vario, init_guess="current", nugget=False) +model.plot("vario_yadrenko", ax=ax, label="fitted", x_max=0.3) +ax.scatter(bin_center, emp_vario, color="k") +print(model) + +############################################################################### +# The resulting field can also be easily visualized with the aid of +# `cartopy `_. + +import matplotlib.pyplot as plt +import cartopy.crs as ccrs + +fig, ax = plt.subplots(subplot_kw={"projection": ccrs.Orthographic(-45, 45)}) +cont = ax.contourf(lon, lat, field, transform=ccrs.PlateCarree()) +ax.set_title("lat-lon random field generated with GSTools") +ax.coastlines() +ax.set_global() +fig.colorbar(cont) diff --git a/examples/08_geo_coordinates/README.rst b/examples/08_geo_coordinates/README.rst new file mode 100644 index 00000000..6655bf2d --- /dev/null +++ b/examples/08_geo_coordinates/README.rst @@ -0,0 +1,71 @@ +Tutorial 8: Geographic Coordinates +================================== + +GSTools provides support for +`geographic coordinates `_ +given by: + +- latitude ``lat``: specifies the north–south position of a point on the Earth's surface +- longitude ``lon``: specifies the east–west position of a point on the Earth's surface + +If you want to use this feature for field generation or Kriging, you +have to setup a geographical covariance Model by setting ``latlon=True`` +in your desired model (see :any:`CovModel`): + +.. code-block:: python + + import numpy as np + import gstools as gs + + model = gs.Gaussian(latlon=True, var=2, len_scale=np.pi / 16) + +By doing so, the model will use the associated `Yadrenko` model on a sphere +(see `here `_). +The `len_scale` is given in radians to scale the arc-length. +In order to have a more meaningful length scale, one can use the ``rescale`` +argument: + +.. code-block:: python + + import gstools as gs + + model = gs.Gaussian(latlon=True, var=2, len_scale=500, rescale=gs.EARTH_RADIUS) + +Then ``len_scale`` can be interpreted as given in km. + +A `Yadrenko` model :math:`C` is derived from a valid +isotropic covariance model in 3D :math:`C_{3D}` by the following relation: + +.. math:: + C(\zeta)=C_{3D}\left(2 \cdot \sin\left(\frac{\zeta}{2}\right)\right) + +Where :math:`\zeta` is the +`great-circle distance `_. + +.. note:: + + ``lat`` and ``lon`` are given in degree, whereas the great-circle distance + :math:`zeta` is given in radians. + +Note, that :math:`2 \cdot \sin(\frac{\zeta}{2})` is the +`chordal distance `_ +of two points on a sphere, which means we simply think of the earth surface +as a sphere, that is cut out of the surrounding three dimensional space, +when using the `Yadrenko` model. + +.. note:: + + Anisotropy is not available with the geographical models, since their + geometry is not euclidean. When passing values for :any:`CovModel.anis` + or :any:`CovModel.angles`, they will be ignored. + + Since the Yadrenko model comes from a 3D model, the model dimension will + be 3 (see :any:`CovModel.dim`) but the `plot_dim` will be 2 in this case + (see :any:`CovModel.plot_dim`). + +.. only:: html + + Gallery + ------- + + Below is a gallery of examples From a6dd3d3f23949bed37d1eb6e9cc832e7db7e338c Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Wed, 25 Nov 2020 00:18:56 +0100 Subject: [PATCH 11/39] Examples: minor fixes --- examples/01_random_field/04_srf_merge.py | 1 + examples/02_cov_model/README.rst | 5 ++++- 2 files changed, 5 insertions(+), 1 deletion(-) diff --git a/examples/01_random_field/04_srf_merge.py b/examples/01_random_field/04_srf_merge.py index 366558da..13890280 100644 --- 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/02_cov_model/README.rst b/examples/02_cov_model/README.rst index 7c9c413c..90811b5c 100644 --- a/examples/02_cov_model/README.rst +++ b/examples/02_cov_model/README.rst @@ -70,6 +70,7 @@ 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 @@ -77,7 +78,9 @@ As a special feature, we also provide truncated power law (TPL) covariance model TPLGaussian TPLExponential TPLStable - TPLSimple + +These models provide a lower and upper length scale truncation +for superpositioned models. .. only:: html From 074f2229c4cef3effef9406b00e2cc836f31e232 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Wed, 25 Nov 2020 00:19:24 +0100 Subject: [PATCH 12/39] Doc: add new examples to doc --- docs/requirements_doc.txt | 1 + docs/source/conf.py | 3 +++ docs/source/tutorials.rst | 1 + 3 files changed, 5 insertions(+) diff --git a/docs/requirements_doc.txt b/docs/requirements_doc.txt index cdaea2ce..fb32d98b 100755 --- a/docs/requirements_doc.txt +++ b/docs/requirements_doc.txt @@ -4,3 +4,4 @@ matplotlib pyvista pykrige meshzoo +cartopy diff --git a/docs/source/conf.py b/docs/source/conf.py index 25b2dc41..7b856dd0 100644 --- 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/tutorials.rst b/docs/source/tutorials.rst index 5c2b6787..fab93145 100644 --- 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 From 5779bbd0492cd503d02b085b9b2d0297e2b7e8a3 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Wed, 25 Nov 2020 12:40:28 +0100 Subject: [PATCH 13/39] Examples: remove cartopy plot (not working with RTD) --- .../08_geo_coordinates/00_field_generation.py | 20 +++++++------------ 1 file changed, 7 insertions(+), 13 deletions(-) diff --git a/examples/08_geo_coordinates/00_field_generation.py b/examples/08_geo_coordinates/00_field_generation.py index bb3e787e..46406273 100755 --- a/examples/08_geo_coordinates/00_field_generation.py +++ b/examples/08_geo_coordinates/00_field_generation.py @@ -14,7 +14,6 @@ To generate the field, we simply pass ``(lat, lon)`` as position tuple to the :any:`SRF` class. """ -# sphinx_gallery_thumbnail_number = 3 import gstools as gs model = gs.Gaussian(latlon=True, var=1, len_scale=777, rescale=gs.EARTH_RADIUS) @@ -49,15 +48,10 @@ print(model) ############################################################################### -# The resulting field can also be easily visualized with the aid of -# `cartopy `_. - -import matplotlib.pyplot as plt -import cartopy.crs as ccrs - -fig, ax = plt.subplots(subplot_kw={"projection": ccrs.Orthographic(-45, 45)}) -cont = ax.contourf(lon, lat, field, transform=ccrs.PlateCarree()) -ax.set_title("lat-lon random field generated with GSTools") -ax.coastlines() -ax.set_global() -fig.colorbar(cont) +# .. note:: +# +# Note, that the estimated variogram coincides with the yadrenko variogram, +# which means it depends on the great-circle distance. +# +# Keep that in mind when defining bins: The range is at most +# :math:`\pi\approx 3.14`. From 17cfc466d9e42b95cf3af4ca2529be4bb0a10afd Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Wed, 25 Nov 2020 12:41:11 +0100 Subject: [PATCH 14/39] Docs: move cartopy plot to README; update pyvista plot; add new tutorials section in README --- README.md | 32 +++++++++++++++++++++++++++++--- docs/requirements_doc.txt | 1 - docs/source/index.rst | 31 ++++++++++++++++++++++++++++--- 3 files changed, 57 insertions(+), 7 deletions(-) diff --git a/README.md b/README.md index 15de8d9f..d51839b9 100644 --- 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() Random field

+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) +``` + +

+lat-lon random field +

+ 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() ```

-3d Random field +3d Random field

@@ -335,6 +360,7 @@ You can contact us via . [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/requirements_doc.txt b/docs/requirements_doc.txt index fb32d98b..cdaea2ce 100755 --- a/docs/requirements_doc.txt +++ b/docs/requirements_doc.txt @@ -4,4 +4,3 @@ matplotlib pyvista pykrige meshzoo -cartopy diff --git a/docs/source/index.rst b/docs/source/index.rst index f76f299f..7ca91ae0 100644 --- 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 `__ - `Conditioned random field generation `__ - `Field transformations `__ +- `Field transformations `__ - `Miscellaneous examples `__ 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 `_. +This works perfectly well with `cartopy `_. + +.. 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 `__ file, which can be visualized with `ParaView `_ 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 From bf2efc3cc9f31940712275c75b90a68bcaef9aaa Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Wed, 25 Nov 2020 20:11:29 +0100 Subject: [PATCH 15/39] Docs: minor changes for latlon --- docs/source/index.rst | 2 +- examples/08_geo_coordinates/00_field_generation.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/docs/source/index.rst b/docs/source/index.rst index 7ca91ae0..296de012 100644 --- a/docs/source/index.rst +++ b/docs/source/index.rst @@ -97,7 +97,7 @@ showing the most important use cases of GSTools, which are - `Kriging `__ - `Conditioned random field generation `__ - `Field transformations `__ -- `Field transformations `__ +- `Geographic Coordinates `__ - `Miscellaneous examples `__ Some more examples are provided in the examples folder. diff --git a/examples/08_geo_coordinates/00_field_generation.py b/examples/08_geo_coordinates/00_field_generation.py index 46406273..91583075 100755 --- a/examples/08_geo_coordinates/00_field_generation.py +++ b/examples/08_geo_coordinates/00_field_generation.py @@ -54,4 +54,4 @@ # which means it depends on the great-circle distance. # # Keep that in mind when defining bins: The range is at most -# :math:`\pi\approx 3.14`. +# :math:`\pi\approx 3.14`, which corresponds to the half globe. From 5cdc80ecff2ddba170d8f2db7d781924e00c8592 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Wed, 25 Nov 2020 20:11:55 +0100 Subject: [PATCH 16/39] Docs: workaround for sphinx-gallery bug with html only directive --- examples/00_misc/README.rst | 8 ++------ examples/01_random_field/README.rst | 8 ++------ examples/02_cov_model/README.rst | 8 ++------ examples/03_variogram/README.rst | 8 ++------ examples/04_vector_field/README.rst | 7 ++----- examples/05_kriging/README.rst | 8 ++------ examples/06_conditioned_fields/README.rst | 8 ++------ examples/07_transformations/README.rst | 8 ++------ examples/08_geo_coordinates/README.rst | 8 ++------ 9 files changed, 18 insertions(+), 53 deletions(-) diff --git a/examples/00_misc/README.rst b/examples/00_misc/README.rst index 15492cbd..69d8fc0d 100644 --- a/examples/00_misc/README.rst +++ b/examples/00_misc/README.rst @@ -3,9 +3,5 @@ Miscellaneous A few miscellaneous examples -.. only:: html - - Gallery - ------- - - Below is a gallery of examples +Examples +-------- diff --git a/examples/01_random_field/README.rst b/examples/01_random_field/README.rst index 0682d23e..de516df9 100644 --- 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/02_cov_model/README.rst b/examples/02_cov_model/README.rst index 90811b5c..39a79f27 100644 --- a/examples/02_cov_model/README.rst +++ b/examples/02_cov_model/README.rst @@ -82,9 +82,5 @@ As a special feature, we also provide truncated power law (TPL) covariance model These models provide a lower and upper length scale truncation for superpositioned models. -.. only:: html - - Gallery - ------- - - Below is a gallery of examples +Examples +-------- diff --git a/examples/03_variogram/README.rst b/examples/03_variogram/README.rst index 5ba19bd7..a0380df6 100644 --- a/examples/03_variogram/README.rst +++ b/examples/03_variogram/README.rst @@ -10,9 +10,5 @@ The same `(semi-)variogram Date: Wed, 25 Nov 2020 20:12:30 +0100 Subject: [PATCH 17/39] Variogram: 'h' as flag for haversine --- gstools/variogram/variogram.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/gstools/variogram/variogram.py b/gstools/variogram/variogram.py index c6a2249d..6f63b336 100644 --- a/gstools/variogram/variogram.py +++ b/gstools/variogram/variogram.py @@ -293,7 +293,8 @@ def vario_estimate( cython_estimator = _set_estimator(estimator) # run if dir_no == 0: - distance_type = "s" if latlon else "e" + # "h"aversine or "e"uclidean distance type + distance_type = "h" if latlon else "e" estimates, counts = unstructured( dim, field, From 4fcc0c80680c6838d5a031c6b71e6cc21a169ca8 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Wed, 25 Nov 2020 20:30:19 +0100 Subject: [PATCH 18/39] CovModel: rename 'plot_dim' to 'field_dim' --- examples/08_geo_coordinates/README.rst | 4 ++-- gstools/covmodel/base.py | 7 +++---- gstools/field/base.py | 2 +- gstools/field/plot.py | 6 +++--- 4 files changed, 9 insertions(+), 10 deletions(-) diff --git a/examples/08_geo_coordinates/README.rst b/examples/08_geo_coordinates/README.rst index 5059bc21..0fdaeb87 100644 --- a/examples/08_geo_coordinates/README.rst +++ b/examples/08_geo_coordinates/README.rst @@ -60,8 +60,8 @@ when using the `Yadrenko` model. or :any:`CovModel.angles`, they will be ignored. Since the Yadrenko model comes from a 3D model, the model dimension will - be 3 (see :any:`CovModel.dim`) but the `plot_dim` will be 2 in this case - (see :any:`CovModel.plot_dim`). + be 3 (see :any:`CovModel.dim`) but the `field_dim` will be 2 in this case + (see :any:`CovModel.field_dim`). Examples -------- diff --git a/gstools/covmodel/base.py b/gstools/covmodel/base.py index e909636a..7dd0c9f8 100644 --- a/gstools/covmodel/base.py +++ b/gstools/covmodel/base.py @@ -531,8 +531,7 @@ def _has_ppf(self): def isometrize(self, pos): """Make a position tuple ready for isotropic operations.""" - dim = 2 if self.latlon else self.dim - pos = np.array(pos, dtype=np.double).reshape((dim, -1)) + pos = np.array(pos, dtype=np.double).reshape((self.field_dim, -1)) if self.latlon: return latlon2pos(pos) return np.dot(matrix_isometrize(self.dim, self.angles, self.anis), pos) @@ -867,8 +866,8 @@ def latlon(self): return self._latlon @property - def plot_dim(self): - """:class:`int`: The plotting dimension of the model.""" + def field_dim(self): + """:class:`int`: The field dimension of the model.""" return 2 if self.latlon else self.dim # standard parameters diff --git a/gstools/field/base.py b/gstools/field/base.py index 937e7160..9ddb8e35 100755 --- a/gstools/field/base.py +++ b/gstools/field/base.py @@ -181,7 +181,7 @@ def pre_pos(self, pos, mesh_type="unstructured"): """ # save mesh-type self.mesh_type = mesh_type - dim = 2 if self.model.latlon else self.model.dim + dim = self.model.field_dim # save pos tuple if mesh_type != "unstructured": pos, shape = format_struct_pos_dim(pos, dim) diff --git a/gstools/field/plot.py b/gstools/field/plot.py index 87c7efc3..edfbd6ee 100644 --- a/gstools/field/plot.py +++ b/gstools/field/plot.py @@ -42,13 +42,13 @@ def plot_field(fld, field="field", fig=None, ax=None): # pragma: no cover """ plot_field = getattr(fld, field) assert not (fld.pos is None or plot_field is None) - if fld.model.plot_dim == 1: + if fld.model.field_dim == 1: ax = _plot_1d(fld.pos, plot_field, fig, ax) - elif fld.model.plot_dim == 2: + elif fld.model.field_dim == 2: ax = _plot_2d( fld.pos, plot_field, fld.mesh_type, fig, ax, fld.model.latlon ) - elif fld.model.plot_dim == 3: + elif fld.model.field_dim == 3: ax = _plot_3d(fld.pos, plot_field, fld.mesh_type, fig, ax) else: raise ValueError("Field.plot: only possible for dim=1,2,3!") From f82b5cb1db7ea386d5ae9d07943b7cff2fd50916 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Fri, 27 Nov 2020 02:15:33 +0100 Subject: [PATCH 19/39] CovModel.fit_variogram: use rescale factor as default init guess for len_scale --- gstools/covmodel/fit.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/gstools/covmodel/fit.py b/gstools/covmodel/fit.py index 4f1121a6..0fa17075 100755 --- a/gstools/covmodel/fit.py +++ b/gstools/covmodel/fit.py @@ -328,7 +328,7 @@ def _init_curve_fit_para(model, para, init_guess, constrain_sill, sill, anis): _init_guess( bounds=[low_bounds[-1], top_bounds[-1]], current=getattr(model, par), - default=1.0, + default=model.rescale if par == "len_scale" else 1.0, typ=init_guess, para_name=par, ) From 7ac7262cdd30f5ae6e4b8e52e3b2f2639c8e30e6 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Fri, 27 Nov 2020 02:16:01 +0100 Subject: [PATCH 20/39] Krige: use correct dimension for drifts and cond_pos values when dealing with lat-lon data --- gstools/krige/base.py | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/gstools/krige/base.py b/gstools/krige/base.py index 4fab9d23..591e38e9 100755 --- a/gstools/krige/base.py +++ b/gstools/krige/base.py @@ -331,7 +331,10 @@ def _post_field(self, field, krige_var): self.field = field else: self.field = field + eval_func( - self.trend_function, self.pos, self.model.dim, self.mesh_type + self.trend_function, + self.pos, + self.model.field_dim, + self.mesh_type, ) self.field += self.mean self.krige_var = self.model.sill - krige_var @@ -402,7 +405,7 @@ def set_condition( """ # correctly format cond_pos and cond_val self._cond_pos, self._cond_val = set_condition( - cond_pos, cond_val, self.model.dim + cond_pos, cond_val, self.model.field_dim ) # set the measurement errors self.cond_err = cond_err @@ -435,7 +438,7 @@ def set_drift_functions(self, drift_functions=None): self._drift_functions = [] elif isinstance(drift_functions, (str, int)): self._drift_functions = get_drift_functions( - self.model.dim, drift_functions + self.model.field_dim, drift_functions ) else: if isinstance(drift_functions, collections.abc.Iterator): From 5e238727a953258588cfbd5c3f49ddb00bee0b8d Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Fri, 27 Nov 2020 02:16:29 +0100 Subject: [PATCH 21/39] Examples: add kriging with lat-lon example with DWD data --- examples/08_geo_coordinates/01_dwd_krige.py | 170 +++++++ examples/08_geo_coordinates/de_borders.txt | 492 +++++++++++++++++++ examples/08_geo_coordinates/temp_obs.txt | 494 ++++++++++++++++++++ 3 files changed, 1156 insertions(+) create mode 100755 examples/08_geo_coordinates/01_dwd_krige.py create mode 100644 examples/08_geo_coordinates/de_borders.txt create mode 100644 examples/08_geo_coordinates/temp_obs.txt diff --git a/examples/08_geo_coordinates/01_dwd_krige.py b/examples/08_geo_coordinates/01_dwd_krige.py new file mode 100755 index 00000000..f4dfb6fe --- /dev/null +++ b/examples/08_geo_coordinates/01_dwd_krige.py @@ -0,0 +1,170 @@ +""" +Kriging geographical data +------------------------- + +In this example we are going to interpolate actual temperature data from +the german weather service `DWD `_. + +Data is retrieved utilizing the beatiful package +`wetterdienst `_, +which serves as an API for the DWD data. + +For better visualization, we also download a simple shapefile of the german +borderline with `cartopy `_. + +In order to hold the number of dependecies low, the calls of both functions +shown beneath are commented out. +""" +# sphinx_gallery_thumbnail_number = 2 +import numpy as np +import matplotlib.pyplot as plt +import gstools as gs + + +def get_borders_germany(): + """Download simple german shape file with cartopy.""" + from cartopy.io import shapereader as shp_read # version 0.18.0 + import geopandas as gp # 0.8.1 + + shpfile = shp_read.natural_earth("50m", "cultural", "admin_0_countries") + df = gp.read_file(shpfile) # only use the simples polygon + poly = df.loc[df["ADMIN"] == "Germany"]["geometry"].values[0][0] + np.savetxt("de_borders.txt", list(poly.exterior.coords)) + + +def get_dwd_temperature(): + """Get air temperature from german weather stations from 9.6.20 12:00.""" + from wetterdienst.dwd import observations as obs # version 0.10.1 + + sites = obs.DWDObservationSites( + parameter_set=obs.DWDObservationParameterSet.TEMPERATURE_AIR, + resolution=obs.DWDObservationResolution.HOURLY, + period=obs.DWDObservationPeriod.RECENT, + start_date="2020-06-09 12:00:00", + end_date="2020-06-09 12:00:00", + ) + df0 = sites.all() + ids, lat, lon = map(np.array, [df0.STATION_ID, df0.LAT, df0.LON]) + observations = obs.DWDObservationData( + station_ids=ids, + parameters=obs.DWDObservationParameter.HOURLY.TEMPERATURE_AIR_200, + resolution=obs.DWDObservationResolution.HOURLY, + start_date="2020-06-09 12:00:00", + end_date="2020-06-09 12:00:00", + ) + df1 = observations.collect_safe() + temp, ids1 = map(np.array, [df1.VALUE, df1.STATION_ID]) + select = np.isfinite(temp) # care about missing values + sorter = np.argsort(ids) # care about missing stations + sort = sorter[np.searchsorted(ids, ids1[select], sorter=np.argsort(ids))] + ids, lat, lon, temp = ids[sort], lat[sort], lon[sort], temp[select] + head = "id, lat, lon, temp" # add a header to the file + np.savetxt("temp_obs.txt", np.array([ids, lat, lon, temp]).T, header=head) + + +############################################################################### +# If you want to download the data again, +# uncomment the two following lines. We will simply load the resulting +# files to gain the border polygon and the observed temperature along with +# the station locations given by lat-lon values. + +# get_borders_germany() +# get_dwd_temperature() + +border = np.loadtxt("de_borders.txt") +ids, lat, lon, temp = np.loadtxt("temp_obs.txt").T + +############################################################################### +# First we will estimate the variogram of our temperature data. +# As the maximal bin distance we choose 8 degrees, which corresponds to a +# chordal length of about 900 km. + +bin_max = np.deg2rad(8) +bins = np.linspace(0, bin_max, 20) +bin_c, vario = gs.vario_estimate((lat, lon), temp, bins, latlon=True) + +############################################################################### +# Now we can use this estimated variogram to fit a model to it. +# Here we will use a :any:`Spherical` model. We select the ``latlon`` option +# to use the `Yadrenko` variant of the model to gain a vaild model for lat-lon +# coordinates and we rescale it to the earth-radius. Otherwise the length +# scale would be given in radians representing the great-circle distance. +# +# We deselect the nugget from fitting and plot the result afterwards. +# +# .. note:: +# +# You need to plot the yadrenko variogram, since the standard variogram +# still holds the ordinary routine that is not respecting the great-circle +# distance. + +model = gs.Spherical(latlon=True, rescale=gs.EARTH_RADIUS) +model.fit_variogram(bin_c, vario, nugget=False) +ax = model.plot("vario_yadrenko", x_max=bin_max) +ax.scatter(bin_c, vario) +print(model) + +############################################################################### +# As we see, we have a rather large correlation length of 600 km. +# +# Now we want to interpolate the data using :any:`Universal` kriging. +# In order to tinker around with the data, we will use a north-south drift +# by assuming a linear correlation to the latitude. +# This can be done as follows: + +north_south_drift = lambda lat, lon: lat + +uk = gs.krige.Universal( + model=model, + cond_pos=(lat, lon), + cond_val=temp, + drift_functions=north_south_drift, +) + +############################################################################### +# Now we generate the kriging field, by defining a lat-lon grid that covers +# whole germany. The :any:`Krige` class provides the option to only krige the +# mean field, so one can have a glimpse at the estimated drift. + +g_lat = np.arange(47, 56.1, 0.1) +g_lon = np.arange(5, 16.1, 0.1) + +field, k_var = uk((g_lat, g_lon), mesh_type="structured") +mean, m_var = uk((g_lat, g_lon), mesh_type="structured", only_mean=True) + +############################################################################### +# And that's it. Now let's have a look at the generated field and the input +# data along with the estimated mean: + +levels = np.linspace(5, 23, 64) +fig, ax = plt.subplots(1, 3, figsize=[10, 5], sharey=True) +sca = ax[0].scatter(lon, lat, c=temp, vmin=5, vmax=23, cmap="coolwarm") +co1 = ax[1].contourf(g_lon, g_lat, field, levels, cmap="coolwarm") +co2 = ax[2].contourf(g_lon, g_lat, mean, levels, cmap="coolwarm") + +[ax[i].plot(border[:, 0], border[:, 1], color="k") for i in range(3)] +[ax[i].set_xlim([5, 16]) for i in range(3)] +[ax[i].set_xlabel("Lon in deg") for i in range(3)] +ax[0].set_ylabel("Lat in deg") + +ax[0].set_title("Temperature observations at 2m\nfrom DWD (2020-06-09 12:00)") +ax[1].set_title("Interpolated temperature\nwith North-South drift") +ax[2].set_title("Estimated mean drift\nfrom Universal Kriging") + +fmt = dict(orientation="horizontal", shrink=0.5, fraction=0.1, pad=0.2) +fig.colorbar(co2, ax=ax, **fmt).set_label("T in [°C]") + +############################################################################### +# To get a better impression of the estimated north-south drift, well take +# a look at a cross-section at a longitude of 10 degree: + +fig, ax = plt.subplots() +ax.plot(g_lat, field[:, 50], label="Interpolated temperature") +ax.plot(g_lat, mean[:, 50], label="North-South mean drift") +ax.set_xlabel("Lat in deg") +ax.set_ylabel("T in [°C]") +ax.set_title("North-South cross-section at 10°") +ax.legend() + +############################################################################### +# Interpretion of the results is now up to you! ;-) diff --git a/examples/08_geo_coordinates/de_borders.txt b/examples/08_geo_coordinates/de_borders.txt new file mode 100644 index 00000000..c8cdb5a8 --- /dev/null +++ b/examples/08_geo_coordinates/de_borders.txt @@ -0,0 +1,492 @@ +9.524023437500005684e+00 4.752421874999999574e+01 +9.350000000000022737e+00 4.759892578124999574e+01 +9.182812500000011369e+00 4.767070312499999574e+01 +9.127539062500005684e+00 4.767070312499999574e+01 +8.881152343750017053e+00 4.765639648437499432e+01 +8.874023437500000000e+00 4.766269531249999858e+01 +8.831152343750005684e+00 4.770361328125000000e+01 +8.793066406250005684e+00 4.771655273437500000e+01 +8.770117187500005684e+00 4.770991210937499716e+01 +8.754785156250022737e+00 4.769804687499999574e+01 +8.728320312500017053e+00 4.770004882812499858e+01 +8.617871093750011369e+00 4.776611328125000000e+01 +8.572656250000022737e+00 4.777563476562500000e+01 +8.509863281250005684e+00 4.776689453124999574e+01 +8.435742187500011369e+00 4.773134765624999432e+01 +8.403417968750005684e+00 4.768779296874999574e+01 +8.413281250000011369e+00 4.766269531249999858e+01 +8.451757812500005684e+00 4.765180664062499716e+01 +8.552343750000005684e+00 4.765913085937499716e+01 +8.567089843750011369e+00 4.765190429687499574e+01 +8.570507812500011369e+00 4.763779296874999858e+01 +8.559472656250022737e+00 4.762402343750000000e+01 +8.477636718750005684e+00 4.761269531249999432e+01 +8.454003906250022737e+00 4.759619140625000000e+01 +8.430078125000022737e+00 4.759213867187499858e+01 +8.414746093750011369e+00 4.758959960937500000e+01 +8.327832031250011369e+00 4.760693359375000000e+01 +8.198242187500000000e+00 4.760693359375000000e+01 +8.093750000000000000e+00 4.757617187500000000e+01 +7.927050781250017053e+00 4.756386718749999432e+01 +7.698046875000017053e+00 4.756987304687499574e+01 +7.615625000000022737e+00 4.759272460937499716e+01 +7.565429687500000000e+00 4.760654296874999858e+01 +7.529394531250005684e+00 4.767387695312499574e+01 +7.538574218750000000e+00 4.777363281249999716e+01 +7.593261718750000000e+00 4.790566406249999432e+01 +7.608496093750005684e+00 4.800258789062499432e+01 +7.584179687500011369e+00 4.806430664062499858e+01 +7.616601562500022737e+00 4.815678710937499574e+01 +7.705664062500005684e+00 4.828002929687500000e+01 +7.765136718750000000e+00 4.841000976562499858e+01 +7.794824218750022737e+00 4.854682617187499716e+01 +7.837988281250005684e+00 4.863603515624999574e+01 +7.922753906250022737e+00 4.869853515624999574e+01 +8.124023437500000000e+00 4.887329101562500000e+01 +8.140332031250011369e+00 4.888642578124999716e+01 +8.134863281250005684e+00 4.897358398437499716e+01 +8.080664062500005684e+00 4.898588867187499574e+01 +8.001269531250017053e+00 4.901093749999999716e+01 +7.799218750000022737e+00 4.904189453124999432e+01 +7.610937500000005684e+00 4.906176757812500000e+01 +7.525488281250005684e+00 4.908637695312499716e+01 +7.450585937500022737e+00 4.915219726562499858e+01 +7.404199218750022737e+00 4.915307617187500000e+01 +7.313378906250022737e+00 4.912954101562499432e+01 +7.199902343750011369e+00 4.911362304687499858e+01 +7.117382812500011369e+00 4.912753906249999858e+01 +7.065722656250017053e+00 4.912485351562499858e+01 +7.036718750000005684e+00 4.911269531249999432e+01 +7.022167968750011369e+00 4.912343749999999432e+01 +7.001464843750000000e+00 4.917988281249999716e+01 +6.958300781250017053e+00 4.919462890624999574e+01 +6.891210937500005684e+00 4.920751953125000000e+01 +6.849511718750022737e+00 4.920195312499999574e+01 +6.820703125000022737e+00 4.917392578124999858e+01 +6.776269531250022737e+00 4.915415039062499858e+01 +6.735449218750005684e+00 4.916059570312499716e+01 +6.607617187500011369e+00 4.929086914062499858e+01 +6.574707031250000000e+00 4.931967773437499858e+01 +6.566308593750022737e+00 4.934619140625000000e+01 +6.534277343750005684e+00 4.939467773437499432e+01 +6.458105468750005684e+00 4.944287109375000000e+01 +6.382226562500022737e+00 4.945815429687499432e+01 +6.344335937500005684e+00 4.945273437499999858e+01 +6.348437500000017053e+00 4.951269531250000000e+01 +6.378320312500022737e+00 4.959960937500000000e+01 +6.406738281250000000e+00 4.964497070312499716e+01 +6.444628906250017053e+00 4.968203124999999432e+01 +6.484765625000022737e+00 4.970781249999999574e+01 +6.493750000000005684e+00 4.975439453125000000e+01 +6.487304687500000000e+00 4.979848632812499432e+01 +6.440917968750000000e+00 4.980532226562499432e+01 +6.324609375000022737e+00 4.983789062500000000e+01 +6.256054687500011369e+00 4.987216796874999858e+01 +6.204882812500017053e+00 4.991513671874999858e+01 +6.138183593750000000e+00 4.997431640624999716e+01 +6.109765625000022737e+00 5.003437499999999716e+01 +6.108300781250022737e+00 5.009423828125000000e+01 +6.116503906250017053e+00 5.012099609374999432e+01 +6.121289062500011369e+00 5.013935546874999716e+01 +6.175097656250017053e+00 5.023266601562500000e+01 +6.364453125000011369e+00 5.031616210937500000e+01 +6.343652343750022737e+00 5.040024414062499858e+01 +6.340917968750005684e+00 5.045175781249999858e+01 +6.294921875000000000e+00 5.048549804687499432e+01 +6.203027343750022737e+00 5.049912109374999858e+01 +6.178710937500000000e+00 5.052250976562499574e+01 +6.168457031250000000e+00 5.054536132812499716e+01 +6.235937500000005684e+00 5.059667968750000000e+01 +6.154492187500011369e+00 5.063725585937499574e+01 +6.119433593750017053e+00 5.067924804687499574e+01 +6.005957031250005684e+00 5.073222656249999574e+01 +5.993945312500017053e+00 5.075043945312499716e+01 +6.048437500000005684e+00 5.090488281249999858e+01 +6.006835937500000000e+00 5.094995117187500000e+01 +5.955078125000000000e+00 5.097294921874999574e+01 +5.894726562500011369e+00 5.098422851562499858e+01 +5.867187500000000000e+00 5.100566406249999574e+01 +5.857519531250005684e+00 5.103012695312499858e+01 +5.868359375000011369e+00 5.104531249999999432e+01 +5.939257812500017053e+00 5.104082031249999574e+01 +5.961035156250005684e+00 5.105668945312499574e+01 +6.129980468750005684e+00 5.114741210937499716e+01 +6.136914062500011369e+00 5.116484374999999574e+01 +6.113378906250005684e+00 5.117470703124999432e+01 +6.082421875000022737e+00 5.117998046874999574e+01 +6.074804687500005684e+00 5.119902343749999574e+01 +6.075878906250011369e+00 5.122412109375000000e+01 +6.166210937500011369e+00 5.135483398437499858e+01 +6.192871093750000000e+00 5.141059570312499716e+01 +6.198828125000005684e+00 5.144999999999999574e+01 +6.193261718750022737e+00 5.148891601562499432e+01 +6.141601562500000000e+00 5.155009765624999574e+01 +6.091113281250017053e+00 5.159892578124999574e+01 +6.089355468750000000e+00 5.163779296874999858e+01 +6.052734375000000000e+00 5.165825195312499574e+01 +5.948535156250017053e+00 5.176240234374999716e+01 +5.948730468750000000e+00 5.180268554687499716e+01 +6.007617187500017053e+00 5.183398437500000000e+01 +6.089843750000000000e+00 5.185395507812499716e+01 +6.117187500000000000e+00 5.187041015624999574e+01 +6.166503906250000000e+00 5.188076171875000142e+01 +6.297070312500011369e+00 5.185073242187500142e+01 +6.355664062500011369e+00 5.182465820312499716e+01 +6.372167968750005684e+00 5.183002929687499716e+01 +6.425000000000011369e+00 5.185839843750000000e+01 +6.517578125000000000e+00 5.185395507812499716e+01 +6.741796875000005684e+00 5.191088867187500000e+01 +6.775195312500017053e+00 5.193828124999999574e+01 +6.800390625000005684e+00 5.196738281249999858e+01 +6.802441406250011369e+00 5.198017578124999716e+01 +6.715625000000017053e+00 5.203618164062499574e+01 +6.712988281250005684e+00 5.205688476562500000e+01 +6.724511718750022737e+00 5.208022460937500142e+01 +6.749023437500000000e+00 5.209868164062499574e+01 +6.800390625000005684e+00 5.211123046875000142e+01 +6.855078125000005684e+00 5.213579101562499574e+01 +6.977246093750011369e+00 5.220551757812499716e+01 +7.019628906250005684e+00 5.226601562500000142e+01 +7.032617187500022737e+00 5.233149414062499716e+01 +7.035156250000000000e+00 5.238022460937499858e+01 +7.001855468750022737e+00 5.241899414062499574e+01 +6.968164062500022737e+00 5.244409179687500000e+01 +6.922070312500011369e+00 5.244028320312499858e+01 +6.832519531250000000e+00 5.244228515625000142e+01 +6.748828125000017053e+00 5.246401367187500142e+01 +6.702929687500017053e+00 5.249921874999999716e+01 +6.691601562500011369e+00 5.253017578125000142e+01 +6.712402343750000000e+00 5.254965820312499858e+01 +6.718750000000000000e+00 5.257358398437499858e+01 +6.705371093750017053e+00 5.259765625000000000e+01 +6.710742187500017053e+00 5.261787109374999716e+01 +6.748437500000022737e+00 5.263408203124999574e+01 +7.013183593750000000e+00 5.263354492187500000e+01 +7.033007812500017053e+00 5.265136718750000000e+01 +7.050878906250005684e+00 5.274477539062500142e+01 +7.117089843750022737e+00 5.288701171874999574e+01 +7.179492187500017053e+00 5.296621093750000142e+01 +7.189941406250000000e+00 5.299951171875000000e+01 +7.188964843750000000e+00 5.318720703124999716e+01 +7.197265625000000000e+00 5.328227539062499574e+01 +7.152050781250011369e+00 5.332695312499999574e+01 +7.053320312500005684e+00 5.337583007812499858e+01 +7.074316406250005684e+00 5.347763671874999858e+01 +7.107128906250011369e+00 5.355698242187499858e+01 +7.206445312500022737e+00 5.365454101562500000e+01 +7.285253906250005684e+00 5.368134765624999716e+01 +7.629199218750017053e+00 5.369726562500000000e+01 +8.009277343750000000e+00 5.369072265624999574e+01 +8.167089843750005684e+00 5.354340820312499716e+01 +8.108496093750005684e+00 5.346767578125000142e+01 +8.200781250000005684e+00 5.343242187499999574e+01 +8.245214843750005684e+00 5.344531250000000000e+01 +8.279003906250011369e+00 5.351118164062499716e+01 +8.301562500000017053e+00 5.358413085937500142e+01 +8.333886718750022737e+00 5.360620117187500000e+01 +8.451367187500011369e+00 5.355170898437499716e+01 +8.492675781250000000e+00 5.351435546874999716e+01 +8.495214843750005684e+00 5.339423828124999716e+01 +8.538476562500022737e+00 5.355688476562500000e+01 +8.506250000000022737e+00 5.367075195312499858e+01 +8.528417968750005684e+00 5.378110351562499858e+01 +8.575585937500022737e+00 5.383847656249999858e+01 +8.618945312500017053e+00 5.387500000000000000e+01 +8.897753906250017053e+00 5.383569335937500000e+01 +9.205566406250000000e+00 5.385595703125000000e+01 +9.321972656250011369e+00 5.381347656250000000e+01 +9.585351562500022737e+00 5.360048828125000142e+01 +9.673144531250017053e+00 5.356562499999999716e+01 +9.783984375000017053e+00 5.355463867187499716e+01 +9.631250000000022737e+00 5.360019531249999858e+01 +9.312011718750000000e+00 5.385913085937500000e+01 +9.216406250000005684e+00 5.389121093749999858e+01 +9.069628906250017053e+00 5.390092773437499574e+01 +8.978125000000005684e+00 5.392622070312499716e+01 +8.920410156250000000e+00 5.396533203125000000e+01 +8.903515625000011369e+00 5.400029296874999574e+01 +8.906640625000022737e+00 5.426079101562499574e+01 +8.851562500000000000e+00 5.429956054687500000e+01 +8.780371093750005684e+00 5.431303710937499574e+01 +8.736035156250011369e+00 5.429521484374999574e+01 +8.644921875000022737e+00 5.429497070312499574e+01 +8.625781250000017053e+00 5.435395507812499716e+01 +8.648046875000005684e+00 5.439765624999999716e+01 +8.831152343750005684e+00 5.442753906249999574e+01 +8.951855468750011369e+00 5.446757812499999574e+01 +8.957226562500011369e+00 5.453833007812500000e+01 +8.880957031250005684e+00 5.459394531249999716e+01 +8.789648437500005684e+00 5.469594726562500142e+01 +8.682324218750011369e+00 5.479184570312499858e+01 +8.670312500000022737e+00 5.490341796874999858e+01 +8.670703125000017053e+00 5.490332031250000000e+01 +8.857226562500017053e+00 5.490112304687500000e+01 +8.902929687500005684e+00 5.489692382812499716e+01 +9.185839843750017053e+00 5.484467773437499716e+01 +9.254980468750005684e+00 5.480800781250000142e+01 +9.341992187500011369e+00 5.480629882812500142e+01 +9.498730468750011369e+00 5.484042968749999858e+01 +9.615820312500005684e+00 5.485541992187499716e+01 +9.661230468750005684e+00 5.483437500000000142e+01 +9.725000000000022737e+00 5.482553710937499858e+01 +9.739746093750000000e+00 5.482553710937499858e+01 +9.745898437500017053e+00 5.480717773437499574e+01 +9.892285156250011369e+00 5.478061523437499858e+01 +9.953808593750011369e+00 5.473828125000000000e+01 +1.002216796875001137e+01 5.467392578124999858e+01 +1.002880859375000000e+01 5.458129882812500000e+01 +9.941308593750022737e+00 5.451464843750000000e+01 +9.868652343750000000e+00 5.447246093749999574e+01 +1.014345703125002274e+01 5.448842773437500142e+01 +1.017080078125002274e+01 5.445019531250000000e+01 +1.021240234375000000e+01 5.440893554687500000e+01 +1.036044921875000568e+01 5.443833007812499858e+01 +1.073154296875000568e+01 5.431625976562499858e+01 +1.095595703125002274e+01 5.437568359374999716e+01 +1.101337890625001137e+01 5.437915039062500000e+01 +1.106435546875002274e+01 5.428051757812500000e+01 +1.100859375000001705e+01 5.418115234375000000e+01 +1.081074218750001137e+01 5.407514648437499716e+01 +1.085458984375000568e+01 5.400981445312499574e+01 +1.091777343750001705e+01 5.399531249999999716e+01 +1.110429687500001705e+01 5.400917968750000142e+01 +1.139960937500001137e+01 5.394462890624999574e+01 +1.146113281250001137e+01 5.396474609375000142e+01 +1.170058593750002274e+01 5.411352539062500000e+01 +1.179628906250002274e+01 5.414545898437499716e+01 +1.211132812500000000e+01 5.416831054687499858e+01 +1.216865234375001137e+01 5.422587890624999574e+01 +1.229628906250002274e+01 5.428378906249999858e+01 +1.237851562500000568e+01 5.434702148437499858e+01 +1.257539062500001137e+01 5.446738281249999858e+01 +1.277910156250001705e+01 5.444570312500000142e+01 +1.289804687500000568e+01 5.442265624999999574e+01 +1.302861328125001705e+01 5.441103515625000142e+01 +1.314746093750000000e+01 5.428271484375000000e+01 +1.344804687500001705e+01 5.414086914062500000e+01 +1.372421875000000568e+01 5.415322265625000142e+01 +1.382226562500000000e+01 5.401904296875000000e+01 +1.386552734375001705e+01 5.385336914062499858e+01 +1.395039062500001137e+01 5.380136718749999858e+01 +1.402500000000000568e+01 5.376743164062499858e+01 +1.425000000000000000e+01 5.373188476562499716e+01 +1.425888671875000568e+01 5.372963867187500142e+01 +1.426611328125000000e+01 5.370712890624999858e+01 +1.427988281250000568e+01 5.362475585937500000e+01 +1.429873046875002274e+01 5.355644531249999574e+01 +1.441455078125000000e+01 5.328349609374999574e+01 +1.441230468750001137e+01 5.321674804687499716e+01 +1.441093750000001705e+01 5.319902343749999574e+01 +1.436855468750002274e+01 5.310556640624999858e+01 +1.429316406250001137e+01 5.302675781250000142e+01 +1.419365234375001705e+01 5.298232421875000142e+01 +1.413886718750001137e+01 5.293286132812500000e+01 +1.412861328125001137e+01 5.287822265624999574e+01 +1.425371093750001705e+01 5.278251953124999574e+01 +1.451406250000002274e+01 5.264560546874999858e+01 +1.461943359375001705e+01 5.252851562499999716e+01 +1.456972656250002274e+01 5.243110351562499716e+01 +1.455458984375002274e+01 5.235966796874999574e+01 +1.457392578125001137e+01 5.231416015624999716e+01 +1.461562500000002274e+01 5.227763671874999574e+01 +1.467988281250001137e+01 5.225000000000000000e+01 +1.470537109375001705e+01 5.220747070312499716e+01 +1.469238281250000000e+01 5.215004882812500142e+01 +1.470458984375000000e+01 5.211020507812499858e+01 +1.475253906250000568e+01 5.208183593749999574e+01 +1.474814453125000568e+01 5.207080078125000000e+01 +1.472480468750001137e+01 5.203085937499999858e+01 +1.469296875000000568e+01 5.195800781250000000e+01 +1.467490234375000568e+01 5.190483398437499574e+01 +1.460166015625000568e+01 5.183237304687499858e+01 +1.462392578125002274e+01 5.177080078124999574e+01 +1.468134765625001137e+01 5.169819335937499716e+01 +1.472490234375001705e+01 5.166171874999999858e+01 +1.473867187500002274e+01 5.162714843749999716e+01 +1.471093750000000000e+01 5.154492187500000000e+01 +1.472470703125000568e+01 5.152387695312499716e+01 +1.490595703125001137e+01 5.146333007812499716e+01 +1.493554687500000000e+01 5.143535156249999574e+01 +1.495312500000000000e+01 5.137714843749999716e+01 +1.501660156250000000e+01 5.125273437499999574e+01 +1.496386718750000000e+01 5.109511718749999432e+01 +1.491748046875000000e+01 5.100874023437499716e+01 +1.481425781250001705e+01 5.087163085937499574e+01 +1.480937500000001705e+01 5.085898437499999858e+01 +1.479746093750000568e+01 5.084233398437499574e+01 +1.476650390625002274e+01 5.081831054687499716e+01 +1.472333984375001137e+01 5.081469726562500000e+01 +1.465820312500000000e+01 5.083261718749999858e+01 +1.461357421875001705e+01 5.085556640624999858e+01 +1.462382812500001705e+01 5.091474609374999716e+01 +1.459521484375000000e+01 5.091860351562499432e+01 +1.455966796875000568e+01 5.095493164062499858e+01 +1.454570312500001705e+01 5.099394531249999574e+01 +1.450732421875000000e+01 5.100986328124999858e+01 +1.436728515625000568e+01 5.102626953124999432e+01 +1.431972656250002274e+01 5.103779296874999716e+01 +1.428320312500000000e+01 5.102949218749999716e+01 +1.425585937500000000e+01 5.100185546874999432e+01 +1.427333984375002274e+01 5.097690429687499858e+01 +1.429941406250000568e+01 5.095258789062499716e+01 +1.437705078125000568e+01 5.091406250000000000e+01 +1.436904296875002274e+01 5.089873046874999574e+01 +1.420175781250000568e+01 5.086123046874999432e+01 +1.409648437500001705e+01 5.082275390625000000e+01 +1.399843750000002274e+01 5.080112304687499858e+01 +1.389853515625000568e+01 5.076127929687499574e+01 +1.370136718750001137e+01 5.071650390624999716e+01 +1.355673828125000568e+01 5.070463867187499574e+01 +1.352656250000001137e+01 5.069282226562499716e+01 +1.347255859375002274e+01 5.061694335937500000e+01 +1.343613281250000568e+01 5.060107421875000000e+01 +1.340117187500001705e+01 5.060932617187499716e+01 +1.337460937500000568e+01 5.062172851562499432e+01 +1.334101562500001137e+01 5.061142578124999858e+01 +1.330605468750002274e+01 5.058632812499999432e+01 +1.326953125000000000e+01 5.057641601562500000e+01 +1.323769531250002274e+01 5.057675781249999858e+01 +1.318115234375000000e+01 5.051049804687500000e+01 +1.301640625000001705e+01 5.049038085937499432e+01 +1.299707031250000000e+01 5.045605468750000000e+01 +1.296679687500000000e+01 5.041621093749999716e+01 +1.294267578125001705e+01 5.040644531249999716e+01 +1.286826171875000568e+01 5.042221679687499858e+01 +1.276542968750001705e+01 5.043095703124999574e+01 +1.270644531250002274e+01 5.040913085937499716e+01 +1.263554687500001705e+01 5.039707031249999858e+01 +1.254902343750001137e+01 5.039340820312499858e+01 +1.245263671875000000e+01 5.034980468749999716e+01 +1.235859375000001137e+01 5.027324218749999574e+01 +1.230566406250000000e+01 5.020571289062499432e+01 +1.227734375000000000e+01 5.018144531249999574e+01 +1.223115234375001137e+01 5.024487304687500000e+01 +1.217480468750000000e+01 5.028837890624999574e+01 +1.213486328125000568e+01 5.031093749999999432e+01 +1.209921875000000568e+01 5.031098632812499716e+01 +1.208984375000000000e+01 5.030175781250000000e+01 +1.208974609375002274e+01 5.026855468750000000e+01 +1.212783203125002274e+01 5.021342773437499574e+01 +1.217500000000001137e+01 5.017583007812499574e+01 +1.218251953125002274e+01 5.014804687499999858e+01 +1.220781250000001705e+01 5.009750976562499858e+01 +1.227646484375000568e+01 5.004233398437499858e+01 +1.238417968750002274e+01 4.999858398437499574e+01 +1.245761718750000568e+01 4.995551757812499716e+01 +1.251201171875001705e+01 4.989580078124999574e+01 +1.251250000000001705e+01 4.987744140625000000e+01 +1.249755859375000000e+01 4.985307617187499574e+01 +1.247187500000001137e+01 4.983007812500000000e+01 +1.245019531250000000e+01 4.980014648437499858e+01 +1.239052734375002274e+01 4.973964843749999432e+01 +1.240820312500000000e+01 4.971318359374999574e+01 +1.245703125000000000e+01 4.967978515624999858e+01 +1.250029296875001705e+01 4.963969726562499574e+01 +1.255576171875000568e+01 4.957485351562499432e+01 +1.263203125000001137e+01 4.946123046874999574e+01 +1.268115234375000000e+01 4.941450195312499716e+01 +1.274785156250001705e+01 4.936621093750000000e+01 +1.281337890625002274e+01 4.932934570312500000e+01 +1.291669921875001137e+01 4.933046874999999432e+01 +1.302373046875001705e+01 4.926010742187499858e+01 +1.314052734375002274e+01 4.915834960937499432e+01 +1.322783203125001705e+01 4.911166992187499858e+01 +1.328876953125001137e+01 4.909746093749999574e+01 +1.333906250000001137e+01 4.906079101562500000e+01 +1.338369140625002274e+01 4.900810546874999574e+01 +1.340117187500001705e+01 4.897758789062499574e+01 +1.344072265625001705e+01 4.895556640625000000e+01 +1.354765625000001705e+01 4.895966796874999716e+01 +1.368496093750002274e+01 4.887670898437500000e+01 +1.376992187500002274e+01 4.881596679687499574e+01 +1.381474609375001705e+01 4.876694335937499858e+01 +1.380292968750001137e+01 4.874750976562499716e+01 +1.379746093750000568e+01 4.868642578124999432e+01 +1.379882812500000000e+01 4.862167968749999858e+01 +1.378535156250001137e+01 4.858745117187499574e+01 +1.372392578125001705e+01 4.854238281249999432e+01 +1.369218750000001705e+01 4.853276367187499574e+01 +1.367519531250002274e+01 4.852304687499999858e+01 +1.348662109375001705e+01 4.858183593749999574e+01 +1.347167968750000000e+01 4.857182617187499574e+01 +1.345986328125002274e+01 4.856455078124999858e+01 +1.340937500000001137e+01 4.839414062499999858e+01 +1.337460937500000568e+01 4.836137695312499574e+01 +1.332285156250000568e+01 4.833124999999999716e+01 +1.321523437500002274e+01 4.830190429687499432e+01 +1.314042968750001705e+01 4.828994140624999432e+01 +1.308212890625000568e+01 4.827509765624999716e+01 +1.289746093750000000e+01 4.820371093749999858e+01 +1.281425781250001705e+01 4.816083984374999716e+01 +1.276035156250000568e+01 4.810698242187499574e+01 +1.276005859375001705e+01 4.807597656249999574e+01 +1.284990234375001705e+01 4.798481445312499716e+01 +1.295351562500002274e+01 4.789062500000000000e+01 +1.295419921875000568e+01 4.780776367187499432e+01 +1.290830078125000568e+01 4.774580078124999716e+01 +1.289765625000001137e+01 4.772187499999999716e+01 +1.292812500000002274e+01 4.771284179687499716e+01 +1.298554687500001137e+01 4.770942382812499716e+01 +1.303359375000002274e+01 4.769873046875000000e+01 +1.305410156250002274e+01 4.765512695312499858e+01 +1.304794921875000568e+01 4.757915039062499574e+01 +1.303154296875001705e+01 4.750800781249999716e+01 +1.301435546875001137e+01 4.747807617187499574e+01 +1.296806640625001705e+01 4.747568359374999858e+01 +1.287890625000000000e+01 4.750644531249999858e+01 +1.280937500000001705e+01 4.754218749999999716e+01 +1.278281250000000568e+01 4.756416015624999716e+01 +1.278115234375002274e+01 4.759042968749999858e+01 +1.279619140625001705e+01 4.760703124999999858e+01 +1.277138671875002274e+01 4.763940429687500000e+01 +1.268583984375001705e+01 4.766933593749999432e+01 +1.259423828125000000e+01 4.765629882812499574e+01 +1.252656250000001137e+01 4.763613281249999432e+01 +1.248291015625000000e+01 4.763730468749999858e+01 +1.243574218750001137e+01 4.766611328124999858e+01 +1.236318359375002274e+01 4.768818359374999716e+01 +1.226835937500001705e+01 4.770273437499999858e+01 +1.220927734375001705e+01 4.771826171875000000e+01 +1.219687500000000568e+01 4.770908203124999858e+01 +1.220380859375001137e+01 4.764672851562500000e+01 +1.218564453125000568e+01 4.761953124999999432e+01 +1.171679687500000000e+01 4.758349609375000000e+01 +1.157392578125001137e+01 4.754975585937499716e+01 +1.146992187500001137e+01 4.750610351562500000e+01 +1.139296875000002274e+01 4.748715820312499858e+01 +1.137412109375000568e+01 4.746025390624999574e+01 +1.129794921875000568e+01 4.742490234374999858e+01 +1.121191406250000000e+01 4.741362304687499574e+01 +1.119121093750001705e+01 4.742519531249999432e+01 +1.113603515625001705e+01 4.740888671874999716e+01 +1.104199218750000000e+01 4.739311523437499574e+01 +1.098085937500002274e+01 4.739814453124999716e+01 +1.095214843750000000e+01 4.742670898437499716e+01 +1.089394531250002274e+01 4.747045898437500000e+01 +1.087060546875000000e+01 4.750078124999999574e+01 +1.087304687500000000e+01 4.752021484374999716e+01 +1.074160156250002274e+01 4.752412109374999716e+01 +1.065869140625000000e+01 4.754721679687499858e+01 +1.048281250000002274e+01 4.754179687499999574e+01 +1.043945312500000000e+01 4.755156249999999574e+01 +1.043037109375001137e+01 4.754106445312499574e+01 +1.040390625000000568e+01 4.741699218750000000e+01 +1.036914062500000000e+01 4.736606445312499858e+01 +1.031279296875001705e+01 4.731342773437499716e+01 +1.024062500000002274e+01 4.728413085937499716e+01 +1.018300781250002274e+01 4.727880859375000000e+01 +1.018574218750001137e+01 4.731718749999999574e+01 +1.020029296875000568e+01 4.736342773437499432e+01 +1.015878906250000568e+01 4.737426757812500000e+01 +1.009648437500001705e+01 4.737958984374999716e+01 +1.006630859375002274e+01 4.739335937499999574e+01 +1.007421875000000000e+01 4.742851562499999574e+01 +1.005986328125001705e+01 4.744907226562499858e+01 +1.003408203125002274e+01 4.747358398437499716e+01 +9.971582031250022737e+00 4.750532226562499716e+01 +9.839160156250017053e+00 4.755229492187499574e+01 +9.748925781250022737e+00 4.757553710937499858e+01 +9.715136718750017053e+00 4.755078125000000000e+01 +9.650585937500011369e+00 4.752587890625000000e+01 +9.548925781250005684e+00 4.753403320312499858e+01 +9.524023437500005684e+00 4.752421874999999574e+01 diff --git a/examples/08_geo_coordinates/temp_obs.txt b/examples/08_geo_coordinates/temp_obs.txt new file mode 100644 index 00000000..aa8e60fc --- /dev/null +++ b/examples/08_geo_coordinates/temp_obs.txt @@ -0,0 +1,494 @@ +# id, lat, lon, temp +4.400000000000000000e+01 5.293359999999999843e+01 8.237000000000000099e+00 1.569999999999999929e+01 +7.300000000000000000e+01 4.861590000000000344e+01 1.305059999999999931e+01 1.390000000000000036e+01 +7.800000000000000000e+01 5.248530000000000229e+01 7.912600000000000300e+00 1.509999999999999964e+01 +9.100000000000000000e+01 5.074459999999999837e+01 9.345000000000000639e+00 1.700000000000000000e+01 +9.600000000000000000e+01 5.294369999999999976e+01 1.285180000000000078e+01 2.189999999999999858e+01 +1.020000000000000000e+02 5.386330000000000240e+01 8.127499999999999503e+00 1.190000000000000036e+01 +1.250000000000000000e+02 4.783420000000000272e+01 1.086669999999999980e+01 1.140000000000000036e+01 +1.310000000000000000e+02 5.108809999999999718e+01 1.293260000000000076e+01 1.719999999999999929e+01 +1.420000000000000000e+02 4.840599999999999881e+01 1.131170000000000009e+01 1.290000000000000036e+01 +1.500000000000000000e+02 4.972729999999999961e+01 8.116400000000000503e+00 1.719999999999999929e+01 +1.510000000000000000e+02 4.946909999999999741e+01 1.185459999999999958e+01 1.340000000000000036e+01 +1.540000000000000000e+02 4.801970000000000027e+01 1.229250000000000043e+01 1.390000000000000036e+01 +1.610000000000000000e+02 5.042369999999999663e+01 7.420200000000000351e+00 1.810000000000000142e+01 +1.640000000000000000e+02 5.303159999999999741e+01 1.399080000000000013e+01 2.130000000000000071e+01 +1.670000000000000000e+02 5.384120000000000061e+01 1.368459999999999965e+01 2.130000000000000071e+01 +1.830000000000000000e+02 5.467920000000000158e+01 1.343430000000000035e+01 1.739999999999999858e+01 +1.910000000000000000e+02 4.996940000000000026e+01 9.911400000000000432e+00 1.860000000000000142e+01 +1.980000000000000000e+02 5.137449999999999761e+01 1.129199999999999982e+01 2.019999999999999929e+01 +2.170000000000000000e+02 4.787740000000000151e+01 1.136430000000000007e+01 1.269999999999999929e+01 +2.220000000000000000e+02 5.059080000000000155e+01 1.271390000000000065e+01 1.580000000000000071e+01 +2.320000000000000000e+02 4.842530000000000001e+01 1.094170000000000087e+01 1.340000000000000036e+01 +2.570000000000000000e+02 4.872699999999999676e+01 8.245699999999999363e+00 1.359999999999999964e+01 +2.590000000000000000e+02 4.780639999999999645e+01 7.638700000000000045e+00 1.440000000000000036e+01 +2.820000000000000000e+02 4.987429999999999808e+01 1.092060000000000031e+01 1.580000000000000071e+01 +2.940000000000000000e+02 5.231989999999999696e+01 9.429999999999999716e+00 2.150000000000000000e+01 +2.980000000000000000e+02 5.434060000000000201e+01 1.271080000000000076e+01 1.769999999999999929e+01 +3.030000000000000000e+02 5.206139999999999901e+01 1.349959999999999916e+01 2.119999999999999929e+01 +3.140000000000000000e+02 5.116040000000000276e+01 1.450420000000000087e+01 2.039999999999999858e+01 +3.200000000000000000e+02 4.996670000000000300e+01 1.151970000000000027e+01 1.469999999999999929e+01 +3.300000000000000000e+02 4.956170000000000186e+01 8.967299999999999827e+00 1.409999999999999964e+01 +3.420000000000000000e+02 5.231700000000000017e+01 8.169399999999999551e+00 1.639999999999999858e+01 +3.680000000000000000e+02 5.281519999999999726e+01 9.924799999999999400e+00 1.919999999999999929e+01 +3.770000000000000000e+02 4.910699999999999932e+01 7.996699999999999697e+00 1.669999999999999929e+01 +3.790000000000000000e+02 5.090740000000000265e+01 1.126650000000000063e+01 1.810000000000000142e+01 +3.900000000000000000e+02 5.098369999999999891e+01 8.368299999999999628e+00 1.480000000000000071e+01 +4.000000000000000000e+02 5.263089999999999691e+01 1.350220000000000020e+01 2.250000000000000000e+01 +4.030000000000000000e+02 5.245369999999999777e+01 1.330170000000000030e+01 1.919999999999999929e+01 +4.100000000000000000e+02 5.240400000000000347e+01 1.373090000000000011e+01 2.080000000000000071e+01 +4.200000000000000000e+02 5.254469999999999885e+01 1.355979999999999919e+01 2.169999999999999929e+01 +4.270000000000000000e+02 5.238069999999999737e+01 1.353059999999999974e+01 2.250000000000000000e+01 +4.300000000000000000e+02 5.256439999999999912e+01 1.330879999999999974e+01 2.060000000000000142e+01 +4.330000000000000000e+02 5.246750000000000114e+01 1.340210000000000079e+01 2.060000000000000142e+01 +4.450000000000000000e+02 5.182180000000000319e+01 1.171100000000000030e+01 2.200000000000000000e+01 +4.600000000000000000e+02 4.926409999999999911e+01 6.686799999999999855e+00 1.519999999999999929e+01 +5.350000000000000000e+02 5.003719999999999857e+01 7.307900000000000063e+00 1.769999999999999929e+01 +5.910000000000000000e+02 5.339110000000000156e+01 1.068779999999999930e+01 1.889999999999999858e+01 +5.960000000000000000e+02 5.400280000000000058e+01 1.119079999999999941e+01 1.630000000000000071e+01 +6.030000000000000000e+02 5.072930000000000206e+01 7.203999999999999737e+00 1.760000000000000142e+01 +6.170000000000000000e+02 5.187299999999999756e+01 6.886300000000000310e+00 1.559999999999999964e+01 +6.560000000000000000e+02 5.172339999999999804e+01 1.060210000000000008e+01 1.580000000000000071e+01 +6.620000000000000000e+02 5.229149999999999920e+01 1.044640000000000057e+01 1.939999999999999858e+01 +6.910000000000000000e+02 5.304500000000000171e+01 8.797900000000000276e+00 1.650000000000000000e+01 +7.010000000000000000e+02 5.353320000000000078e+01 8.576100000000000279e+00 1.380000000000000071e+01 +7.040000000000000000e+02 5.344509999999999650e+01 9.138999999999999346e+00 1.610000000000000142e+01 +7.220000000000000000e+02 5.179860000000000042e+01 1.061829999999999963e+01 1.009999999999999964e+01 +7.550000000000000000e+02 4.951820000000000022e+01 9.321300000000000807e+00 1.600000000000000000e+01 +7.570000000000000000e+02 4.796249999999999858e+01 7.998300000000000409e+00 1.280000000000000071e+01 +7.600000000000000000e+02 5.336290000000000333e+01 9.943500000000000227e+00 1.750000000000000000e+01 +7.660000000000000000e+02 5.017459999999999809e+01 7.059499999999999886e+00 1.680000000000000071e+01 +7.690000000000000000e+02 5.228170000000000073e+01 9.088900000000000645e+00 1.889999999999999858e+01 +8.170000000000000000e+02 5.103059999999999974e+01 8.814600000000000435e+00 1.850000000000000000e+01 +8.400000000000000000e+02 5.043130000000000024e+01 1.261139999999999972e+01 1.080000000000000071e+01 +8.500000000000000000e+02 5.259590000000000032e+01 1.002960000000000029e+01 1.939999999999999858e+01 +8.530000000000000000e+02 5.079129999999999967e+01 1.287199999999999989e+01 1.710000000000000142e+01 +8.560000000000000000e+02 4.788430000000000319e+01 1.254039999999999999e+01 1.330000000000000071e+01 +8.670000000000000000e+02 5.030660000000000309e+01 1.096790000000000020e+01 1.739999999999999858e+01 +8.800000000000000000e+02 5.177600000000000335e+01 1.431680000000000064e+01 1.989999999999999858e+01 +8.910000000000000000e+02 5.387129999999999797e+01 8.705799999999999983e+00 1.469999999999999929e+01 +8.960000000000000000e+02 5.107780000000000342e+01 1.086190000000000033e+01 1.869999999999999929e+01 +9.170000000000000000e+02 4.988089999999999691e+01 8.677899999999999281e+00 1.639999999999999858e+01 +9.530000000000000000e+02 4.976189999999999714e+01 7.054199999999999804e+00 1.660000000000000142e+01 +9.540000000000000000e+02 5.417960000000000065e+01 7.458700000000000330e+00 1.250000000000000000e+01 +9.630000000000000000e+02 5.258809999999999718e+01 8.342399999999999594e+00 1.660000000000000142e+01 +9.790000000000000000e+02 5.073640000000000327e+01 8.267200000000000770e+00 2.060000000000000142e+01 +9.830000000000000000e+02 4.855619999999999692e+01 1.055990000000000073e+01 1.350000000000000000e+01 +9.910000000000000000e+02 5.091159999999999997e+01 1.370870000000000033e+01 1.869999999999999929e+01 +1.001000000000000000e+03 5.164509999999999934e+01 1.357469999999999999e+01 2.110000000000000142e+01 +1.048000000000000000e+03 5.112800000000000011e+01 1.375430000000000064e+01 1.869999999999999929e+01 +1.050000000000000000e+03 5.102210000000000178e+01 1.384699999999999953e+01 2.000000000000000000e+01 +1.051000000000000000e+03 5.102479999999999905e+01 1.377500000000000036e+01 1.980000000000000071e+01 +1.052000000000000000e+03 5.221739999999999782e+01 1.216409999999999947e+01 2.150000000000000000e+01 +1.072000000000000000e+03 4.947189999999999799e+01 8.192899999999999849e+00 1.669999999999999929e+01 +1.078000000000000000e+03 5.129599999999999937e+01 6.768600000000000172e+00 1.660000000000000142e+01 +1.103000000000000000e+03 4.810029999999999717e+01 1.198719999999999963e+01 1.340000000000000036e+01 +1.107000000000000000e+03 4.985199999999999676e+01 1.049910000000000032e+01 1.559999999999999964e+01 +1.161000000000000000e+03 4.887769999999999726e+01 1.123489999999999966e+01 1.450000000000000000e+01 +1.197000000000000000e+03 4.898949999999999960e+01 1.013119999999999976e+01 1.459999999999999964e+01 +1.200000000000000000e+03 5.406909999999999883e+01 9.010500000000000398e+00 1.519999999999999929e+01 +1.207000000000000000e+03 5.027049999999999841e+01 1.227420000000000044e+01 1.450000000000000000e+01 +1.214000000000000000e+03 4.820120000000000005e+01 8.108800000000000452e+00 1.250000000000000000e+01 +1.224000000000000000e+03 4.813779999999999859e+01 7.835099999999999731e+00 1.400000000000000000e+01 +1.228000000000000000e+03 5.416510000000000247e+01 6.346000000000000085e+00 1.180000000000000071e+01 +1.246000000000000000e+03 5.184179999999999922e+01 8.060700000000000642e+00 1.819999999999999929e+01 +1.262000000000000000e+03 4.834770000000000323e+01 1.181339999999999968e+01 1.350000000000000000e+01 +1.266000000000000000e+03 5.429919999999999902e+01 9.316200000000000259e+00 1.559999999999999964e+01 +1.270000000000000000e+03 5.098290000000000077e+01 1.096080000000000076e+01 1.639999999999999858e+01 +1.279000000000000000e+03 4.964970000000000283e+01 1.100740000000000052e+01 1.480000000000000071e+01 +1.297000000000000000e+03 5.120409999999999684e+01 1.001379999999999981e+01 1.639999999999999858e+01 +1.300000000000000000e+03 5.125399999999999778e+01 8.156499999999999417e+00 1.309999999999999964e+01 +1.303000000000000000e+03 5.140409999999999968e+01 6.967699999999999783e+00 1.559999999999999964e+01 +1.327000000000000000e+03 5.071189999999999998e+01 6.790499999999999758e+00 1.440000000000000036e+01 +1.332000000000000000e+03 4.848319999999999652e+01 1.272409999999999997e+01 1.290000000000000036e+01 +1.339000000000000000e+03 5.291570000000000107e+01 1.018849999999999945e+01 1.869999999999999929e+01 +1.346000000000000000e+03 4.787489999999999668e+01 8.003800000000000026e+00 4.700000000000000178e+00 +1.357000000000000000e+03 4.998069999999999879e+01 1.183760000000000012e+01 1.340000000000000036e+01 +1.358000000000000000e+03 5.042830000000000013e+01 1.295350000000000001e+01 9.400000000000000355e+00 +1.411000000000000000e+03 5.053090000000000259e+01 1.004800000000000004e+01 1.300000000000000000e+01 +1.420000000000000000e+03 5.002590000000000003e+01 8.521300000000000097e+00 1.800000000000000000e+01 +1.424000000000000000e+03 5.012689999999999912e+01 8.669399999999999551e+00 1.819999999999999929e+01 +1.443000000000000000e+03 4.802320000000000277e+01 7.834299999999999820e+00 1.380000000000000071e+01 +1.451000000000000000e+03 5.382770000000000010e+01 9.249299999999999855e+00 1.619999999999999929e+01 +1.468000000000000000e+03 4.845380000000000109e+01 8.409000000000000696e+00 9.599999999999999645e+00 +1.503000000000000000e+03 5.306430000000000291e+01 7.902199999999999669e+00 1.559999999999999964e+01 +1.504000000000000000e+03 5.111899999999999977e+01 9.279899999999999594e+00 1.869999999999999929e+01 +1.526000000000000000e+03 5.056680000000000064e+01 9.653299999999999770e+00 1.819999999999999929e+01 +1.544000000000000000e+03 5.251290000000000191e+01 1.139409999999999989e+01 2.100000000000000000e+01 +1.550000000000000000e+03 4.748299999999999699e+01 1.106209999999999916e+01 1.440000000000000036e+01 +1.580000000000000000e+03 4.998590000000000089e+01 7.954799999999999649e+00 1.989999999999999858e+01 +1.584000000000000000e+03 4.792419999999999902e+01 8.647299999999999542e+00 1.130000000000000071e+01 +1.587000000000000000e+03 4.894809999999999661e+01 1.142890000000000050e+01 1.309999999999999964e+01 +1.590000000000000000e+03 5.149419999999999931e+01 6.246299999999999741e+00 1.680000000000000071e+01 +1.602000000000000000e+03 4.843299999999999983e+01 7.993000000000000327e+00 1.450000000000000000e+01 +1.605000000000000000e+03 5.238750000000000284e+01 1.216009999999999991e+01 2.230000000000000071e+01 +1.612000000000000000e+03 5.088130000000000308e+01 1.212889999999999979e+01 1.830000000000000071e+01 +1.639000000000000000e+03 5.060170000000000101e+01 8.643900000000000361e+00 1.760000000000000142e+01 +1.645000000000000000e+03 5.096560000000000201e+01 9.050000000000000711e+00 1.930000000000000071e+01 +1.666000000000000000e+03 5.482730000000000103e+01 9.505800000000000693e+00 1.469999999999999929e+01 +1.684000000000000000e+03 5.116219999999999857e+01 1.495059999999999967e+01 1.869999999999999929e+01 +1.691000000000000000e+03 5.150019999999999953e+01 9.950699999999999434e+00 1.580000000000000071e+01 +1.694000000000000000e+03 5.360600000000000165e+01 1.210330000000000084e+01 1.960000000000000142e+01 +1.721000000000000000e+03 4.966400000000000148e+01 1.122390000000000043e+01 1.259999999999999964e+01 +1.735000000000000000e+03 4.878940000000000055e+01 1.362899999999999956e+01 1.280000000000000071e+01 +1.736000000000000000e+03 5.357309999999999661e+01 1.067970000000000041e+01 1.800000000000000000e+01 +1.757000000000000000e+03 5.409669999999999845e+01 1.340559999999999974e+01 1.919999999999999929e+01 +1.759000000000000000e+03 5.424369999999999692e+01 1.391019999999999968e+01 1.850000000000000000e+01 +1.766000000000000000e+03 5.213439999999999941e+01 7.696900000000000297e+00 1.610000000000000142e+01 +1.832000000000000000e+03 4.911290000000000333e+01 1.313380000000000081e+01 7.000000000000000000e+00 +1.863000000000000000e+03 5.026670000000000016e+01 9.185399999999999565e+00 1.730000000000000071e+01 +1.869000000000000000e+03 5.331530000000000058e+01 1.393379999999999974e+01 2.019999999999999929e+01 +1.886000000000000000e+03 4.848780000000000001e+01 1.026079999999999970e+01 1.290000000000000036e+01 +1.964000000000000000e+03 4.994449999999999790e+01 6.382100000000000328e+00 1.730000000000000071e+01 +1.975000000000000000e+03 5.363320000000000221e+01 9.988099999999999312e+00 1.739999999999999858e+01 +1.981000000000000000e+03 5.347769999999999868e+01 9.895699999999999719e+00 1.860000000000000142e+01 +2.014000000000000000e+03 5.246439999999999770e+01 9.677899999999999281e+00 2.039999999999999858e+01 +2.023000000000000000e+03 4.879180000000000206e+01 1.070620000000000083e+01 1.390000000000000036e+01 +2.039000000000000000e+03 5.190019999999999811e+01 1.056990000000000052e+01 1.789999999999999858e+01 +2.044000000000000000e+03 5.165200000000000102e+01 1.113669999999999938e+01 1.750000000000000000e+01 +2.074000000000000000e+03 4.837519999999999953e+01 8.980000000000000426e+00 1.130000000000000071e+01 +2.110000000000000000e+03 5.104110000000000014e+01 6.104199999999999626e+00 1.490000000000000036e+01 +2.115000000000000000e+03 5.417499999999999716e+01 7.892000000000000348e+00 1.400000000000000000e+01 +2.171000000000000000e+03 5.085199999999999676e+01 9.737700000000000244e+00 1.869999999999999929e+01 +2.174000000000000000e+03 5.162550000000000239e+01 1.036950000000000038e+01 1.619999999999999929e+01 +2.201000000000000000e+03 5.457500000000000284e+01 1.310440000000000005e+01 1.680000000000000071e+01 +2.211000000000000000e+03 5.073709999999999809e+01 7.652800000000000047e+00 1.430000000000000071e+01 +2.252000000000000000e+03 5.089900000000000091e+01 1.474569999999999936e+01 1.880000000000000071e+01 +2.261000000000000000e+03 5.031230000000000047e+01 1.187599999999999945e+01 1.590000000000000036e+01 +2.290000000000000000e+03 4.780089999999999861e+01 1.101079999999999970e+01 9.800000000000000711e+00 +2.303000000000000000e+03 5.431459999999999866e+01 9.538999999999999702e+00 1.630000000000000071e+01 +2.306000000000000000e+03 5.431940000000000168e+01 1.067319999999999958e+01 1.480000000000000071e+01 +2.315000000000000000e+03 5.176570000000000249e+01 1.316660000000000075e+01 1.930000000000000071e+01 +2.319000000000000000e+03 4.788230000000000075e+01 1.169609999999999950e+01 1.309999999999999964e+01 +2.323000000000000000e+03 5.185289999999999822e+01 9.495300000000000296e+00 1.960000000000000142e+01 +2.362000000000000000e+03 5.056510000000000105e+01 7.484300000000000175e+00 1.650000000000000000e+01 +2.385000000000000000e+03 4.969270000000000209e+01 7.326399999999999579e+00 1.630000000000000071e+01 +2.410000000000000000e+03 4.871119999999999806e+01 1.153619999999999912e+01 1.350000000000000000e+01 +2.429000000000000000e+03 5.398969999999999914e+01 9.569599999999999440e+00 1.689999999999999858e+01 +2.437000000000000000e+03 5.445700000000000074e+01 9.520300000000000651e+00 1.580000000000000071e+01 +2.444000000000000000e+03 5.092510000000000048e+01 1.158300000000000018e+01 2.019999999999999929e+01 +2.480000000000000000e+03 5.006430000000000291e+01 8.993000000000000327e+00 1.810000000000000142e+01 +2.483000000000000000e+03 5.118030000000000257e+01 8.489100000000000534e+00 1.380000000000000071e+01 +2.485000000000000000e+03 4.891700000000000159e+01 9.687099999999999156e+00 1.380000000000000071e+01 +2.486000000000000000e+03 4.942620000000000147e+01 7.755700000000000038e+00 1.669999999999999929e+01 +2.497000000000000000e+03 5.050139999999999674e+01 6.526399999999999757e+00 1.269999999999999929e+01 +2.559000000000000000e+03 4.772330000000000183e+01 1.033479999999999954e+01 1.240000000000000036e+01 +2.564000000000000000e+03 5.437760000000000105e+01 1.014240000000000030e+01 1.530000000000000071e+01 +2.575000000000000000e+03 4.918039999999999878e+01 9.980000000000000426e+00 1.469999999999999929e+01 +2.578000000000000000e+03 5.399949999999999761e+01 1.143410000000000082e+01 1.639999999999999858e+01 +2.597000000000000000e+03 5.022399999999999665e+01 1.007920000000000016e+01 1.730000000000000071e+01 +2.600000000000000000e+03 4.973629999999999995e+01 1.017810000000000059e+01 1.769999999999999929e+01 +2.601000000000000000e+03 5.022180000000000177e+01 8.446899999999999409e+00 1.259999999999999964e+01 +2.618000000000000000e+03 5.084579999999999700e+01 1.048029999999999973e+01 1.409999999999999964e+01 +2.627000000000000000e+03 5.155539999999999878e+01 1.388449999999999918e+01 2.030000000000000071e+01 +2.629000000000000000e+03 5.176120000000000232e+01 6.095399999999999707e+00 1.580000000000000071e+01 +2.638000000000000000e+03 4.810540000000000305e+01 8.754799999999999471e+00 9.000000000000000000e+00 +2.641000000000000000e+03 5.151850000000000307e+01 1.290649999999999942e+01 2.000000000000000000e+01 +2.667000000000000000e+03 5.086460000000000292e+01 7.157499999999999751e+00 1.830000000000000071e+01 +2.680000000000000000e+03 5.028399999999999892e+01 1.044560000000000066e+01 1.800000000000000000e+01 +2.700000000000000000e+03 4.883019999999999783e+01 1.148719999999999963e+01 1.380000000000000071e+01 +2.704000000000000000e+03 5.175110000000000099e+01 1.200939999999999941e+01 2.069999999999999929e+01 +2.708000000000000000e+03 4.766519999999999868e+01 1.108050000000000068e+01 1.269999999999999929e+01 +2.712000000000000000e+03 4.769519999999999982e+01 9.130699999999999150e+00 1.450000000000000000e+01 +2.750000000000000000e+03 5.025229999999999819e+01 1.132089999999999996e+01 1.719999999999999929e+01 +2.773000000000000000e+03 4.942830000000000013e+01 1.190160000000000018e+01 1.269999999999999929e+01 +2.794000000000000000e+03 5.293630000000000280e+01 1.240930000000000000e+01 2.150000000000000000e+01 +2.796000000000000000e+03 5.391559999999999775e+01 1.227899999999999991e+01 1.880000000000000071e+01 +2.812000000000000000e+03 4.836469999999999914e+01 7.828000000000000291e+00 1.450000000000000000e+01 +2.814000000000000000e+03 4.851209999999999667e+01 9.764499999999999957e+00 1.119999999999999929e+01 +2.856000000000000000e+03 5.191729999999999734e+01 1.308779999999999966e+01 1.900000000000000000e+01 +2.878000000000000000e+03 5.139090000000000202e+01 1.187860000000000049e+01 1.900000000000000000e+01 +2.886000000000000000e+03 4.821759999999999735e+01 9.909700000000000841e+00 1.200000000000000000e+01 +2.905000000000000000e+03 4.818489999999999895e+01 1.085069999999999979e+01 1.230000000000000071e+01 +2.907000000000000000e+03 5.479030000000000200e+01 8.951399999999999579e+00 1.340000000000000036e+01 +2.925000000000000000e+03 5.139330000000000354e+01 1.031230000000000047e+01 1.739999999999999858e+01 +2.928000000000000000e+03 5.131510000000000105e+01 1.244619999999999926e+01 2.110000000000000142e+01 +2.932000000000000000e+03 5.143480000000000274e+01 1.223959999999999937e+01 2.030000000000000071e+01 +2.947000000000000000e+03 5.113329999999999842e+01 8.034800000000000608e+00 1.540000000000000036e+01 +2.951000000000000000e+03 5.310070000000000334e+01 1.148639999999999972e+01 1.989999999999999858e+01 +2.953000000000000000e+03 4.785969999999999658e+01 8.230800000000000338e+00 9.599999999999999645e+00 +2.961000000000000000e+03 5.449960000000000093e+01 1.027369999999999983e+01 1.380000000000000071e+01 +2.968000000000000000e+03 5.098940000000000339e+01 6.977699999999999569e+00 1.730000000000000071e+01 +2.985000000000000000e+03 5.093829999999999814e+01 1.420930000000000071e+01 1.939999999999999858e+01 +3.015000000000000000e+03 5.220850000000000080e+01 1.411800000000000033e+01 2.010000000000000142e+01 +3.028000000000000000e+03 5.178540000000000276e+01 8.838800000000000878e+00 1.780000000000000071e+01 +3.031000000000000000e+03 5.163360000000000127e+01 8.394500000000000739e+00 1.900000000000000000e+01 +3.032000000000000000e+03 5.501100000000000279e+01 8.412499999999999645e+00 1.340000000000000036e+01 +3.034000000000000000e+03 5.045049999999999812e+01 1.163499999999999979e+01 1.639999999999999858e+01 +3.042000000000000000e+03 5.056170000000000186e+01 8.238599999999999923e+00 1.930000000000000071e+01 +3.083000000000000000e+03 5.192669999999999675e+01 1.387969999999999970e+01 2.119999999999999929e+01 +3.086000000000000000e+03 5.380250000000000199e+01 1.069890000000000008e+01 1.789999999999999858e+01 +3.093000000000000000e+03 5.297240000000000038e+01 1.113739999999999952e+01 1.860000000000000142e+01 +3.098000000000000000e+03 5.124519999999999698e+01 7.642500000000000071e+00 1.530000000000000071e+01 +3.126000000000000000e+03 5.210289999999999822e+01 1.158270000000000088e+01 2.050000000000000000e+01 +3.137000000000000000e+03 4.996560000000000201e+01 8.213900000000000645e+00 1.689999999999999858e+01 +3.147000000000000000e+03 4.877250000000000085e+01 1.221790000000000020e+01 1.330000000000000071e+01 +3.155000000000000000e+03 5.010150000000000148e+01 6.800900000000000389e+00 1.610000000000000142e+01 +3.158000000000000000e+03 5.254679999999999751e+01 1.454519999999999946e+01 2.089999999999999858e+01 +3.164000000000000000e+03 5.084920000000000329e+01 8.774599999999999511e+00 1.960000000000000142e+01 +3.166000000000000000e+03 5.065100000000000335e+01 1.314690000000000047e+01 1.450000000000000000e+01 +3.167000000000000000e+03 5.066210000000000235e+01 7.960300000000000153e+00 1.559999999999999964e+01 +3.196000000000000000e+03 5.332229999999999848e+01 1.193190000000000062e+01 2.180000000000000071e+01 +3.204000000000000000e+03 5.073349999999999937e+01 1.088150000000000084e+01 1.660000000000000142e+01 +3.226000000000000000e+03 5.172590000000000288e+01 1.151089999999999947e+01 2.050000000000000000e+01 +3.231000000000000000e+03 5.056119999999999948e+01 1.037710000000000043e+01 1.540000000000000036e+01 +3.234000000000000000e+03 5.112939999999999685e+01 1.343280000000000030e+01 1.880000000000000071e+01 +3.244000000000000000e+03 4.798199999999999932e+01 1.013840000000000074e+01 1.200000000000000000e+01 +3.254000000000000000e+03 5.271560000000000201e+01 7.317599999999999660e+00 1.610000000000000142e+01 +3.257000000000000000e+03 4.947729999999999961e+01 9.762199999999999989e+00 1.669999999999999929e+01 +3.268000000000000000e+03 4.816940000000000310e+01 8.943300000000000693e+00 9.500000000000000000e+00 +3.271000000000000000e+03 4.885479999999999734e+01 1.291890000000000072e+01 1.469999999999999929e+01 +3.278000000000000000e+03 4.853770000000000095e+01 9.273400000000000531e+00 1.240000000000000036e+01 +3.284000000000000000e+03 4.966910000000000025e+01 9.008499999999999730e+00 1.630000000000000071e+01 +3.287000000000000000e+03 4.971759999999999735e+01 9.099700000000000344e+00 1.490000000000000036e+01 +3.289000000000000000e+03 5.072809999999999775e+01 1.178379999999999939e+01 1.789999999999999858e+01 +3.307000000000000000e+03 4.747789999999999822e+01 1.126529999999999987e+01 1.250000000000000000e+01 +3.319000000000000000e+03 4.976440000000000197e+01 9.253000000000000114e+00 1.630000000000000071e+01 +3.340000000000000000e+03 5.043829999999999814e+01 7.806099999999999817e+00 1.789999999999999858e+01 +3.362000000000000000e+03 4.897209999999999752e+01 8.873400000000000176e+00 1.390000000000000036e+01 +3.366000000000000000e+03 4.827900000000000347e+01 1.250239999999999974e+01 1.400000000000000000e+01 +3.376000000000000000e+03 5.251760000000000161e+01 1.412320000000000064e+01 2.069999999999999929e+01 +3.379000000000000000e+03 4.816320000000000334e+01 1.154289999999999949e+01 1.390000000000000036e+01 +3.402000000000000000e+03 4.838510000000000133e+01 9.483700000000000685e+00 1.059999999999999964e+01 +3.426000000000000000e+03 5.156600000000000250e+01 1.470079999999999920e+01 1.989999999999999858e+01 +3.442000000000000000e+03 5.035739999999999839e+01 8.750600000000000378e+00 1.610000000000000142e+01 +3.484000000000000000e+03 4.870859999999999701e+01 1.121470000000000056e+01 1.580000000000000071e+01 +3.485000000000000000e+03 4.831150000000000233e+01 1.037729999999999997e+01 1.340000000000000036e+01 +3.490000000000000000e+03 5.053459999999999752e+01 7.085300000000000153e+00 1.730000000000000071e+01 +3.509000000000000000e+03 5.310199999999999676e+01 1.304209999999999958e+01 2.069999999999999929e+01 +3.513000000000000000e+03 5.050019999999999953e+01 1.113439999999999941e+01 1.280000000000000071e+01 +3.527000000000000000e+03 5.089229999999999876e+01 9.404999999999999361e+00 1.590000000000000036e+01 +3.540000000000000000e+03 5.084459999999999980e+01 7.371999999999999886e+00 1.789999999999999858e+01 +3.545000000000000000e+03 4.934400000000000119e+01 7.229700000000000237e+00 1.839999999999999858e+01 +3.571000000000000000e+03 4.981739999999999924e+01 1.186379999999999946e+01 1.369999999999999929e+01 +3.591000000000000000e+03 5.067430000000000234e+01 6.424000000000000377e+00 1.240000000000000036e+01 +3.603000000000000000e+03 4.938949999999999818e+01 9.966699999999999449e+00 1.459999999999999964e+01 +3.612000000000000000e+03 5.267110000000000269e+01 9.222899999999999210e+00 1.919999999999999929e+01 +3.621000000000000000e+03 4.882529999999999859e+01 1.050670000000000037e+01 1.390000000000000036e+01 +3.623000000000000000e+03 5.082939999999999969e+01 6.660199999999999676e+00 1.459999999999999964e+01 +3.631000000000000000e+03 5.371229999999999905e+01 7.151900000000000368e+00 1.380000000000000071e+01 +3.639000000000000000e+03 5.376469999999999771e+01 8.658300000000000551e+00 1.480000000000000071e+01 +3.660000000000000000e+03 5.036019999999999897e+01 6.869699999999999918e+00 1.469999999999999929e+01 +3.667000000000000000e+03 4.942580000000000240e+01 1.125380000000000003e+01 1.340000000000000036e+01 +3.668000000000000000e+03 4.950300000000000011e+01 1.105489999999999995e+01 1.419999999999999929e+01 +3.679000000000000000e+03 4.761869999999999692e+01 1.216649999999999920e+01 1.569999999999999929e+01 +3.730000000000000000e+03 4.739840000000000231e+01 1.027590000000000003e+01 1.450000000000000000e+01 +3.734000000000000000e+03 4.912800000000000011e+01 9.352499999999999147e+00 1.700000000000000000e+01 +3.739000000000000000e+03 4.945210000000000150e+01 1.243650000000000055e+01 1.180000000000000071e+01 +3.761000000000000000e+03 4.920700000000000074e+01 9.517599999999999838e+00 1.660000000000000142e+01 +3.811000000000000000e+03 5.129599999999999937e+01 1.309280000000000044e+01 1.950000000000000000e+01 +3.821000000000000000e+03 5.108729999999999905e+01 1.192919999999999980e+01 1.950000000000000000e+01 +3.836000000000000000e+03 5.045380000000000109e+01 1.022109999999999985e+01 1.700000000000000000e+01 +3.857000000000000000e+03 4.763620000000000232e+01 1.038920000000000066e+01 1.159999999999999964e+01 +3.875000000000000000e+03 4.915100000000000335e+01 1.168960000000000043e+01 1.230000000000000071e+01 +3.897000000000000000e+03 5.408930000000000149e+01 1.087729999999999997e+01 1.669999999999999929e+01 +3.904000000000000000e+03 4.953540000000000276e+01 6.378899999999999793e+00 1.880000000000000071e+01 +3.925000000000000000e+03 4.893289999999999651e+01 8.697300000000000253e+00 1.290000000000000036e+01 +3.927000000000000000e+03 4.793449999999999989e+01 9.286899999999999267e+00 1.130000000000000071e+01 +3.939000000000000000e+03 4.919120000000000203e+01 7.587900000000000311e+00 1.569999999999999929e+01 +3.946000000000000000e+03 5.048190000000000310e+01 1.213000000000000078e+01 1.719999999999999929e+01 +3.975000000000000000e+03 4.947769999999999868e+01 1.153570000000000029e+01 1.269999999999999929e+01 +3.987000000000000000e+03 5.238130000000000308e+01 1.306220000000000070e+01 1.989999999999999858e+01 +4.024000000000000000e+03 5.436430000000000007e+01 1.347710000000000008e+01 1.800000000000000000e+01 +4.032000000000000000e+03 5.179529999999999745e+01 1.113199999999999967e+01 1.950000000000000000e+01 +4.036000000000000000e+03 5.138949999999999818e+01 1.154119999999999990e+01 1.930000000000000071e+01 +4.039000000000000000e+03 5.373310000000000031e+01 9.877599999999999270e+00 1.719999999999999929e+01 +4.063000000000000000e+03 5.244610000000000127e+01 8.590600000000000236e+00 1.730000000000000071e+01 +4.094000000000000000e+03 4.780619999999999692e+01 9.620599999999999596e+00 1.390000000000000036e+01 +4.104000000000000000e+03 4.904249999999999687e+01 1.210190000000000055e+01 1.430000000000000071e+01 +4.127000000000000000e+03 5.099060000000000059e+01 7.695800000000000196e+00 1.639999999999999858e+01 +4.160000000000000000e+03 4.874249999999999972e+01 8.923999999999999488e+00 1.130000000000000071e+01 +4.169000000000000000e+03 4.867029999999999745e+01 7.993900000000000006e+00 1.440000000000000036e+01 +4.175000000000000000e+03 4.755899999999999750e+01 7.772100000000000009e+00 1.430000000000000071e+01 +4.177000000000000000e+03 4.897259999999999991e+01 8.330099999999999838e+00 1.480000000000000071e+01 +4.189000000000000000e+03 4.814789999999999992e+01 9.459600000000000009e+00 1.209999999999999964e+01 +4.261000000000000000e+03 4.787530000000000285e+01 1.212800000000000011e+01 1.469999999999999929e+01 +4.271000000000000000e+03 5.418030000000000257e+01 1.208079999999999998e+01 1.669999999999999929e+01 +4.275000000000000000e+03 5.312879999999999825e+01 9.339800000000000324e+00 1.719999999999999929e+01 +4.280000000000000000e+03 4.921620000000000061e+01 1.110350000000000037e+01 1.380000000000000071e+01 +4.287000000000000000e+03 4.938479999999999848e+01 1.017319999999999958e+01 1.500000000000000000e+01 +4.300000000000000000e+03 4.818139999999999645e+01 8.635600000000000165e+00 1.140000000000000036e+01 +4.301000000000000000e+03 4.985020000000000095e+01 7.871000000000000441e+00 2.089999999999999858e+01 +4.323000000000000000e+03 4.964679999999999893e+01 7.883700000000000152e+00 1.509999999999999964e+01 +4.336000000000000000e+03 4.921280000000000143e+01 7.107700000000000351e+00 1.700000000000000000e+01 +4.349000000000000000e+03 4.895689999999999742e+01 9.070999999999999730e+00 1.369999999999999929e+01 +4.354000000000000000e+03 4.878320000000000078e+01 1.331460000000000043e+01 1.400000000000000000e+01 +4.371000000000000000e+03 5.210419999999999874e+01 8.752100000000000435e+00 1.860000000000000142e+01 +4.377000000000000000e+03 5.035179999999999723e+01 1.000339999999999918e+01 1.550000000000000000e+01 +4.393000000000000000e+03 5.432789999999999964e+01 8.603099999999999525e+00 1.469999999999999929e+01 +4.411000000000000000e+03 4.991949999999999932e+01 8.967100000000000293e+00 1.789999999999999858e+01 +4.445000000000000000e+03 5.176579999999999870e+01 1.065329999999999977e+01 1.569999999999999929e+01 +4.464000000000000000e+03 5.056790000000000163e+01 1.180410000000000004e+01 1.580000000000000071e+01 +4.466000000000000000e+03 5.452750000000000341e+01 9.548700000000000188e+00 1.569999999999999929e+01 +4.480000000000000000e+03 5.034470000000000312e+01 9.553399999999999892e+00 1.810000000000000142e+01 +4.501000000000000000e+03 5.065460000000000207e+01 1.076929999999999943e+01 1.180000000000000071e+01 +4.508000000000000000e+03 5.029679999999999751e+01 6.419400000000000439e+00 1.300000000000000000e+01 +4.548000000000000000e+03 5.018469999999999942e+01 1.207910000000000039e+01 1.490000000000000036e+01 +4.559000000000000000e+03 4.916440000000000055e+01 1.261749999999999972e+01 1.359999999999999964e+01 +4.560000000000000000e+03 5.049249999999999972e+01 9.122600000000000264e+00 1.730000000000000071e+01 +4.592000000000000000e+03 4.932780000000000342e+01 1.208709999999999951e+01 1.330000000000000071e+01 +4.605000000000000000e+03 5.064410000000000167e+01 1.119359999999999999e+01 1.830000000000000071e+01 +4.625000000000000000e+03 5.364249999999999829e+01 1.138719999999999999e+01 1.860000000000000142e+01 +4.642000000000000000e+03 5.289110000000000156e+01 1.172969999999999935e+01 2.050000000000000000e+01 +4.651000000000000000e+03 5.190400000000000347e+01 1.018849999999999945e+01 1.839999999999999858e+01 +4.703000000000000000e+03 4.807189999999999941e+01 9.194300000000000139e+00 1.200000000000000000e+01 +4.706000000000000000e+03 4.827179999999999893e+01 1.302730000000000032e+01 1.400000000000000000e+01 +4.709000000000000000e+03 4.999960000000000093e+01 7.598099999999999632e+00 1.540000000000000036e+01 +4.745000000000000000e+03 5.296039999999999992e+01 9.792999999999999261e+00 1.819999999999999929e+01 +4.763000000000000000e+03 5.106069999999999709e+01 9.926600000000000534e+00 1.789999999999999858e+01 +4.841000000000000000e+03 5.369460000000000122e+01 8.873499999999999943e+00 1.509999999999999964e+01 +4.857000000000000000e+03 5.355340000000000344e+01 9.609700000000000131e+00 1.750000000000000000e+01 +4.878000000000000000e+03 5.166460000000000008e+01 1.088109999999999999e+01 1.680000000000000071e+01 +4.887000000000000000e+03 4.866559999999999775e+01 9.864800000000000679e+00 1.140000000000000036e+01 +4.896000000000000000e+03 5.466539999999999822e+01 9.804999999999999716e+00 1.519999999999999929e+01 +4.911000000000000000e+03 4.882750000000000057e+01 1.255969999999999942e+01 1.350000000000000000e+01 +4.928000000000000000e+03 4.882809999999999917e+01 9.199999999999999289e+00 1.269999999999999929e+01 +4.931000000000000000e+03 4.868829999999999814e+01 9.223499999999999588e+00 1.230000000000000071e+01 +4.978000000000000000e+03 5.063900000000000290e+01 1.002280000000000015e+01 1.610000000000000142e+01 +4.997000000000000000e+03 5.097710000000000008e+01 1.234190000000000076e+01 1.989999999999999858e+01 +5.009000000000000000e+03 5.376100000000000279e+01 1.255739999999999945e+01 1.810000000000000142e+01 +5.014000000000000000e+03 5.327579999999999671e+01 8.985699999999999577e+00 1.650000000000000000e+01 +5.017000000000000000e+03 5.040019999999999811e+01 1.138889999999999958e+01 1.440000000000000036e+01 +5.029000000000000000e+03 4.947370000000000090e+01 7.038499999999999979e+00 1.639999999999999858e+01 +5.046000000000000000e+03 4.985759999999999792e+01 1.235420000000000051e+01 1.380000000000000071e+01 +5.064000000000000000e+03 5.128970000000000340e+01 6.443699999999999761e+00 1.680000000000000071e+01 +5.097000000000000000e+03 5.406609999999999872e+01 1.276750000000000007e+01 1.919999999999999929e+01 +5.099000000000000000e+03 4.973259999999999792e+01 6.613100000000000200e+00 1.880000000000000071e+01 +5.100000000000000000e+03 4.974790000000000134e+01 6.658299999999999663e+00 1.860000000000000142e+01 +5.109000000000000000e+03 5.359969999999999857e+01 1.330390000000000050e+01 2.000000000000000000e+01 +5.111000000000000000e+03 4.803110000000000213e+01 1.253960000000000008e+01 1.359999999999999964e+01 +5.133000000000000000e+03 5.133440000000000225e+01 8.913199999999999790e+00 1.860000000000000142e+01 +5.142000000000000000e+03 5.374439999999999884e+01 1.406969999999999921e+01 1.950000000000000000e+01 +5.146000000000000000e+03 5.294140000000000157e+01 1.052890000000000015e+01 2.010000000000000142e+01 +5.149000000000000000e+03 4.957410000000000139e+01 1.019149999999999956e+01 1.630000000000000071e+01 +5.158000000000000000e+03 5.216009999999999991e+01 1.117590000000000039e+01 1.919999999999999929e+01 +5.229000000000000000e+03 4.804529999999999745e+01 8.460800000000000765e+00 1.059999999999999964e+01 +5.275000000000000000e+03 4.924450000000000216e+01 8.537399999999999878e+00 1.710000000000000142e+01 +5.279000000000000000e+03 5.161939999999999884e+01 9.574899999999999523e+00 1.500000000000000000e+01 +5.280000000000000000e+03 5.392240000000000322e+01 1.022669999999999924e+01 1.780000000000000071e+01 +5.300000000000000000e+03 5.025959999999999894e+01 8.360699999999999577e+00 1.660000000000000142e+01 +5.335000000000000000e+03 5.089629999999999654e+01 1.054840000000000089e+01 1.739999999999999858e+01 +5.347000000000000000e+03 5.150390000000000157e+01 9.111800000000000566e+00 1.789999999999999858e+01 +5.349000000000000000e+03 5.351959999999999695e+01 1.266539999999999999e+01 2.130000000000000071e+01 +5.371000000000000000e+03 5.049730000000000274e+01 9.942700000000000315e+00 1.159999999999999964e+01 +5.397000000000000000e+03 4.966629999999999967e+01 1.218449999999999989e+01 1.369999999999999929e+01 +5.404000000000000000e+03 4.840240000000000009e+01 1.169459999999999944e+01 1.350000000000000000e+01 +5.424000000000000000e+03 5.101769999999999783e+01 1.135440000000000005e+01 1.980000000000000071e+01 +5.426000000000000000e+03 4.937579999999999814e+01 8.121299999999999741e+00 1.300000000000000000e+01 +5.433000000000000000e+03 4.955340000000000344e+01 6.812000000000000277e+00 1.780000000000000071e+01 +5.440000000000000000e+03 4.901149999999999807e+01 1.093079999999999963e+01 1.400000000000000000e+01 +5.480000000000000000e+03 5.157630000000000337e+01 7.887900000000000134e+00 1.800000000000000000e+01 +5.490000000000000000e+03 5.184539999999999793e+01 1.076859999999999928e+01 1.810000000000000142e+01 +5.516000000000000000e+03 5.452830000000000155e+01 1.106060000000000088e+01 1.500000000000000000e+01 +5.538000000000000000e+03 4.788269999999999982e+01 1.115760000000000041e+01 1.269999999999999929e+01 +5.541000000000000000e+03 5.013199999999999790e+01 8.317000000000000171e+00 1.689999999999999858e+01 +5.546000000000000000e+03 5.212069999999999936e+01 1.245850000000000080e+01 1.910000000000000142e+01 +5.562000000000000000e+03 4.865160000000000196e+01 8.680099999999999483e+00 1.109999999999999964e+01 +5.629000000000000000e+03 5.188920000000000243e+01 1.264450000000000074e+01 2.060000000000000142e+01 +5.640000000000000000e+03 5.355040000000000333e+01 7.667200000000000237e+00 1.409999999999999964e+01 +5.643000000000000000e+03 5.318639999999999901e+01 1.249489999999999945e+01 2.200000000000000000e+01 +5.664000000000000000e+03 4.829529999999999745e+01 8.239100000000000534e+00 1.340000000000000036e+01 +5.676000000000000000e+03 5.239620000000000033e+01 1.068919999999999959e+01 1.969999999999999929e+01 +5.688000000000000000e+03 4.770029999999999859e+01 8.105700000000000571e+00 1.059999999999999964e+01 +5.692000000000000000e+03 4.960510000000000019e+01 8.365899999999999892e+00 1.689999999999999858e+01 +5.705000000000000000e+03 4.977040000000000219e+01 9.957599999999999341e+00 1.719999999999999929e+01 +5.715000000000000000e+03 5.246050000000000324e+01 9.431100000000000705e+00 2.000000000000000000e+01 +5.717000000000000000e+03 5.122560000000000002e+01 7.105199999999999960e+00 1.639999999999999858e+01 +5.731000000000000000e+03 4.767830000000000013e+01 8.380100000000000549e+00 1.390000000000000036e+01 +5.745000000000000000e+03 5.296640000000000015e+01 1.332680000000000042e+01 2.100000000000000000e+01 +5.750000000000000000e+03 5.103139999999999787e+01 1.214949999999999974e+01 2.019999999999999929e+01 +5.779000000000000000e+03 5.073140000000000072e+01 1.375159999999999982e+01 1.309999999999999964e+01 +5.792000000000000000e+03 4.742099999999999937e+01 1.098479999999999990e+01 2.799999999999999822e+00 +5.797000000000000000e+03 5.068789999999999907e+01 1.243290000000000006e+01 1.750000000000000000e+01 +5.800000000000000000e+03 4.902799999999999869e+01 1.323850000000000016e+01 1.250000000000000000e+01 +5.822000000000000000e+03 5.286299999999999955e+01 8.698800000000000310e+00 1.689999999999999858e+01 +5.825000000000000000e+03 5.261979999999999791e+01 1.278669999999999973e+01 2.089999999999999858e+01 +5.839000000000000000e+03 5.338810000000000144e+01 7.228699999999999903e+00 1.509999999999999964e+01 +5.856000000000000000e+03 4.854509999999999792e+01 1.335319999999999929e+01 1.319999999999999929e+01 +5.871000000000000000e+03 4.994619999999999749e+01 7.264499999999999957e+00 1.650000000000000000e+01 +5.906000000000000000e+03 4.950619999999999976e+01 8.558500000000000441e+00 1.630000000000000071e+01 +5.930000000000000000e+03 5.464099999999999824e+01 1.002379999999999960e+01 1.490000000000000036e+01 +5.941000000000000000e+03 4.767540000000000333e+01 1.246979999999999933e+01 1.440000000000000036e+01 +6.093000000000000000e+03 5.321390000000000242e+01 1.047039999999999971e+01 1.839999999999999858e+01 +6.105000000000000000e+03 5.431940000000000168e+01 9.805099999999999483e+00 1.569999999999999929e+01 +6.109000000000000000e+03 5.338369999999999749e+01 1.437279999999999980e+01 2.089999999999999858e+01 +6.129000000000000000e+03 5.105930000000000035e+01 1.442660000000000053e+01 1.939999999999999858e+01 +6.157000000000000000e+03 5.364099999999999824e+01 8.080799999999999983e+00 1.430000000000000071e+01 +6.158000000000000000e+03 4.922469999999999857e+01 1.060839999999999961e+01 1.369999999999999929e+01 +6.159000000000000000e+03 5.295420000000000016e+01 7.319600000000000328e+00 1.519999999999999929e+01 +6.163000000000000000e+03 5.416539999999999822e+01 1.035190000000000055e+01 1.619999999999999929e+01 +6.170000000000000000e+03 5.201919999999999789e+01 1.472540000000000049e+01 1.960000000000000142e+01 +6.197000000000000000e+03 5.186639999999999873e+01 9.271000000000000796e+00 1.719999999999999929e+01 +6.199000000000000000e+03 5.424839999999999662e+01 1.304180000000000028e+01 1.930000000000000071e+01 +6.217000000000000000e+03 4.924060000000000059e+01 6.935100000000000264e+00 1.860000000000000142e+01 +6.258000000000000000e+03 4.768449999999999989e+01 9.440899999999999181e+00 1.530000000000000071e+01 +6.259000000000000000e+03 4.902100000000000080e+01 9.603300000000000836e+00 1.430000000000000071e+01 +6.260000000000000000e+03 4.933279999999999887e+01 9.704000000000000625e+00 1.559999999999999964e+01 +6.262000000000000000e+03 4.876950000000000074e+01 9.873699999999999477e+00 1.430000000000000071e+01 +6.263000000000000000e+03 4.777380000000000138e+01 8.821899999999999409e+00 1.359999999999999964e+01 +6.264000000000000000e+03 5.141400000000000148e+01 8.650000000000000355e+00 1.580000000000000071e+01 +6.265000000000000000e+03 5.236129999999999995e+01 1.238669999999999938e+01 2.219999999999999929e+01 +6.266000000000000000e+03 5.203040000000000020e+01 1.096260000000000012e+01 1.860000000000000142e+01 +6.272000000000000000e+03 5.084259999999999735e+01 1.025179999999999936e+01 1.630000000000000071e+01 +6.273000000000000000e+03 5.250750000000000028e+01 1.185510000000000019e+01 2.050000000000000000e+01 +6.275000000000000000e+03 4.867049999999999699e+01 9.462699999999999889e+00 1.369999999999999929e+01 +6.305000000000000000e+03 5.120609999999999928e+01 1.049779999999999980e+01 1.930000000000000071e+01 +6.310000000000000000e+03 5.410490000000000066e+01 1.382390000000000008e+01 1.869999999999999929e+01 +6.312000000000000000e+03 4.953139999999999787e+01 1.064179999999999993e+01 1.480000000000000071e+01 +6.314000000000000000e+03 5.105069999999999908e+01 1.330030000000000001e+01 1.769999999999999929e+01 +6.336000000000000000e+03 5.001319999999999766e+01 9.653999999999999915e+00 1.590000000000000036e+01 +6.337000000000000000e+03 5.176630000000000109e+01 7.519400000000000084e+00 1.730000000000000071e+01 +6.344000000000000000e+03 5.039399999999999835e+01 8.142300000000000537e+00 1.950000000000000000e+01 +6.346000000000000000e+03 4.820700000000000074e+01 1.120350000000000001e+01 1.290000000000000036e+01 +6.347000000000000000e+03 5.005789999999999651e+01 1.029720000000000013e+01 1.710000000000000142e+01 +7.075000000000000000e+03 4.870199999999999818e+01 1.184929999999999950e+01 1.309999999999999964e+01 +7.099000000000000000e+03 5.201109999999999900e+01 1.039659999999999940e+01 1.660000000000000142e+01 +7.105000000000000000e+03 4.783500000000000085e+01 1.265479999999999983e+01 1.350000000000000000e+01 +7.106000000000000000e+03 5.207139999999999702e+01 8.456500000000000128e+00 1.789999999999999858e+01 +7.187000000000000000e+03 4.976359999999999673e+01 9.406299999999999883e+00 1.710000000000000142e+01 +7.298000000000000000e+03 5.452680000000000149e+01 9.042500000000000426e+00 1.509999999999999964e+01 +7.319000000000000000e+03 4.873740000000000094e+01 1.073930000000000007e+01 1.450000000000000000e+01 +7.321000000000000000e+03 5.115070000000000050e+01 1.133210000000000051e+01 1.939999999999999858e+01 +7.329000000000000000e+03 5.054670000000000130e+01 1.228630000000000067e+01 1.600000000000000000e+01 +7.330000000000000000e+03 5.146329999999999671e+01 7.977999999999999758e+00 1.700000000000000000e+01 +7.331000000000000000e+03 4.860990000000000322e+01 1.026740000000000030e+01 1.269999999999999929e+01 +7.341000000000000000e+03 5.009000000000000341e+01 8.786199999999999122e+00 1.750000000000000000e+01 +7.343000000000000000e+03 5.061990000000000123e+01 1.348160000000000025e+01 1.390000000000000036e+01 +7.350000000000000000e+03 4.910880000000000223e+01 1.282310000000000016e+01 1.190000000000000036e+01 +7.351000000000000000e+03 5.331750000000000256e+01 1.341750000000000043e+01 2.110000000000000142e+01 +7.364000000000000000e+03 5.168200000000000216e+01 1.230419999999999980e+01 2.010000000000000142e+01 +7.367000000000000000e+03 5.196430000000000149e+01 9.807199999999999918e+00 1.889999999999999858e+01 +7.368000000000000000e+03 5.100070000000000192e+01 1.036209999999999987e+01 1.600000000000000000e+01 +7.369000000000000000e+03 4.916230000000000189e+01 1.036609999999999943e+01 1.319999999999999929e+01 +7.370000000000000000e+03 4.939099999999999824e+01 1.268379999999999974e+01 1.219999999999999929e+01 +7.373000000000000000e+03 5.359839999999999804e+01 6.702399999999999913e+00 1.430000000000000071e+01 +7.374000000000000000e+03 5.208129999999999882e+01 6.940900000000000070e+00 1.469999999999999929e+01 +7.389000000000000000e+03 5.274609999999999843e+01 1.384270000000000067e+01 2.139999999999999858e+01 +7.393000000000000000e+03 5.144930000000000092e+01 1.425329999999999941e+01 1.989999999999999858e+01 +7.394000000000000000e+03 5.003150000000000119e+01 1.197450000000000081e+01 1.359999999999999964e+01 +7.395000000000000000e+03 4.865950000000000131e+01 1.253880000000000017e+01 1.380000000000000071e+01 +7.396000000000000000e+03 5.050840000000000174e+01 9.224700000000000344e+00 1.350000000000000000e+01 +7.403000000000000000e+03 4.779549999999999699e+01 1.003240000000000087e+01 1.230000000000000071e+01 +7.410000000000000000e+03 5.075130000000000052e+01 9.022399999999999309e+00 1.650000000000000000e+01 +7.412000000000000000e+03 5.000829999999999842e+01 9.423799999999999955e+00 1.509999999999999964e+01 +7.419000000000000000e+03 5.066100000000000136e+01 1.207559999999999967e+01 1.810000000000000142e+01 +7.420000000000000000e+03 5.110439999999999827e+01 1.171119999999999983e+01 1.960000000000000142e+01 +7.424000000000000000e+03 4.777239999999999753e+01 1.290729999999999933e+01 1.639999999999999858e+01 +7.427000000000000000e+03 5.401879999999999882e+01 9.925499999999999545e+00 1.660000000000000142e+01 +7.428000000000000000e+03 5.041669999999999874e+01 1.081559999999999988e+01 1.700000000000000000e+01 +7.431000000000000000e+03 4.801299999999999812e+01 1.155240000000000045e+01 1.290000000000000036e+01 +7.432000000000000000e+03 5.264229999999999876e+01 1.066269999999999918e+01 1.919999999999999929e+01 +1.367000000000000000e+04 5.150880000000000081e+01 6.701800000000000423e+00 1.730000000000000071e+01 +1.367400000000000000e+04 4.929429999999999978e+01 8.905300000000000438e+00 1.550000000000000000e+01 +1.367500000000000000e+04 5.208180000000000121e+01 9.407700000000000173e+00 2.039999999999999858e+01 +1.369600000000000000e+04 5.159660000000000224e+01 7.404799999999999827e+00 1.710000000000000142e+01 +1.370000000000000000e+04 5.133290000000000219e+01 7.341099999999999959e+00 1.580000000000000071e+01 +1.371000000000000000e+04 4.857339999999999947e+01 1.225760000000000005e+01 1.269999999999999929e+01 +1.371100000000000000e+04 5.068200000000000216e+01 1.151500000000000057e+01 1.839999999999999858e+01 +1.371300000000000000e+04 5.108990000000000009e+01 7.628899999999999793e+00 1.650000000000000000e+01 +1.377700000000000000e+04 5.224669999999999703e+01 1.095919999999999916e+01 1.980000000000000071e+01 +1.396500000000000000e+04 4.826389999999999958e+01 8.813399999999999679e+00 1.040000000000000036e+01 +1.500000000000000000e+04 5.079829999999999757e+01 6.024399999999999977e+00 1.290000000000000036e+01 +1.520700000000000000e+04 5.128349999999999653e+01 9.358999999999999986e+00 1.710000000000000142e+01 +1.544400000000000000e+04 4.844180000000000064e+01 9.921599999999999753e+00 1.169999999999999929e+01 +1.555500000000000000e+04 4.787610000000000099e+01 1.058489999999999931e+01 1.080000000000000071e+01 From 537b96fdf624547acda66477c6176a86a4264850 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Thu, 17 Dec 2020 16:43:17 +0100 Subject: [PATCH 22/39] Field: use field_dim for mesh output options --- gstools/field/base.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/gstools/field/base.py b/gstools/field/base.py index 02908f98..ef4fde7f 100755 --- a/gstools/field/base.py +++ b/gstools/field/base.py @@ -122,15 +122,15 @@ def mesh( pass if isinstance(direction, str) and direction == "all": - select = list(range(self.model.dim)) + select = list(range(self.model.field_dim)) elif isinstance(direction, str): - select = _get_select(direction)[: self.model.dim] + select = _get_select(direction)[: self.model.field_dim] else: - select = direction[: self.model.dim] - if len(select) < self.model.dim: + select = direction[: self.model.field_dim] + if len(select) < self.model.field_dim: raise ValueError( "Field.mesh: need at least {} direction(s), got '{}'".format( - self.model.dim, direction + self.model.field_dim, direction ) ) # convert pyvista mesh @@ -158,7 +158,7 @@ def mesh( offset = [] length = [] mesh_dim = mesh.points.shape[1] - if mesh_dim < self.model.dim: + if mesh_dim < self.model.field_dim: raise ValueError("Field.mesh: mesh dimension too low!") pnts = np.empty((0, mesh_dim), dtype=np.double) for cell in mesh.cells: From f1b36c74794419ac61eac86eceb40b6456cdbba3 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Mon, 11 Jan 2021 01:53:15 +0100 Subject: [PATCH 23/39] Generator: mc sampling now around inverse of 'rescaled' length --- gstools/field/generator.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/gstools/field/generator.py b/gstools/field/generator.py index 213f3841..2a66552e 100644 --- a/gstools/field/generator.py +++ b/gstools/field/generator.py @@ -212,7 +212,7 @@ def reset_seed(self, seed=np.nan): rad = self._rng.sample_ln_pdf( ln_pdf=self.model.ln_spectral_rad_pdf, size=self._mode_no, - sample_around=1.0 / self.model.len_scale, + sample_around=1.0 / self.model.len_rescaled, ) # get fully spatial samples by multiplying sphere samples and radii self._cov_sample = rad * sphere_coord From 303c7618910183b87d7ca14f57155f7d7f4e3d47 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Mon, 11 Jan 2021 01:53:47 +0100 Subject: [PATCH 24/39] SRF: conditioning with latlon needs field_dim --- gstools/field/srf.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/gstools/field/srf.py b/gstools/field/srf.py index f63857a9..fd984936 100644 --- a/gstools/field/srf.py +++ b/gstools/field/srf.py @@ -188,7 +188,7 @@ def set_condition( """ if cond_pos is not None: self._cond_pos, self._cond_val = set_condition( - cond_pos, cond_val, self.model.dim + cond_pos, cond_val, self.model.field_dim ) else: self._cond_pos = self._cond_val = None From 42d7a33457a25efbe8ce5ad81329054ddc426400 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Mon, 11 Jan 2021 01:54:14 +0100 Subject: [PATCH 25/39] Examples: apply black --- examples/00_misc/03_precipitation.py | 72 ++++++++++++++-------------- 1 file changed, 36 insertions(+), 36 deletions(-) diff --git a/examples/00_misc/03_precipitation.py b/examples/00_misc/03_precipitation.py index 83055b5e..429d6971 100644 --- 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) From 54d48f580b3f97dbafbee808f4354f8325c1c400 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Mon, 11 Jan 2021 01:55:07 +0100 Subject: [PATCH 26/39] Example: polish latlon examples --- examples/01_random_field/mesh_ensemble.vtk | Bin 0 -> 19209 bytes .../08_geo_coordinates/00_field_generation.py | 2 +- examples/08_geo_coordinates/01_dwd_krige.py | 2 +- 3 files changed, 2 insertions(+), 2 deletions(-) create mode 100644 examples/01_random_field/mesh_ensemble.vtk diff --git a/examples/01_random_field/mesh_ensemble.vtk b/examples/01_random_field/mesh_ensemble.vtk new file mode 100644 index 0000000000000000000000000000000000000000..0ea0d7bc11249afb93cfec8632e0008350fc391b GIT binary patch literal 19209 zcmeIac|29$7e9Q>Lz7fe##BO)j3u>Bg=VQFiIO5ol8`1z8c0Yq2#Jy@G#H9(nF*O^ znXY-dTz7C?{LcAAeLmmk_x%0*^LV|^>#X6u_t|H!z4zH`t$nSOc6K~3t#iupl-@a8 zOKBrZhs)WI$79S3X=X`{D#OEc{$`y={?(7rt*Z~Y&+*YD>-&wu6U$3N|V2Lu1eIcBw}@gI5A44(6k zoLtbPw0+YpApaxzBZuGhM4mqB9{skHnGcikgxC8MlArzx(rl=gaIm@E1Qb-n1iUTmJh#^S*8~?c|@;&dO){lk!LY z`+ehQ`t5H&{I37~%sTV8-+tp0(C_~IgP$3Pk?8lh`qR!C`P-TK_~*T5#??F2`cF&# z;(tFkIz4Q*g}?7JR8Tl*x6EuOJGs6>lLN|*i^n`8}Psj;c0l$d)2%6B({TENrvwGrA z@Pw6S-=Aoc%nVP^M4d`AJV6s>*3a++O_UXx;R!r}oh?t`3G8fn0#9IP%M&<4w-UdI z=MyxcQ=Z`on!qoe;R!vVo5d4)LN|*in$U~>iznz=J#nAa5qRQ{Xp$v*igqft@W+;0T>Gei6?jXhNqt!xJ=t zm!07WJ)x7H;R!vVo5d4yf?k1NM12HJ=-2#2>+49mjO<-rs6LoyH7L&;)+|3{U6@ z-7KEa6S`SEAtz`t{37ZjXhOf>Upzt2>WMqS6JKcpPy7*WT7+K&PtZi2nln5>6J<8e z@B~ei6`J7*Jb|4pPv8mcY$v*igqft@W+;0f$(c>+i1`0-CSCujlu zBI+Y(LNER=o}g#-#GT+d{`nf5G(2XYOMDUpaS>tG|Nr{`s1JU-StkiM40p8z^>*v- z(~;IYXK8zeaOb2qs0i9CS`n&l|^-pzf>cVj?Xv_eMR-2tBUeLJ+e zhJ<=F9Hp0TS3(UY#!(7uJy5ru>a)jVvZ$?Qo%g6&5lS+7H!dgV280Ar}nYNv^+6ppQWB*@2q9Orb~JzKlHt`;GP~I!FU{--JOqk(5V7|BcA=%SwnM!K>v+}|r7hPR%W^EfXS*p2!tUZ@p8%v4&}I!-_MY3TK| z<*+d@14nPV=!wFF_P&USk6m!f^^KIbI|&WAXYq1%(J>C!!}@9Tr!Yx9?z+rA97ZcE zto=e?05wYIqNeCRV6Zo>gxfs8T6>E>^inJ|>SH-tEz{t%rmCA2LVch#LyVqJ7 zT1rnnLg~do$?nZFn1c`Jmw}nNfL?KVz24rp zJxo;z1_sGVo%hieRBrXK!?W5Ra?hz*A5&Qcm3!&36CrAU%&+oa^N*Qrx%xu(In?`G z9zUy<1C$Mc6wzo-H1IuRVo?)6YP~aXMrf4>bU5xww|N!<1K+-AWjVKEB)((Rckz`N z$LmbzMaUS&_g>5Q4oZj4nqdRKW8$i$73GyAb4U8x(>efqD0m?gCcC}9oen%yfkVqt$7vVi0{<@PqK74J=A=NESi(*2 zXUAbUf@^E^!Ot*0bobDeV@AMSeb}v|ayL+~cdg3&ehk%CsN@POe1`I?s^v>hm7;ci z*cyFY8+Ey#W$2d{BWg)&*3&m<&`@FbKIsx+pe3Yfb^jEFnxHH7@#_^Rb#+)#LK`5*JE{@47| zKav)da83g4$~l|dQ*JQ*<6^cjGQl{4Cbpd#@qv-7>$M5-UO@F-`kXVA1ryuvT3mc3 zhjD4VIUCNQjgelk;GOU~V12qd74|F`m_68b*KPhVuDprmAn6Db@o|RS@(-ZwnVq0= zw=fiO>M~l^gh5Xs*UnX%hR~betu~&10!E&ymOi#)z~sTyJ9n@4VO&1F(MDsP80TH1 zcUx-F8+S z7>^gmSdHF>>1~&|PTGoN+}`qopJV(mo}eoUu8R5?$y;ti2R|EW7&q(`Guyuo2HKeun>O5mYP6OStz!Xg zM^hdee!78?IB4G1jiKjM(pW-O5Hz~Ief!xv0CGbbjpNvH zz}{>r*U-ERI+apI))Xv+506@{6^~=6L#0)n&Q*h|Ru&dF-%Urk1(hww^(d%yo>bt6 zxjisE{+asg5#Id~!#m6P@4>JGIe}02DJoSvODSIx1jz1@CYmM;>Y=J8HtT>bZUsS5cUrGC?!FgbfEjJ z*|4Oe75c2`W_>h;gfit82W2u6VFYg0C2u|p?1{Foy|>~q?w02}b+maKta7pw~>6@sBkVXm^w zgEHts0#|+WpF&+2N8-BjUvS-yQ$abw4c)$WeYn5h4K+XPF!;J{75b^9E4*m0G_YQ6 zQFfF0?QtkB= zuyzGzww%L&^7WOh7F7vG#=ST?#X{lJ64}cA2}~GW!l>yq@CVxB4IhKs=L4H^SL#!* zD71toiPm2GyY*B5ulboqLI-VDP+(|PTOnf}3#P2UNlk2c3al+U8r~d{C@epKWTTD2 zpnH^y%Kc>+X(~AR#i$RkRymSu*FAs<&4~4v=lcMoShi~EDtDmCxxGkv@D?T%+u|KZ z8h~A_uGaaa6?F)@H~adsQQn=vfMP2*=&rx5>pNx#?FOe?o_~G=?P&)zZ%sS}+J{)q z7Xh-svTrc$ahw3U`yxlSpax8sJ1KBo!)IAuZzB6A6XOY!coB6V72^;(rnIg_7np~_ z{A(*oK>Oyire0ndMqOM@$&p>a2!G!mwOJn6cSyJEt2;5$e9=cMMKoZr`n&$y)3zvW z^VgFJdr7FdoJZS7m>+0o`o_zNbybLmkd_x{X&jsKV*RVh6&y<0T~V|%`7svc&*=n5f!0ew&CE}Un>?f3@e zru02lzKeHvfver;hXXLUb;$O}!A?~9by%F%7Y5}X63(@V8-`g&qC1)ffwAD`{^1&H zpw1hLf1;BD1KOW#)6!m`v4uAiOKp~;k)gR9Y8#?}{8dpTMt%(AYcUORIjDt^zGV)% zV7Fjcvx)CI*A_I?p0-y}h>R#nJKR>~_MuAYlU2TRxPfjZe=#{R7N{LN1e6w3!oZoR zIlQj>5IN1|CQpCg+SLa)(Vim z4wL(iq_A&X1FBRKWrfr{_;udg{LrVx(5uc0Y@MS7!}-cJzqHQ*d33vE$c9Bg70z(q zQuPif`6sSP=+(pURjs{7bvt0V-r_@`z7LEFbZdXtH3Rao)Zl9I5ny;2WxgVjVdP2d z^80!pQ17z)c~+HWh+h)J!v=(H-#DG_YEAI(tEnMloO!W;K zzJB1l1{FhnnwbX|%r!#cKc1%7-3$9;e$D@wADD0J9!g6J0;@x&WUIac#=FyN_}4~7 z^g({P|BCPCklguNU1qm23@h*crMcV)Xv}3ME}lskX`zsFzmyujIBhQ6WRMKg+cmq? z>ZD;b!D#K^byFCZQP0i2)PeDwm^g8{Hwd}~*qbgz^1~QQZbI}b{;dn%V{$>dYUi$r0vQOG>R#R}sRKK97ScA-S+wGYB4cLn5W+tP&7-@a!sVMnS7^>1<9e5xB=&1?5WZ`cxvZcGmuWB`9`!CuO z8*B|_^wG~htO}s&uCauInh4~rOnZHfTO9AV+t-%HS;FL@w6%}e(HLi$NdimX4%pjc zytH1ghlwO^x}W#cwqp# zwEP#Wox2a;3dgKIa=;4T%DH)#tc?NoTGOk;ayA%gXxr}0d3P|N}n#3}N!+Ah)H=F47 zQH5GJ?Zn=|0O%r>7wh_jLU%K#(t>SKkZu;I+4d?G^|v?P3(&oXCSx6#m%8-P#E-?? zf-$D3@1{eb$7u!VHc9AnIPd_PyB|7zx|a#L8~NfHKh)rRdd~r=^Xs4@`-c0^lm1Zu zQACd8S34@Y@Lu2Zp)qRW6KJaPvPRvk6(>yjbPyx#b4qHJ7V7s4Iwxjki#nZ|*4xJn zp?g|uj^A4~yuIZDCOZ;gAR5cnu+Tys^7$){b@D>9$przcB|#L8IH+)=FsjS6@C5PZ5t4rLxc!RHZV4Zq&VitlnUhPu?#Z+EV8 zKm#HYOAovWMMd=887DT$p$?0J9rx5@ATMfaK+$Y5q}L1Zv*!sx$6;ZvqsebjWx3nA zMQ1R`IoRyX`O_Hv?BnxnZF`PbM}nv7Z_Gs_6l(yoBm^d_$9LI z-P<{N7@r7#O5d$_P?=^&zTk{I9JkkfGs;s1_qhhKfx2$c?=@#|aNiPCSYL0h>r4Sg zu9iZ`{o5FqD?j~?yE7*6L-)tCffks~x%pCm@p>3*P`tdoHV(SdA-eIJD{3=Z`)FxD z2~v9XnCs^jz_46qUq)&o)DEarFKiw~32X8u59^Xpg=M(LiDOMLNO8NUsAdinHFU0J zsVEeR=q!`Jrv^EKj8A(R1<*?J?#W7Qg3;rP$wHpP5W7LD?o>k<3}@M>z8^P(F`vW; z%M~NoT#coI@s+-KyRKJw$n-#!MxV_f$7j^5n_JD;k%6e!=lRoHoKf%V#(neIUx7UL z`lO(@Ba9T^mbSqxVIoP;F*NTC(4v};j?}ZEqPqrJtZhOCi$x#r{rMc|HGzEnQl&65 z)fKt4t`MsF_l5{G;d`PtE#J3G2V*4m!8)oD8K@1Lm~G9sVXWp1*-}3h8XFJ)HWI>AOQD(Y8jF5bYQGgs!Bhi3F9p)&`sM=^~d}h z{u}lK+6Ml)Z^NsBu}F@`Etdh)6?+2;sd&7@uFg8c7WX`!8&utVdIae0gUW&`Bp8zm z>|VGm8S0&6ZG3iaM|s5u_1<1~fmZ#?P2KKafxY28Q#YOuw1=Bw14aU22XzZaW{@z1JW3LY3t7C(F z&bz1XFsZ*Po?pZTdbLgncz7*$`Jr#TMu_Wo~(k-RW-8ql82y0 zaau>D!Wb9_m2Yes*$C66E!EpvgD|eD>X9$2$QVcST>E>M!(ieA&)|)LVi>M`Q@+`? z4Z2Ge2RhIVXkYlfBHW%Is5X2r&ykHF^Wl`~{1gnjRxkJOeTi?5-lYv`uIK|=?)rfA zcut^y>?_h)KL_IsJAZqr*c_mpJ|Zu7)*2=|N9!JJ(uCoSbk&1)RxrwPUPzLC2|a9S zlhC^-VBBiZbYr{@#~`{AZ&a57`c-h>k!S{#l(3Z(>y}~MQ%WZVH1A{FA@=?xfech7 zet47X&=_i8^eW`a;UO35brlSTtxyeSR&XX5^56-diG<7 z8Vc654nXD@Os`42`1(^JjGnPN@#_9V$b1&|dEZNG)PeN<8{524zLZNxTU zi|qUU=(rKafuyvWGtL5yTQKTkv;ojY!a>f34wLmGi@xM~W8A(h%gsvV7*BD7O7i0X zjJI=6V4ag8Q15@L@7+=hgT;-s`ym<7yGe6(wE2AKqzLl$`g=k4rBiqG-3Y zFbFy09)_>8ozMJJZKHo@Y5vIy;7OYQOIG0Yzp?^&W`E|078E`n^B9RCL+2M}cw4Rx z43u81tx_FK$(XE^SP#B*zOzRnuZ_w!{Nvj^N(AvI#J}Qu;0b*Fj*`4DHKV znhi7$Uh{;7JAlUC2J(zwzx@8I9&ZUDFl6P<9_srYnvQN>#$>y2y5?Fi)6TwNdq|<^aBY z1h}94rhsv-lxu5)2%48W|13UBRg{0lSs$F`mQD2h2ZM zLv_b~A=O_pzz9_k-hG(|)cA+taTdpz`74*yg&_3Tf3jTQz zG5frc6AgK*HnV)(5M%HBt>h#})b;CZ(8sQ?Xh3JVCB{xt@X2SP-$I8$xm!!r=8%R(udKO8{zlZ=`#AJx1L-l#WT*>=?`E!1W3ZCcJ-2Q{nr z%X~VKipXTy=M{bNsNpg*Tdyw|WpAkEqa1NTRUU=f4_+3d_R@Z_?yJwy*WH6Lfi@KQ zZl?9+o7Zzxs}d7zJT``Yv?|s5Dknn4*SDt4^;VGVStpmj@d2R4d|E9(>VO>K>t3T3 z1UXCUj-1GIfp)uLaXT(oXuZRTT(o}^^uE}vE~#%0o%?S|K2%nLPMZ41Pw72^v(tiD zHUBj&T0FyrPdf}icF8x=D2X5X`5rzuj1^hZSqy{Xg-Zf zsxn&G3gji18;D6Yffhjyp(88R?VEwKhibC^!_)f732XMt*`Nv&QkcvxK^r+^ptIhJpA z;9lm0?-p)zmq^tn(LKNkbozHbs%2h{+xc?oy1 zSP?Lqf`>UhTYB18nsspa#?#$7P_rLd;ST?l%9=7%`EumVPE{~f@08MJ_|4P3=#Fq*M6 zc5>GtVB9fXzHf9RjIN4sYW>j;%uHvlOG+`o)b_g|(Qy`7ZU?+mt>d9G$Z4cDpa?$K z+kK)$F`*--PUi8p7+~6c0=9I<gn z!(17s6^1Y>eyC`(?hf>^=|=A_hCO^Wv8HqJW6+mmT_SRN9FZfUZb;hZp(_2`qrWCM zLgxz2QyCiWsAu2Rig(g}@K{2d^7I)AQ8=o0J*YK7#Wn}Y5z_Xk_nhLB7gv+fu&5$? zva1STO`do^j_O3+D;CR>m1Ph$q{R9@?-UxjRS=vqM-Js^X#}5>uty^&=LMzvoIupv zv71q1Whli-o<_>Reb0gJWoi1-h;n_6@4D6VQRfrO`zg()Xk^uvU#leQP^XfUw*B}X z)MvGL@VbB=>iha;%ia7M)V1DAvUn2TwCu9-`*wIPzNy_iUra3-Tlo zNo8&{HguOHmrA!f0&6n(cHH%3=ug}}|EK3sJhFB7&6Cm$=+8NQ>D*6FnAU$aS!S(- zak`y1q72jkTlUEXUT!AzDc^r5>Uj*~tiH6`KduzxkTg8LJHv40ol~%HD za=I{8*cwtB+XEx%uFnLf4&mcS{ob^aG0?9Zk$SXCAI8&MR-Yi>h0!W5E`D8Qpi1|f z4vR>_@T6Rd^Y`m8ExIdJYsFrSGyCAq_k%lO;<;9MUdj-Rg!*NZ`s^`IU((JG{6{bj zh4<;_8}xutb4a|*70;qQ-J*C}Dgxu6ZmmsniG&F)ca03UB8)Rrfqiw$8jQ4aRZ!>> zJOEtXBl$>63MMUIcs`{L02}eGvX;l4zX{9kOxff=(}KvZ{hbyVN6w)mhdORyoXv~g z2(2-IX}|Rs9kgU&eDj5o%#2+a7oTLv*83YV&LvI@vAI`aa!*2IYDgT$+4@A|?c79+ zbN}j)jXrn;>|Cpl=d-66*V+%-PV?8nbncv_<(0VporbSCg$0v>Xxd`+T#WQ??v=4g z3?uDh^rwdk0=qi)42O$9#vw(*9JN>&XDn%9-IofCLu_JC_m9uO&J6cIcJ(X95ich{ zxkVA<$P;eSyoCF(GH14#?H$5MpN-wyRU%=U)6}i`MFLC@eaVtpjjx{1?ULswMPZ!M zyN9+oMgn`8y;6)Q-o5(Q)B^-Z@b{MQ-du$Nqd!>8&a@g|#k!AdX*!5;VU50Y${UR5 zx@GlzAvd75=;nJIqrz0!YB{&;RE(#5K-qoT7DiSy^P4X_2<$xaE*^>~#x3!~NrJWm z#y00)2@6!fNX6|RXv{&G4$K1p5$GA3d#*j!yUIShCbRO&D1{Jr)MfE^4FY!uD*Nu z#xkd1IC+1P?M4}l_53b11Q)OWnu3i#%bI|iAH1iW_jQ^i65Un z*M7{gww49EyOa3g1B#%&Gp>71rN|~7#@s;5};=Xx~!`X{WC3y ztmEHl0lm2rH$LC2f}#W5p98k9LQPL}>TZ7Bh&o!9wuC93fNGy-8yBT_p&v+Ai01VG zb?lH@A7e2F`tSRW%6JB4v%tqsuxBaHKT;q=iuJgxRrUo}0%3gh5}{?XOPI!+0%uFP^%c0o1a`C3B+;Fg~}6u@UNqhtoPB~7_VV`|9)+5VBO22yew*f>3wfp_TO`azV^HV_-F+q zlEW`d7RN$ciK#j@d=*UH6UfdR+=FpPbY1sLa{GHtavmh0i>!OH~#iz--Cjwh}IVbt(9T;5L+5TS11Q;0}(k82Ffz`iVdgFK| zj9gJDdBi;oEW6%DbSxAj9jQ}ZqC5&yZi|=ZAFYMag7)gmHOeq~_Cnr0_bV`69T9$H zEgy98F50y@z!=7Aoo)G6sNrc+gCwD`Fd%O|{8K|(2}YCr1@G_UK}bTC+9`}dWAzQQ z;7tbTN1y3i_pSoee(mY)=+VEMC?XsCuW7;5-D|V<%^0v&O9=3OW&!JD-qUICQlQ?Z zcweot1g3ZnU*l0dV8vE}=f^9+JW77G`aT|!B_}%n>{7tbU~CrLeGEU@H0`!0u_Fk^ zLze76W=nyIoFV_65%@jb{P^8lI>1(vt>0m%4D21p$o?U3@cf0x)kAg{fbMdpzE0W& z#vd&5ja&Q~kKAaRSQ78k#jCGUlCc$^^z@CO6yd(U8| zFJPSN%9)2oOfhb|7@j3u@fhisps4FoGmJZozVvX*N{r){kiKv|FUGAKLw0&Gg>f1t zAGvlH|C(_We2;dR3uC7%6od0cfaYITmazB(O!2K>lQj>6sq9Y?+&A2SvQ5XQ`$Y!O zjcuY5s9G?U9JkD^pc7g)ZX3&6Yz_lITHG2c@%Q^8AH1MW7%Eo0JHo7y1}bL;a}UQ9 zDtcRDwNf3Aq%v+*U(-7e>B;Nj7iwZC=z*el5g&%9`IK@3?%G2A{^o+%Q3`4cusNK5 z+yMPDrIn3sSp}3z8QIM;3MglD?0MIQ9@KTYNq?Hp ze0O;;arP>V882Nlej@>TUe7V%zibM$rAsKMb=-k+`1=vx!b+eY7!;+n;$Ka%aHDH` z*8s!DB>i~kPM8cxvsiM`1z6td`aEB7-=j+cow(o*%!XnPsKvh~r}GYfx@`xn8ht~R zA9%!S;5&c*hGv*l*N;4lpY560RKwh6!GsC<@ybWiUNCtzA_i+ug^34GdI~Zufqh^~ z{_I|S+`iY4=U3o|2|V8oi63EHb%9peZ5bH5<7N0}&-kCO1(BuxJ1tP8NxnV%moyZ) zTXy-J#eIsJGf!IH=HM=?-XeeDC@8nqeJ9@%1U)4ov>K6e7<(ok-jr+&WII}e%v~#> znD35X>F5fye81E?bo}AV&PzCyBm*_rp-$~61_Rj_&+f4FhY?LHyQDe$fg&`K#(lO5 zD3_1F-+Dh0%CbZQYA$(0PtiBSjUz5Fu)JEu!f6Dkx1!?Dl%+%WD~}z~FCM}`kkYL= z{6a9Itk$wzLLY{xo7KI};Q{cGJ;tmo5>R#Ydkvpu!svR3#m;Odj8u_b1+vs(bV!G5 zPUjem91>Z(%~b;^dEI(X>#SfP;BASFrvys4|5!+zVGg~$2da%_@v|czH#0VwzlM>^ zymtIu>!5gk_>!Yn$6!?6>1bb>2++^3a}#+TkK*2#>P2_qX<3g_Mq5$~uNqusHI(Xz{>m5?;D}G=;UMRQnpdj?-eON5HMhhzM-ue2gy9l+lcyHzv z2u0)rrx+2XS*T3}l4JSs{8`!-n>LmY>Xz7lZTU?r)axBqXi?;a2J>z`=@r#OEt8!u zPtiP3UxuHYZ=?qr8+`vGZb%8$?8`YgH{A!F?O#( zNMs^7SBiqFIb-#=yzD_`_nvF9lysm`f$h6I%N!bV+*dx)<;8a!7W|~0mvH+Qs7Ger z`!g+w?ET+q0TZ!vvswlA!cdN1S&;fQkz>g>wXci11{GZ zK8uD4$He!&$Pd`UB+rBiYhbVINc!q8158_#cVPf`X($G(X~H@%DG|J9xz;p(vhO)` zNZi0UE`{#MU2TGKWw_HtOz0TzyZx0D)h@u~eUi0)dmnzbrCgFO`xvHqpX_|iV~ue> zI&JmZ1CQ9G&{u4(;m!{8Md-QySBG7yu-+Y1(yW#Wgf#1D% z6VqT}VU-N@MK*r24_*u8khh=yQ`( zDXl>JcxTDW%pf3#xu}DfB$%Oc2ZMDaV0L(ZRdD%^PqX}wj!T|GdCSVZ29a_wrFl{Gqau@%VlIpz0k@Q9R^!EJ}9gV2=Zcv{6J9J(sD zz0q+hg35vwA5OXzz<9jA*TohZ#$m>L<+IcWm{>U|?|u9+w9a?#sN}#;<;i9Q-$ zjGrVt*PEPFn+4L`1|)5$I(ciz1^9>{%Ux>-1I$0s{Hn7$qF5q8XMj6+Pn;y z#ay2jjO@oqYGGC6$N^v}T1O_;Xv4HsBr~YN92kk;Qm@wY!2~9muc=lB<0+x6J2^5i znb%hAQRR(sa^-3>jQnBxj;8^qO#m>|<>YS{4CBkKw@J(T zVH}GzjXG}N5j!`f8xtuXfG(X`%$~HxNC)}XX_Y77yPNR$`3n1BdfpVjbW00Nh5VT2 z;=PBFES-hJRgYmDOMbQLtxEheEr?SMf2Rdvd7q8xy6TTuFQW{WDbmr{t)*MS?-W4! zxmS1GmsBCT?AcJim->jdU1k62k|xx!@|AvYE1otwt!zqjwnbgLjTl=zxFKEUU>83v z5!K#1>9+C5bMzyor$x(}39my}yTqCrK+;_Qv(&Ln@O*x5eb^@qzNvm`F!nM=v^1;RkyxuMtl>M=x3}J7jiea z%%9xv2<=BCnbaUbsIPLEa}bYIbx4SRa9%tG&82Vlej_PCOL~i#rbjxIGp$3<1S%mK zP2!ZoxjOWLrrKTeI1w6TMLp_N3LsA5N8HCqA=Gz2E|2sK@BbRxk1J%tAoKf#KvVZd zG^i6I=QKwey}i6{SNJJ6gg%AVrPT|eCK(GGzDpSTV*hMeyIUD*{hGSg^NJ*@?|EI< zanT&Ls4t(6xz~fb8!l5fNSs0JsX%EzB?m-b*kw0Q*ADeL_AN@TwL|?>6K-eGQ8dZ3 zEh6sA14LD2sa-hlj|Ph_j3>WoM^mEH%gA$Ch~Cj86*DgfQKfX1_|B5i#M@oLtn0>T z%qxR_)X*6*H`I-#eiud~HKxo!!yrUCy>+n@o#LNLCrQn^eTU(H0OL&n#K;(PRMc4;~}z?6Tj z^6?$M$_IF;>1pB-f(G_OH+m}2_qwloA8;H zL+WL(kf8ddVfL@%p)e?QO{{!ZJ^rzE4*`{N zHF@xw2~2Dr3FeFVrL$>&# zEgR(o8$7^zKVy9iDbOP&R}H4+4O7c(tEzoCUCk zvg_;JV_>RxJ-KOL8%)nvK5M+d7~@b(EmLY|V4UThQJ=&OF`ifDPLlR`dVBl96}ug+ zFz#y~1r|&?;J!3hTgzN8U}=9|yjHFn|G7wRw}Gt>#yPO}Xo(npMzd2|Wn=dP{7hs2 zhPpumjH9!5WvSW;pqzLa@m#VKXe?D-=OvXeouR#~YmGAuDr_(fS!??0~+$=QOF;277ZYe$v z7?%{E>){d&JmwOboFaldftTB9<7amQE$>j4&ATRGKeO6XUeCf49%(ubtN`e4yn9E} zXFiOoYYU6(gaMPRce7mhJVbPK)H^&Ag8D0GPdzink6|2EuGs&&6@A7Aq{w6ON?y!f9yH+7}qtsOf(p)akC zCb`NK+CzH95|5Fgzi6{T8*c@kfV<9CR6PK#JD<&y7jJ>Oh*l}R0DOLw$A?B0^Fe99 z@V@2@LzH9gHov9jIkXFsw|LX%qC8Cz(*TMebVyD1{^0FHWm?ZS?G5IJ^35lIrXIq7 zEf-X1aXc{{#(Gf0WA88+?aE%XE0P0-+fFS>vNHin>gzWPWT`L~>d(V(xCN-B!z!Yw zbui*?cwuT$0*t-Szd3#^4#wDb3U&biRrV=puA$4agTrq*i7 zp@VoVCFhEVnC>Th_kZTSV}%4p`dMsqcZDUM87Mw%?@z&hZYAZ^vP=Y@ADA(h`C5EF zS-vcIvJU90rG&KRs{vEuhWSq2RT$}`0;SDpE->3e*^edm;iuJJo_5c8gmHoW{+ zyJjD8)sIfbIE95$FQwoi$NR6pEc>j2aWgqaONH=E(6bB5LSL&f?&!%KOTOd50N2&5 zRy>i)kYAd5i0cbXt({n=>0*NazQ`pY?}j8$3cSq6RurF7ONcJmwPKLVuPj-G!#gEB+k&@<%3rA)3yE)b$b%fT9;&uM} zyivAhy;cRO0M#YBY<;|%j62g!m!=uXC^ydHa`|h1)Ggies;|}vk)!pFhDT^Y1ZDXz z3{Q-98~5{zOlYCrY-7Pv7(>+$Evuus?^I=SBh4{KC6Fhjr K;L^N+3H~4Bcjtuw literal 0 HcmV?d00001 diff --git a/examples/08_geo_coordinates/00_field_generation.py b/examples/08_geo_coordinates/00_field_generation.py index 91583075..6a94ad1b 100755 --- a/examples/08_geo_coordinates/00_field_generation.py +++ b/examples/08_geo_coordinates/00_field_generation.py @@ -42,7 +42,7 @@ ) ax = model.plot("vario_yadrenko", x_max=0.3) -model.fit_variogram(bin_center, emp_vario, init_guess="current", nugget=False) +model.fit_variogram(bin_center, emp_vario, nugget=False) model.plot("vario_yadrenko", ax=ax, label="fitted", x_max=0.3) ax.scatter(bin_center, emp_vario, color="k") print(model) diff --git a/examples/08_geo_coordinates/01_dwd_krige.py b/examples/08_geo_coordinates/01_dwd_krige.py index f4dfb6fe..fa7bf565 100755 --- a/examples/08_geo_coordinates/01_dwd_krige.py +++ b/examples/08_geo_coordinates/01_dwd_krige.py @@ -27,7 +27,7 @@ def get_borders_germany(): import geopandas as gp # 0.8.1 shpfile = shp_read.natural_earth("50m", "cultural", "admin_0_countries") - df = gp.read_file(shpfile) # only use the simples polygon + df = gp.read_file(shpfile) # only use the simplest polygon poly = df.loc[df["ADMIN"] == "Germany"]["geometry"].values[0][0] np.savetxt("de_borders.txt", list(poly.exterior.coords)) From c30532b67e01bd71063cbdc4064d225d66d4bc9b Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Mon, 11 Jan 2021 01:55:42 +0100 Subject: [PATCH 27/39] Tests: SRF + PyVista dim fix --- tests/test_srf.py | 1 + 1 file changed, 1 insertion(+) diff --git a/tests/test_srf.py b/tests/test_srf.py index d9606c61..efd9ceea 100644 --- a/tests/test_srf.py +++ b/tests/test_srf.py @@ -243,6 +243,7 @@ def test_calls(self): def test_mesh_pyvista(self): """Test the `.mesh` call with various PyVista meshes.""" # Create model + self.cov_model.dim = 3 srf = SRF(self.cov_model, mean=self.mean, mode_no=self.mode_no) # Get the field the normal way for comparison field = srf((self.x_tuple, self.y_tuple, self.z_tuple), seed=self.seed) From 17468730aecdc6f1f94e8b06078855eae7fce9f6 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Mon, 11 Jan 2021 02:10:54 +0100 Subject: [PATCH 28/39] CovModel: skip coverage for repr --- gstools/covmodel/tools.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/gstools/covmodel/tools.py b/gstools/covmodel/tools.py index 70b1ec6b..9c1d6c1f 100644 --- a/gstools/covmodel/tools.py +++ b/gstools/covmodel/tools.py @@ -610,7 +610,7 @@ def compare(this, that): return equal -def model_repr(model): +def model_repr(model): # pragma: no cover """ Generate the model string representation. From 3930aa4f930c73c70c4d6fd803a515d69df29d70 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Mon, 11 Jan 2021 02:11:39 +0100 Subject: [PATCH 29/39] Tests: add tests for latlon --- tests/test_latlon.py | 115 +++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 115 insertions(+) create mode 100644 tests/test_latlon.py diff --git a/tests/test_latlon.py b/tests/test_latlon.py new file mode 100644 index 00000000..6cfbcdfd --- /dev/null +++ b/tests/test_latlon.py @@ -0,0 +1,115 @@ +# -*- coding: utf-8 -*- +""" +This is the unittest of CovModel class. +""" + +import numpy as np +import unittest +import gstools as gs + + +def _rel_err(a, b): + return np.abs(a / ((a + b) / 2) - 1) + + +class TestCondition(unittest.TestCase): + def setUp(self): + self.cov_model = gs.Gaussian( + latlon=True, var=2, len_scale=777, rescale=gs.EARTH_RADIUS + ) + self.lat = self.lon = range(-80, 81) + + self.data = np.array( + [ + [52.9336, 8.237, 15.7], + [48.6159, 13.0506, 13.9], + [52.4853, 7.9126, 15.1], + [50.7446, 9.345, 17.0], + [52.9437, 12.8518, 21.9], + [53.8633, 8.1275, 11.9], + [47.8342, 10.8667, 11.4], + [51.0881, 12.9326, 17.2], + [48.406, 11.3117, 12.9], + [49.7273, 8.1164, 17.2], + [49.4691, 11.8546, 13.4], + [48.0197, 12.2925, 13.9], + [50.4237, 7.4202, 18.1], + [53.0316, 13.9908, 21.3], + [53.8412, 13.6846, 21.3], + [54.6792, 13.4343, 17.4], + [49.9694, 9.9114, 18.6], + [51.3745, 11.292, 20.2], + [47.8774, 11.3643, 12.7], + [50.5908, 12.7139, 15.8], + ] + ) + + def test_cov_model(self): + self.assertAlmostEqual( + self.cov_model.vario_yadrenko(1.234), + self.cov_model.sill - self.cov_model.cov_yadrenko(1.234), + ) + self.assertAlmostEqual( + self.cov_model.cov_yadrenko(1.234), + self.cov_model.var * self.cov_model.cor_yadrenko(1.234), + ) + + def test_vario_est(self): + + srf = gs.SRF(self.cov_model, seed=12345) + field = srf.structured((self.lat, self.lon)) + + bin_edges = [0.01 * i for i in range(30)] + bin_center, emp_vario = gs.vario_estimate( + *((self.lat, self.lon), field, bin_edges), + latlon=True, + mesh_type="structured", + sampling_size=2000, + sampling_seed=12345, + ) + mod = gs.Gaussian(latlon=True, rescale=gs.EARTH_RADIUS) + mod.fit_variogram(bin_center, emp_vario, nugget=False) + # allow 10 percent relative error + self.assertLess(_rel_err(mod.var, self.cov_model.var), 0.1) + self.assertLess(_rel_err(mod.len_scale, self.cov_model.len_scale), 0.1) + + def test_krige(self): + bin_max = np.deg2rad(8) + bin_edges = np.linspace(0, bin_max, 5) + emp_vario = gs.vario_estimate( + (self.data[:, 0], self.data[:, 1]), + self.data[:, 2], + bin_edges, + latlon=True, + ) + mod = gs.Spherical(latlon=True, rescale=gs.EARTH_RADIUS) + mod.fit_variogram(*emp_vario, nugget=False) + kri = gs.krige.Ordinary( + mod, + (self.data[:, 0], self.data[:, 1]), + self.data[:, 2], + ) + field, var = kri((self.data[:, 0], self.data[:, 1])) + for i, dat in enumerate(self.data[:, 2]): + self.assertAlmostEqual(field[i], dat) + + def test_cond_srf(self): + bin_max = np.deg2rad(8) + bin_edges = np.linspace(0, bin_max, 5) + emp_vario = gs.vario_estimate( + (self.data[:, 0], self.data[:, 1]), + self.data[:, 2], + bin_edges, + latlon=True, + ) + mod = gs.Spherical(latlon=True, rescale=gs.EARTH_RADIUS) + mod.fit_variogram(*emp_vario, nugget=False) + srf = gs.SRF(mod) + srf.set_condition((self.data[:, 0], self.data[:, 1]), self.data[:, 2]) + field = srf((self.data[:, 0], self.data[:, 1])) + for i, dat in enumerate(self.data[:, 2]): + self.assertAlmostEqual(field[i], dat) + + +if __name__ == "__main__": + unittest.main() From 120756daf1243ac1a8ee313b55ad06af54f43fe9 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Mon, 11 Jan 2021 02:23:22 +0100 Subject: [PATCH 30/39] Tests: fix reference values for modes, when sampled with mc and new rescaled length as start --- tests/test_incomprrandmeth.py | 7 ++++--- tests/test_randmeth.py | 10 ++++++---- 2 files changed, 10 insertions(+), 7 deletions(-) diff --git a/tests/test_incomprrandmeth.py b/tests/test_incomprrandmeth.py index 1fd6f7ee..a10fa9eb 100644 --- a/tests/test_incomprrandmeth.py +++ b/tests/test_incomprrandmeth.py @@ -38,9 +38,10 @@ def test_unstruct_2d(self): def test_unstruct_3d(self): modes = self.rm_3d((self.x_tuple, self.y_tuple, self.z_tuple)) - self.assertAlmostEqual(modes[0, 0], 1.49469700) - self.assertAlmostEqual(modes[0, 1], 1.38687858) - self.assertAlmostEqual(modes[1, 0], -0.27245271) + # print(modes[0, 0], modes[0, 1], modes[1, 0]) + self.assertAlmostEqual(modes[0, 0], 0.7924546333550331) + self.assertAlmostEqual(modes[0, 1], 1.660747056686244) + self.assertAlmostEqual(modes[1, 0], -0.28049855754819514) def test_assertions(self): cov_model_1d = Gaussian(dim=1, var=1.5, len_scale=2.5) diff --git a/tests/test_randmeth.py b/tests/test_randmeth.py index 422567b5..61e6299e 100644 --- a/tests/test_randmeth.py +++ b/tests/test_randmeth.py @@ -41,8 +41,9 @@ def test_unstruct_2d(self): def test_unstruct_3d(self): modes = self.rm_3d((self.x_tuple, self.y_tuple, self.z_tuple)) - self.assertAlmostEqual(modes[0], 0.55488481) - self.assertAlmostEqual(modes[1], 1.18506639) + # print(modes[:2]) + self.assertAlmostEqual(modes[0], 1.3240234883187239) + self.assertAlmostEqual(modes[1], 1.6367244277732766) def test_reset(self): modes = self.rm_2d((self.x_tuple, self.y_tuple)) @@ -61,8 +62,9 @@ def test_reset(self): self.rm_1d.model = self.cov_model_3d modes = self.rm_1d((self.x_tuple, self.y_tuple, self.z_tuple)) - self.assertAlmostEqual(modes[0], 0.55488481) - self.assertAlmostEqual(modes[1], 1.18506639) + # print(modes[:2]) + self.assertAlmostEqual(modes[0], 1.3240234883187239) + self.assertAlmostEqual(modes[1], 1.6367244277732766) self.rm_2d.mode_no = 800 modes = self.rm_2d((self.x_tuple, self.y_tuple)) From 943894ccffaf79fb5c4cd5b0531dc8c63dfa8704 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Mon, 11 Jan 2021 03:03:12 +0100 Subject: [PATCH 31/39] Variogram: fix typo in err-msg --- gstools/variogram/variogram.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/gstools/variogram/variogram.py b/gstools/variogram/variogram.py index 6f63b336..d6703d1b 100644 --- a/gstools/variogram/variogram.py +++ b/gstools/variogram/variogram.py @@ -233,7 +233,7 @@ def vario_estimate( pos, field.shape, check_stacked_shape=True ) if latlon and dim != 2: - raise ValueError("Variogram: given field to be 2D for lat-lon.") + raise ValueError("Variogram: given field needs to be 2D for lat-lon.") # prepare the field pnt_cnt = len(pos[0]) field = field.reshape((-1, pnt_cnt)) From e5af2bb7c64ebbc8d4e684eac6cd28076b39489c Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Mon, 11 Jan 2021 03:03:45 +0100 Subject: [PATCH 32/39] Tests: more latlon tests including Error testing --- tests/test_latlon.py | 46 ++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 46 insertions(+) diff --git a/tests/test_latlon.py b/tests/test_latlon.py index 6cfbcdfd..b93658c4 100644 --- a/tests/test_latlon.py +++ b/tests/test_latlon.py @@ -12,6 +12,14 @@ def _rel_err(a, b): return np.abs(a / ((a + b) / 2) - 1) +class ErrMod(gs.CovModel): + def cor(self, h): + return np.exp(-(h ** 2)) + + def fix_dim(self): + return 2 + + class TestCondition(unittest.TestCase): def setUp(self): self.cov_model = gs.Gaussian( @@ -44,6 +52,13 @@ def setUp(self): ] ) + def test_conv(self): + p_ll = gs.tools.geometric.latlon2pos((self.lat, self.lon), 2.56) + ll_p = gs.tools.geometric.pos2latlon(p_ll, 2.56) + for i, v in enumerate(self.lat): + self.assertAlmostEqual(v, ll_p[0, i]) + self.assertAlmostEqual(v, ll_p[1, i]) + def test_cov_model(self): self.assertAlmostEqual( self.cov_model.vario_yadrenko(1.234), @@ -53,6 +68,18 @@ def test_cov_model(self): self.cov_model.cov_yadrenko(1.234), self.cov_model.var * self.cov_model.cor_yadrenko(1.234), ) + self.assertAlmostEqual( + 8, self.cov_model.anisometrize(self.cov_model.isometrize((8, 6)))[0, 0] + ) + self.assertAlmostEqual( + 6, self.cov_model.anisometrize(self.cov_model.isometrize((8, 6)))[1, 0] + ) + self.assertAlmostEqual( + 1, self.cov_model.isometrize(self.cov_model.anisometrize((1, 0, 0)))[0, 0] + ) + # test if callable + self.cov_model.anis = [1, 2] + self.cov_model.angles = [1, 2, 3] def test_vario_est(self): @@ -110,6 +137,25 @@ def test_cond_srf(self): for i, dat in enumerate(self.data[:, 2]): self.assertAlmostEqual(field[i], dat) + def error_test(self): + # try fitting directional variogram + with self.assertRaises(ValueError): + mod = gs.Gaussian(latlon=True) + mod.fit_variogram([0, 1], [[0, 1], [0, 1], [0, 1]]) + # try to use fixed dim=2 with latlon + with self.assertRaises(ValueError): + ErrMod(latlon=True) + # try to estimate latlon vario on wrong dim + with self.assertRaises(ValueError): + gs.vario_estimate([[1], [1], [1]], [1], [0, 1], latlon=True) + # try to estimate directional vario with latlon + with self.assertRaises(ValueError): + gs.vario_estimate([[1], [1]], [1], [0, 1], latlon=True, angles=1) + # try to create a vector field with latlon + with self.assertRaises(ValueError): + srf = gs.SRF(mod, generator="VectorField", mode_no=2) + srf([1, 2]) + if __name__ == "__main__": unittest.main() From f6f9140bf0028e8f5c69a3265a4c73b86195ea8d Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Mon, 11 Jan 2021 03:14:43 +0100 Subject: [PATCH 33/39] Tests: latlon test restructuring --- tests/test_latlon.py | 39 ++++++++++++++++++++++----------------- 1 file changed, 22 insertions(+), 17 deletions(-) diff --git a/tests/test_latlon.py b/tests/test_latlon.py index b93658c4..595dfcef 100644 --- a/tests/test_latlon.py +++ b/tests/test_latlon.py @@ -1,6 +1,6 @@ # -*- coding: utf-8 -*- """ -This is the unittest of CovModel class. +This is the unittest for latlon related routines. """ import numpy as np @@ -22,7 +22,7 @@ def fix_dim(self): class TestCondition(unittest.TestCase): def setUp(self): - self.cov_model = gs.Gaussian( + self.cmod = gs.Gaussian( latlon=True, var=2, len_scale=777, rescale=gs.EARTH_RADIUS ) self.lat = self.lon = range(-80, 81) @@ -58,32 +58,37 @@ def test_conv(self): for i, v in enumerate(self.lat): self.assertAlmostEqual(v, ll_p[0, i]) self.assertAlmostEqual(v, ll_p[1, i]) - - def test_cov_model(self): self.assertAlmostEqual( - self.cov_model.vario_yadrenko(1.234), - self.cov_model.sill - self.cov_model.cov_yadrenko(1.234), + 8, self.cmod.anisometrize(self.cmod.isometrize((8, 6)))[0, 0] ) self.assertAlmostEqual( - self.cov_model.cov_yadrenko(1.234), - self.cov_model.var * self.cov_model.cor_yadrenko(1.234), + 6, self.cmod.anisometrize(self.cmod.isometrize((8, 6)))[1, 0] ) self.assertAlmostEqual( - 8, self.cov_model.anisometrize(self.cov_model.isometrize((8, 6)))[0, 0] + 1, self.cmod.isometrize(self.cmod.anisometrize((1, 0, 0)))[0, 0] ) + + def test_cov_model(self): self.assertAlmostEqual( - 6, self.cov_model.anisometrize(self.cov_model.isometrize((8, 6)))[1, 0] + self.cmod.vario_yadrenko(1.234), + self.cmod.sill - self.cmod.cov_yadrenko(1.234), ) self.assertAlmostEqual( - 1, self.cov_model.isometrize(self.cov_model.anisometrize((1, 0, 0)))[0, 0] + self.cmod.cov_yadrenko(1.234), + self.cmod.var * self.cmod.cor_yadrenko(1.234), ) - # test if callable - self.cov_model.anis = [1, 2] - self.cov_model.angles = [1, 2, 3] + # test if correctly handling tries to set anisotropy + self.cmod.anis = [1, 2] + self.cmod.angles = [1, 2, 3] + self.assertAlmostEqual(self.cmod.anis[0], 1) + self.assertAlmostEqual(self.cmod.anis[1], 1) + self.assertAlmostEqual(self.cmod.angles[0], 0) + self.assertAlmostEqual(self.cmod.angles[1], 0) + self.assertAlmostEqual(self.cmod.angles[2], 0) def test_vario_est(self): - srf = gs.SRF(self.cov_model, seed=12345) + srf = gs.SRF(self.cmod, seed=12345) field = srf.structured((self.lat, self.lon)) bin_edges = [0.01 * i for i in range(30)] @@ -97,8 +102,8 @@ def test_vario_est(self): mod = gs.Gaussian(latlon=True, rescale=gs.EARTH_RADIUS) mod.fit_variogram(bin_center, emp_vario, nugget=False) # allow 10 percent relative error - self.assertLess(_rel_err(mod.var, self.cov_model.var), 0.1) - self.assertLess(_rel_err(mod.len_scale, self.cov_model.len_scale), 0.1) + self.assertLess(_rel_err(mod.var, self.cmod.var), 0.1) + self.assertLess(_rel_err(mod.len_scale, self.cmod.len_scale), 0.1) def test_krige(self): bin_max = np.deg2rad(8) From 16409885c13deb108bb42aeae37caefc30044f0c Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Mon, 11 Jan 2021 13:41:59 +0100 Subject: [PATCH 34/39] Tests: latlon: minor fix --- tests/test_latlon.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_latlon.py b/tests/test_latlon.py index 595dfcef..bc973755 100644 --- a/tests/test_latlon.py +++ b/tests/test_latlon.py @@ -144,8 +144,8 @@ def test_cond_srf(self): def error_test(self): # try fitting directional variogram + mod = gs.Gaussian(latlon=True) with self.assertRaises(ValueError): - mod = gs.Gaussian(latlon=True) mod.fit_variogram([0, 1], [[0, 1], [0, 1], [0, 1]]) # try to use fixed dim=2 with latlon with self.assertRaises(ValueError): From 185c39541acb8f3cca8e3fac479a7dd75efc0a64 Mon Sep 17 00:00:00 2001 From: LSchueler Date: Tue, 12 Jan 2021 22:08:10 +0100 Subject: [PATCH 35/39] Fix some typos in examples and docstrings --- .../08_geo_coordinates/00_field_generation.py | 6 +++--- examples/08_geo_coordinates/01_dwd_krige.py | 20 +++++++++---------- examples/08_geo_coordinates/README.rst | 2 +- gstools/covmodel/base.py | 4 ++-- gstools/variogram/variogram.py | 4 ++-- 5 files changed, 18 insertions(+), 18 deletions(-) diff --git a/examples/08_geo_coordinates/00_field_generation.py b/examples/08_geo_coordinates/00_field_generation.py index 6a94ad1b..0abb1c81 100755 --- a/examples/08_geo_coordinates/00_field_generation.py +++ b/examples/08_geo_coordinates/00_field_generation.py @@ -8,10 +8,10 @@ First we setup a model, with ``latlon=True``, to get the associated Yadrenko model. -In addition we will use the earth radius provided by :any:`EARTH_RADIUS`, +In addition, we will use the earth radius provided by :any:`EARTH_RADIUS`, to have a meaningful length scale in km. -To generate the field, we simply pass ``(lat, lon)`` as position tuple +To generate the field, we simply pass ``(lat, lon)`` as the position tuple to the :any:`SRF` class. """ import gstools as gs @@ -30,7 +30,7 @@ # The :any:`vario_estimate` routine also provides a ``latlon`` switch to # indicate, that the given field is defined on geographical variables. # -# As we will see, everthing went well... Phew! +# As we will see, everthing went well... phew! bin_edges = [0.01 * i for i in range(30)] bin_center, emp_vario = gs.vario_estimate( diff --git a/examples/08_geo_coordinates/01_dwd_krige.py b/examples/08_geo_coordinates/01_dwd_krige.py index fa7bf565..5b7e690f 100755 --- a/examples/08_geo_coordinates/01_dwd_krige.py +++ b/examples/08_geo_coordinates/01_dwd_krige.py @@ -3,16 +3,16 @@ ------------------------- In this example we are going to interpolate actual temperature data from -the german weather service `DWD `_. +the German weather service `DWD `_. -Data is retrieved utilizing the beatiful package +Data is retrieved utilizing the beautiful package `wetterdienst `_, which serves as an API for the DWD data. -For better visualization, we also download a simple shapefile of the german +For better visualization, we also download a simple shapefile of the German borderline with `cartopy `_. -In order to hold the number of dependecies low, the calls of both functions +In order to keep the number of dependecies low, the calls of both functions shown beneath are commented out. """ # sphinx_gallery_thumbnail_number = 2 @@ -86,7 +86,7 @@ def get_dwd_temperature(): ############################################################################### # Now we can use this estimated variogram to fit a model to it. # Here we will use a :any:`Spherical` model. We select the ``latlon`` option -# to use the `Yadrenko` variant of the model to gain a vaild model for lat-lon +# to use the `Yadrenko` variant of the model to gain a valid model for lat-lon # coordinates and we rescale it to the earth-radius. Otherwise the length # scale would be given in radians representing the great-circle distance. # @@ -94,7 +94,7 @@ def get_dwd_temperature(): # # .. note:: # -# You need to plot the yadrenko variogram, since the standard variogram +# You need to plot the Yadrenko variogram, since the standard variogram # still holds the ordinary routine that is not respecting the great-circle # distance. @@ -109,7 +109,7 @@ def get_dwd_temperature(): # # Now we want to interpolate the data using :any:`Universal` kriging. # In order to tinker around with the data, we will use a north-south drift -# by assuming a linear correlation to the latitude. +# by assuming a linear correlation with the latitude. # This can be done as follows: north_south_drift = lambda lat, lon: lat @@ -123,8 +123,8 @@ def get_dwd_temperature(): ############################################################################### # Now we generate the kriging field, by defining a lat-lon grid that covers -# whole germany. The :any:`Krige` class provides the option to only krige the -# mean field, so one can have a glimpse at the estimated drift. +# the whole of Germany. The :any:`Krige` class provides the option to only +# krige the mean field, so one can have a glimpse at the estimated drift. g_lat = np.arange(47, 56.1, 0.1) g_lon = np.arange(5, 16.1, 0.1) @@ -155,7 +155,7 @@ def get_dwd_temperature(): fig.colorbar(co2, ax=ax, **fmt).set_label("T in [°C]") ############################################################################### -# To get a better impression of the estimated north-south drift, well take +# To get a better impression of the estimated north-south drift, we'll take # a look at a cross-section at a longitude of 10 degree: fig, ax = plt.subplots() diff --git a/examples/08_geo_coordinates/README.rst b/examples/08_geo_coordinates/README.rst index 0fdaeb87..4196f3c9 100644 --- a/examples/08_geo_coordinates/README.rst +++ b/examples/08_geo_coordinates/README.rst @@ -9,7 +9,7 @@ given by: - longitude ``lon``: specifies the east–west position of a point on the Earth's surface If you want to use this feature for field generation or Kriging, you -have to setup a geographical covariance Model by setting ``latlon=True`` +have to set up a geographical covariance Model by setting ``latlon=True`` in your desired model (see :any:`CovModel`): .. code-block:: python diff --git a/gstools/covmodel/base.py b/gstools/covmodel/base.py index 7dd0c9f8..4eb8c85d 100644 --- a/gstools/covmodel/base.py +++ b/gstools/covmodel/base.py @@ -96,11 +96,11 @@ class CovModel(metaclass=InitSubclassMeta): Whether the model is describing 2D fields on earths surface described by latitude and longitude. When using this, the model will internally use the associated 'Yadrenko' model to represent a valid model. - This means, the spatial distance :math:`r` will be replace with + This means, the spatial distance :math:`r` will be replaced by :math:`2\sin(\alpha/2)`, where :math:`\alpha` is the great-circle distance, which is equal to the spatial distance of two points in 3D. As a consequence, `dim` will be set to `3` and anisotropy will be - disabled. `rescale` can be set to earths radius, + disabled. `rescale` can be set to e.g. earth's radius, to have a meaningful `len_scale` parameter. Default: False var_raw : :class:`float` or :any:`None`, optional diff --git a/gstools/variogram/variogram.py b/gstools/variogram/variogram.py index d6703d1b..ffef99b6 100644 --- a/gstools/variogram/variogram.py +++ b/gstools/variogram/variogram.py @@ -146,8 +146,8 @@ def vario_estimate( Whether the data is representing 2D fields on earths surface described by latitude and longitude. When using this, the estimator will use great-circle distance for variogram estimation. - Note, that only an isotropic variogram can be estimated and an - ValueError will be raised, if direction were specified. + Note, that only an isotropic variogram can be estimated and a + ValueError will be raised, if a direction was specified. Bin edges need to be given in radians in this case. Default: False direction : :class:`list` of :class:`numpy.ndarray`, optional From f67aaae7912793136426cba9bd11934ae9bf63ca Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Wed, 13 Jan 2021 00:28:05 +0100 Subject: [PATCH 36/39] Examples: latlon - cleaner formatting --- examples/08_geo_coordinates/00_field_generation.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/examples/08_geo_coordinates/00_field_generation.py b/examples/08_geo_coordinates/00_field_generation.py index 0abb1c81..a44c7735 100755 --- a/examples/08_geo_coordinates/00_field_generation.py +++ b/examples/08_geo_coordinates/00_field_generation.py @@ -34,7 +34,9 @@ bin_edges = [0.01 * i for i in range(30)] bin_center, emp_vario = gs.vario_estimate( - *((lat, lon), field, bin_edges), + (lat, lon), + field, + bin_edges, latlon=True, mesh_type="structured", sampling_size=2000, From 717558552ad1ae5095cf0b8a6499941347801197 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Wed, 13 Jan 2021 00:29:31 +0100 Subject: [PATCH 37/39] Examples: latlon - use function instead of lambda definition --- examples/08_geo_coordinates/01_dwd_krige.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/examples/08_geo_coordinates/01_dwd_krige.py b/examples/08_geo_coordinates/01_dwd_krige.py index 5b7e690f..8b3aacd5 100755 --- a/examples/08_geo_coordinates/01_dwd_krige.py +++ b/examples/08_geo_coordinates/01_dwd_krige.py @@ -112,7 +112,10 @@ def get_dwd_temperature(): # by assuming a linear correlation with the latitude. # This can be done as follows: -north_south_drift = lambda lat, lon: lat + +def north_south_drift(lat, lon): + return lat + uk = gs.krige.Universal( model=model, From 7a7bfea6e86c1887240949a3715b363ba70a33c3 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Wed, 13 Jan 2021 00:31:40 +0100 Subject: [PATCH 38/39] CovModel: reference latlon in yadrenko routines --- gstools/covmodel/base.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/gstools/covmodel/base.py b/gstools/covmodel/base.py index 4eb8c85d..56312163 100644 --- a/gstools/covmodel/base.py +++ b/gstools/covmodel/base.py @@ -249,15 +249,15 @@ def cor_axis(self, r, axis=0): return self.correlation(np.abs(r) / self.anis[axis - 1]) def vario_yadrenko(self, zeta): - r"""Yadrenko variogram for great-circle distance from geo-coords.""" + r"""Yadrenko variogram for great-circle distance from latlon-pos.""" return self.variogram(2 * np.sin(zeta / 2)) def cov_yadrenko(self, zeta): - r"""Yadrenko covariance for great-circle distance from geo-coords.""" + r"""Yadrenko covariance for great-circle distance from latlon-pos.""" return self.covariance(2 * np.sin(zeta / 2)) def cor_yadrenko(self, zeta): - r"""Yadrenko correlation for great-circle distance from geo-coords.""" + r"""Yadrenko correlation for great-circle distance from latlon-pos.""" return self.correlation(2 * np.sin(zeta / 2)) def vario_spatial(self, pos): From 606c4ffb06dfc1034276a6d2a7f47ff5e75cd4e9 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Wed, 13 Jan 2021 00:34:05 +0100 Subject: [PATCH 39/39] Tests: remove test commments --- tests/test_incomprrandmeth.py | 1 - tests/test_latlon.py | 1 - tests/test_randmeth.py | 2 -- 3 files changed, 4 deletions(-) diff --git a/tests/test_incomprrandmeth.py b/tests/test_incomprrandmeth.py index a10fa9eb..8d1a67bc 100644 --- a/tests/test_incomprrandmeth.py +++ b/tests/test_incomprrandmeth.py @@ -38,7 +38,6 @@ def test_unstruct_2d(self): def test_unstruct_3d(self): modes = self.rm_3d((self.x_tuple, self.y_tuple, self.z_tuple)) - # print(modes[0, 0], modes[0, 1], modes[1, 0]) self.assertAlmostEqual(modes[0, 0], 0.7924546333550331) self.assertAlmostEqual(modes[0, 1], 1.660747056686244) self.assertAlmostEqual(modes[1, 0], -0.28049855754819514) diff --git a/tests/test_latlon.py b/tests/test_latlon.py index bc973755..7ae44975 100644 --- a/tests/test_latlon.py +++ b/tests/test_latlon.py @@ -87,7 +87,6 @@ def test_cov_model(self): self.assertAlmostEqual(self.cmod.angles[2], 0) def test_vario_est(self): - srf = gs.SRF(self.cmod, seed=12345) field = srf.structured((self.lat, self.lon)) diff --git a/tests/test_randmeth.py b/tests/test_randmeth.py index 61e6299e..54dd0aa3 100644 --- a/tests/test_randmeth.py +++ b/tests/test_randmeth.py @@ -41,7 +41,6 @@ def test_unstruct_2d(self): def test_unstruct_3d(self): modes = self.rm_3d((self.x_tuple, self.y_tuple, self.z_tuple)) - # print(modes[:2]) self.assertAlmostEqual(modes[0], 1.3240234883187239) self.assertAlmostEqual(modes[1], 1.6367244277732766) @@ -62,7 +61,6 @@ def test_reset(self): self.rm_1d.model = self.cov_model_3d modes = self.rm_1d((self.x_tuple, self.y_tuple, self.z_tuple)) - # print(modes[:2]) self.assertAlmostEqual(modes[0], 1.3240234883187239) self.assertAlmostEqual(modes[1], 1.6367244277732766)