diff --git a/README.md b/README.md index a37d8c70..13613f20 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/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/index.rst b/docs/source/index.rst index f76f299f..296de012 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 `__ +- `Geographic Coordinates `__ - `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 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 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) diff --git a/examples/00_misc/README.rst b/examples/00_misc/README.rst index 61122995..bef7ae57 100644 --- a/examples/00_misc/README.rst +++ b/examples/00_misc/README.rst @@ -5,9 +5,5 @@ More examples which do not really fit into other categories. Some are not more than a code snippet, while others are more complex and more than one part of GSTools is involved. -.. only:: html - - Gallery - ------- - - Below is a gallery of examples +Examples +-------- diff --git a/examples/01_random_field/04_srf_merge.py b/examples/01_random_field/04_srf_merge.py 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/01_random_field/README.rst b/examples/01_random_field/README.rst index bb90ead2..6b226b2f 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/01_random_field/mesh_ensemble.vtk b/examples/01_random_field/mesh_ensemble.vtk new file mode 100644 index 00000000..0ea0d7bc Binary files /dev/null and b/examples/01_random_field/mesh_ensemble.vtk differ diff --git a/examples/02_cov_model/README.rst b/examples/02_cov_model/README.rst index 3bde4b84..6f8d3dad 100644 --- a/examples/02_cov_model/README.rst +++ b/examples/02_cov_model/README.rst @@ -71,6 +71,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 @@ -78,11 +79,9 @@ As a special feature, we also provide truncated power law (TPL) covariance model TPLGaussian TPLExponential TPLStable - TPLSimple - -.. only:: html - Gallery - ------- +These models provide a lower and upper length scale truncation +for superpositioned models. - Below is a gallery of examples +Examples +-------- diff --git a/examples/03_variogram/README.rst b/examples/03_variogram/README.rst index f242b71c..8eb42a8a 100644 --- a/examples/03_variogram/README.rst +++ b/examples/03_variogram/README.rst @@ -10,9 +10,5 @@ The same `(semi-)variogram `_. + +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 +borderline with `cartopy `_. + +In order to keep 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 simplest 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 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. +# +# 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 with the latitude. +# This can be done as follows: + + +def north_south_drift(lat, lon): + return 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 +# 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) + +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, we'll 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/README.rst b/examples/08_geo_coordinates/README.rst new file mode 100644 index 00000000..4196f3c9 --- /dev/null +++ b/examples/08_geo_coordinates/README.rst @@ -0,0 +1,67 @@ +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 set up 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 `field_dim` will be 2 in this case + (see :any:`CovModel.field_dim`). + +Examples +-------- 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 diff --git a/gstools/__init__.py b/gstools/__init__.py index a23a8dfd..695dd92c 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 @@ -68,7 +66,6 @@ Truncated Power Law Covariance Models ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -.. currentmodule:: gstools.covmodel.tpl_models .. autosummary:: TPLGaussian TPLExponential @@ -104,12 +101,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, @@ -182,6 +189,7 @@ __all__ += [ "SRF", "rotated_main_axes", + "EARTH_RADIUS", "vtk_export", "vtk_export_structured", "vtk_export_unstructured", diff --git a/gstools/covmodel/base.py b/gstools/covmodel/base.py index 56b5b924..56312163 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 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 e.g. earth's 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 latlon-pos.""" + return self.variogram(2 * np.sin(zeta / 2)) + + def cov_yadrenko(self, zeta): + 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 latlon-pos.""" + 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,16 @@ 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)) + 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) 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 +732,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 +750,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 +760,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 +770,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 +780,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 +791,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 +801,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 +811,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 +821,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 +832,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 +844,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 +858,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 field_dim(self): + """:class:`int`: The field dimension of the model.""" + return 2 if self.latlon else self.dim + # standard parameters @property @@ -900,9 +921,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 +945,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 +960,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 +1129,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 +1141,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/fit.py b/gstools/covmodel/fit.py index 32056063..0fa17075 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 @@ -321,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, ) 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 diff --git a/gstools/covmodel/tools.py b/gstools/covmodel/tools.py index 96a5b0ae..9c1d6c1f 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): # pragma: no cover + """ + 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 diff --git a/gstools/field/base.py b/gstools/field/base.py index 487a2b64..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: @@ -209,18 +209,21 @@ def pre_pos(self, pos, mesh_type="unstructured"): """ # save mesh-type self.mesh_type = mesh_type + dim = self.model.field_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 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 diff --git a/gstools/field/plot.py b/gstools/field/plot.py index 29a77f81..edfbd6ee 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.field_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.field_dim == 2: + ax = _plot_2d( + fld.pos, plot_field, fld.mesh_type, fig, ax, fld.model.latlon + ) + 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!") @@ -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." 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 diff --git a/gstools/krige/base.py b/gstools/krige/base.py index c468bcb6..591e38e9 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) @@ -332,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 @@ -403,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 @@ -436,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): 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) 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..ffef99b6 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 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 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 needs 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,15 @@ def vario_estimate( cython_estimator = _set_estimator(estimator) # run if dir_no == 0: + # "h"aversine or "e"uclidean distance type + distance_type = "h" if latlon else "e" estimates, counts = unstructured( dim, field, bin_edges, pos, estimator_type=cython_estimator, + distance_type=distance_type, ) else: estimates, counts = directional( diff --git a/tests/test_incomprrandmeth.py b/tests/test_incomprrandmeth.py index 1fd6f7ee..8d1a67bc 100644 --- a/tests/test_incomprrandmeth.py +++ b/tests/test_incomprrandmeth.py @@ -38,9 +38,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], 1.49469700) - self.assertAlmostEqual(modes[0, 1], 1.38687858) - self.assertAlmostEqual(modes[1, 0], -0.27245271) + 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_latlon.py b/tests/test_latlon.py new file mode 100644 index 00000000..7ae44975 --- /dev/null +++ b/tests/test_latlon.py @@ -0,0 +1,165 @@ +# -*- coding: utf-8 -*- +""" +This is the unittest for latlon related routines. +""" + +import numpy as np +import unittest +import gstools as gs + + +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.cmod = 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_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]) + self.assertAlmostEqual( + 8, self.cmod.anisometrize(self.cmod.isometrize((8, 6)))[0, 0] + ) + self.assertAlmostEqual( + 6, self.cmod.anisometrize(self.cmod.isometrize((8, 6)))[1, 0] + ) + self.assertAlmostEqual( + 1, self.cmod.isometrize(self.cmod.anisometrize((1, 0, 0)))[0, 0] + ) + + def test_cov_model(self): + self.assertAlmostEqual( + self.cmod.vario_yadrenko(1.234), + self.cmod.sill - self.cmod.cov_yadrenko(1.234), + ) + self.assertAlmostEqual( + self.cmod.cov_yadrenko(1.234), + self.cmod.var * self.cmod.cor_yadrenko(1.234), + ) + # 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.cmod, 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.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) + 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) + + def error_test(self): + # try fitting directional variogram + mod = gs.Gaussian(latlon=True) + with self.assertRaises(ValueError): + 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() diff --git a/tests/test_randmeth.py b/tests/test_randmeth.py index 422567b5..54dd0aa3 100644 --- a/tests/test_randmeth.py +++ b/tests/test_randmeth.py @@ -41,8 +41,8 @@ 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) + 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 +61,8 @@ 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) + 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)) 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)