From 3e0758a73089e2561edf493613231f1ecf6944f5 Mon Sep 17 00:00:00 2001 From: LSchueler Date: Mon, 27 Jul 2020 17:30:52 +0200 Subject: [PATCH 1/4] Add great-circle distance fct. to variogram estim. Now, variogram estimation can be performed for geo-coordinates. --- gstools/variogram/estimator.pyx | 28 ++++++++++++++++++++++++++-- gstools/variogram/variogram.py | 29 ++++++++++++++++++++++++++++- 2 files changed, 54 insertions(+), 3 deletions(-) diff --git a/gstools/variogram/estimator.pyx b/gstools/variogram/estimator.pyx index 41414f6c..65417e4e 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 +from libc.math cimport fabs, sqrt, pow, sin, cos, atan2, M_PI cimport numpy as np @@ -36,6 +36,27 @@ cdef inline double _distance_2d( ) nogil: return sqrt((x[i] - x[j]) * (x[i] - x[j]) + (y[i] - y[j]) * (y[i] - y[j])) +cdef inline double _distance_2d_haversine( + const double[:] x, + const double[:] y, + const double[:] z, + const int i, + const int j +) nogil: + # x is latitude in degree + # y is longitude in degree + cdef double deg_2_rad = M_PI / 180.0 + cdef double diff_lat = (x[j] - x[i]) * deg_2_rad + cdef double diff_lon = (y[j] - y[i]) * deg_2_rad + cdef double arg = ( + pow(sin(diff_lat/2.0), 2) + + cos(x[i]*deg_2_rad) * + cos(x[j]*deg_2_rad) * + pow(sin(diff_lon/2.0), 2) + ) + # earth radius for WGS84 ellipsoid r ~ 6371km + return 12742000.0 * atan2(sqrt(arg), sqrt(1.0-arg)) + cdef inline double _distance_3d( const double[:] x, const double[:] y, @@ -117,7 +138,8 @@ def unstructured( const double[:] x, const double[:] y=None, const double[:] z=None, - str estimator_type='m' + str estimator_type='m', + str distance_type='e' ): if x.shape[0] != f.shape[0]: raise ValueError('len(x) = {0} != len(f) = {1} '. @@ -138,6 +160,8 @@ def unstructured( raise ValueError('len(y) = {0} != len(f) = {1} '. format(y.shape[0], f.shape[0])) distance = _distance_2d + if distance_type == 's': + distance = _distance_2d_haversine # 1d else: distance = _distance_1d diff --git a/gstools/variogram/variogram.py b/gstools/variogram/variogram.py index e13ddf32..748559a3 100644 --- a/gstools/variogram/variogram.py +++ b/gstools/variogram/variogram.py @@ -32,6 +32,18 @@ def _set_estimator(estimator): ) return cython_estimator +def _set_distance(distance): + """Translate the verbose Python distance identifier to single char.""" + if distance.lower() == "euclidean": + cython_distance = "e" + elif distance.lower() == "spherical": + cython_distance = "s" + else: + raise ValueError( + "Unknown variogram distance function " + str(distance) + ) + return cython_distance + def vario_estimate_unstructured( pos, @@ -40,6 +52,7 @@ def vario_estimate_unstructured( sampling_size=None, sampling_seed=None, estimator="matheron", + distance="euclidean", ): r""" Estimates the variogram on a unstructured grid. @@ -92,6 +105,13 @@ def vario_estimate_unstructured( * "cressie": an estimator more robust to outliers Default: "matheron" + distance : :class:`str`, optional + the distance function, possible choices: + + * "euclidean": the distance in Euclidean spaces + * "spherical": the great-circle distance on spheres + + Default: "euclidean" Returns ------- @@ -116,11 +136,18 @@ def vario_estimate_unstructured( z = z[sampled_idx] cython_estimator = _set_estimator(estimator) + cython_distance = _set_distance(distance) return ( bin_centres, unstructured( - field, bin_edges, x, y, z, estimator_type=cython_estimator + field, + bin_edges, + x, + y, + z, + estimator_type=cython_estimator, + distance_type=cython_distance, ), ) From 34a1b6d7197dda68e5ae040a863a9b5d09825f57 Mon Sep 17 00:00:00 2001 From: LSchueler Date: Mon, 27 Jul 2020 17:32:02 +0200 Subject: [PATCH 2/4] Add simple test for great-circle dist. fct. --- tests/test_variogram_unstructured.py | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/tests/test_variogram_unstructured.py b/tests/test_variogram_unstructured.py index c5a59a85..66c90283 100644 --- a/tests/test_variogram_unstructured.py +++ b/tests/test_variogram_unstructured.py @@ -183,6 +183,15 @@ def test_sampling_3d(self): self.assertAlmostEqual(gamma[len(gamma) // 2], var, places=2) self.assertAlmostEqual(gamma[-1], var, places=2) + def test_spherical_dist(self): + lat = np.array((38.725267, 0.313611, 31.205167)) + lon = np.array((-9.150019, 32.581111, -7.862263)) + d = np.array((2, 1200, 2570)) + bins = np.array((0.0, 844_999.0, 6_021_999.0, 10_000_000.0)) + _, gamma = vario_estimate_unstructured((lat, lon), d, bins, distance='spherical') + self.assertAlmostEqual(gamma[0], 3297312.) + self.assertAlmostEqual(gamma[1], 828026.) + def test_assertions(self): x = np.arange(0, 10) x_e = np.arange(0, 11) From 8b1cdac937bc22e56c461f91885da4d69d2db05d Mon Sep 17 00:00:00 2001 From: LSchueler Date: Mon, 27 Jul 2020 17:55:40 +0200 Subject: [PATCH 3/4] Support Py35 --- tests/test_variogram_unstructured.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_variogram_unstructured.py b/tests/test_variogram_unstructured.py index 66c90283..5cfe623d 100644 --- a/tests/test_variogram_unstructured.py +++ b/tests/test_variogram_unstructured.py @@ -187,7 +187,7 @@ def test_spherical_dist(self): lat = np.array((38.725267, 0.313611, 31.205167)) lon = np.array((-9.150019, 32.581111, -7.862263)) d = np.array((2, 1200, 2570)) - bins = np.array((0.0, 844_999.0, 6_021_999.0, 10_000_000.0)) + bins = np.array((0.0, 844999.0, 6021999.0, 10000000.0)) _, gamma = vario_estimate_unstructured((lat, lon), d, bins, distance='spherical') self.assertAlmostEqual(gamma[0], 3297312.) self.assertAlmostEqual(gamma[1], 828026.) From ee4ce648bcfbdd603e22316f4805ec3a57bc8301 Mon Sep 17 00:00:00 2001 From: LSchueler Date: Thu, 20 Aug 2020 10:41:59 +0200 Subject: [PATCH 4/4] Change great-circle distance name --- gstools/variogram/variogram.py | 4 ++-- tests/test_variogram_unstructured.py | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/gstools/variogram/variogram.py b/gstools/variogram/variogram.py index 748559a3..b3ec89e0 100644 --- a/gstools/variogram/variogram.py +++ b/gstools/variogram/variogram.py @@ -36,7 +36,7 @@ def _set_distance(distance): """Translate the verbose Python distance identifier to single char.""" if distance.lower() == "euclidean": cython_distance = "e" - elif distance.lower() == "spherical": + elif distance.lower() == "great-circle": cython_distance = "s" else: raise ValueError( @@ -109,7 +109,7 @@ def vario_estimate_unstructured( the distance function, possible choices: * "euclidean": the distance in Euclidean spaces - * "spherical": the great-circle distance on spheres + * "great-circle": the great-circle distance on spheres Default: "euclidean" diff --git a/tests/test_variogram_unstructured.py b/tests/test_variogram_unstructured.py index 5cfe623d..8faa998e 100644 --- a/tests/test_variogram_unstructured.py +++ b/tests/test_variogram_unstructured.py @@ -188,7 +188,7 @@ def test_spherical_dist(self): lon = np.array((-9.150019, 32.581111, -7.862263)) d = np.array((2, 1200, 2570)) bins = np.array((0.0, 844999.0, 6021999.0, 10000000.0)) - _, gamma = vario_estimate_unstructured((lat, lon), d, bins, distance='spherical') + _, gamma = vario_estimate_unstructured((lat, lon), d, bins, distance='great-circle') self.assertAlmostEqual(gamma[0], 3297312.) self.assertAlmostEqual(gamma[1], 828026.)