Skip to content

Commit

Permalink
updated semivariogram class models
Browse files Browse the repository at this point in the history
  • Loading branch information
SimonMolinsky committed Oct 17, 2021
1 parent 5378b14 commit 4b9bd65
Show file tree
Hide file tree
Showing 15 changed files with 224 additions and 205 deletions.
1 change: 1 addition & 0 deletions changelog.rst
Original file line number Diff line number Diff line change
Expand Up @@ -13,6 +13,7 @@ Changes by date

* `self.points_values` chenged to `self.points_array` in `TheoreticalSemivariogram` class,
* `NaN` values are tested and checked in `calc_semivariance_from_pt_cloud()` function,
* new semivariogram models included in the package: **cubic**, **circular**, **power**,


2021-08-23
Expand Down
5 changes: 4 additions & 1 deletion docs/source/code_documentation/semivariance.ipynb
Original file line number Diff line number Diff line change
Expand Up @@ -707,6 +707,9 @@
" - gaussian_model(distance, nugget, sill, semivar_range)\n",
" - exponential_model(distance, nugget, sill, semivar_range)\n",
" - linear_model(distance, nugget, sill, semivar_range)\n",
" - cubic_model(distance, nugget, sill, semivar_range)\n",
" - circular_model(distance, nugget, sill, semivar_range)\n",
" - power_model(distance, nugget, sill, semivar_range)\n",
"\n",
"\n",
"INITIALIZATION PARAMS:\n",
Expand Down Expand Up @@ -851,4 +854,4 @@
},
"nbformat": 4,
"nbformat_minor": 4
}
}
29 changes: 16 additions & 13 deletions pyinterpolate/semivariance/semivariogram_fit/fit_semivariance.py
Original file line number Diff line number Diff line change
@@ -1,6 +1,7 @@
"""
Authors:
Scott Gallacher | @scottgallacher-3
Szymon Molinski | @szymon-datalions
Contributors:
Expand Down Expand Up @@ -249,7 +250,6 @@ def cubic_model(lags, nugget, sill, semivar_range):

if lags[0] == 0:
gamma[0] = 0


return gamma

Expand All @@ -269,8 +269,10 @@ def circular_model(lags, nugget, sill, semivar_range):
a = lag / range
(Model 1 Source: https://desktop.arcgis.com/en/arcmap/10.3/tools/spatial-analyst-toolbox/how-kriging-works.htm#GUID-94A34A70-DBCF-4B23-A198-BB50FB955DC0)
(Model 2 Source: https://pro.arcgis.com/en/pro-app/latest/tool-reference/spatial-analyst/how-kriging-works.htm#GUID-94A34A70-DBCF-4B23-A198-BB50FB955DC0)
(Model 1 Source:
https://desktop.arcgis.com/en/arcmap/10.3/tools/spatial-analyst-toolbox/how-kriging-works.htm#GUID-94A34A70-DBCF-4B23-A198-BB50FB955DC0)
(Model 2 Source:
https://pro.arcgis.com/en/pro-app/latest/tool-reference/spatial-analyst/how-kriging-works.htm#GUID-94A34A70-DBCF-4B23-A198-BB50FB955DC0)
INPUT:
Expand All @@ -283,18 +285,19 @@ def circular_model(lags, nugget, sill, semivar_range):
:return: an array of modeled values for given range. Values are calculated based on the circular model.
"""
#check conditions:
#apparently, even using np.where uncovers invalid values in the arccos and square root
#but as long as lag <= range this shouldn't happen
#use np.clip on the arrays to be passed
# TODO: check conditions:
# apparently, even using np.where uncovers invalid values in the arccos and square root
# but as long as lag <= range this shouldn't happen
# use np.clip on the arrays to be passed
a = lags / semivar_range

#use np.clip to limit range of values passed into np.arccos and np.sqrt
gamma = np.where((lags <= semivar_range),
(nugget + sill*(1 - 2/np.pi * np.arccos(np.clip(a, -1, 1)) * np.sqrt(1 - np.clip(a**2, -1, 1))) ),
(nugget + sill))
# use np.clip to limit range of values passed into np.arccos and np.sqrt
# gamma = np.where((lags <= semivar_range),
# (nugget + sill*(1 - 2/np.pi * np.arccos(np.clip(a, -1, 1)) *
# np.sqrt(1 - np.clip(a**2, -1, 1))) ),
# (nugget + sill))

#second formula found which seems to fit better, and looks as expected
# second formula found which seems to fit better, and looks as expected
gamma = nugget + (2/np.pi) * sill*(a * np.sqrt(1 - np.clip(a**2, -1, 1)) + np.arcsin(np.clip(a, -1, 1)))

if lags[0] == 0:
Expand Down Expand Up @@ -322,7 +325,7 @@ def fit_semivariance(self, model_type, number_of_ranges=16):
'spherical': self.spherical_model,
'exponential': self.exponential_model,
'linear': self.linear_model,
'gaussian': self.gaussian_model,
'gaussian': self.gaussian_model
}
model = models[model_type]
self.chosen_model_name = model_type
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -105,7 +105,7 @@
"outputs": [
{
"data": {
"text/plain": "array([25019,\n <shapely.geometry.multipolygon.MultiPolygon object at 0x7fcbbb27df90>,\n 2132629.599151873, 557971.1559491967, 192.2], dtype=object)"
"text/plain": "array([25019,\n <shapely.geometry.multipolygon.MultiPolygon object at 0x7fede1a2d310>,\n 2132629.599151873, 557971.1559491967, 192.2], dtype=object)"
},
"execution_count": 2,
"metadata": {},
Expand Down Expand Up @@ -448,11 +448,11 @@
"name": "stderr",
"output_type": "stream",
"text": [
"/Users/szymonmolinski/miniconda3/envs/pyintpk/lib/python3.7/site-packages/numpy/core/fromnumeric.py:3420: RuntimeWarning: Mean of empty slice.\n",
"/home/szymon/miniconda3/envs/pyintpk/lib/python3.7/site-packages/numpy/core/fromnumeric.py:3420: RuntimeWarning: Mean of empty slice.\n",
" out=out, **kwargs)\n",
"/Users/szymonmolinski/miniconda3/envs/pyintpk/lib/python3.7/site-packages/numpy/core/_methods.py:190: RuntimeWarning: invalid value encountered in long_scalars\n",
"/home/szymon/miniconda3/envs/pyintpk/lib/python3.7/site-packages/numpy/core/_methods.py:190: RuntimeWarning: invalid value encountered in long_scalars\n",
" ret = ret / rcount\n",
"/Users/szymonmolinski/Documents/OpenSource/pyinterpolate-environment/pyinterpolate/pyinterpolate/semivariance/semivariogram_fit/fit_semivariance.py:313: UserWarning: WARNING: linear model fitted to the experimental variogram is better than the core models!\n",
"/home/szymon/Documents/a01_repos/pyinterpolate-environment/pyinterpolate/pyinterpolate/semivariance/semivariogram_fit/fit_semivariance.py:426: UserWarning: WARNING: linear model fitted to the experimental variogram is better than the core models!\n",
" warnings.warn(warning_msg)\n"
]
}
Expand Down Expand Up @@ -627,7 +627,15 @@
"execution_count": 12,
"id": "04d61d71",
"metadata": {},
"outputs": [],
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"CPLE_NotSupported in driver ESRI Shapefile does not support creation option ENCODING\n"
]
}
],
"source": [
"# Save interpolation results\n",
"\n",
Expand Down

Large diffs are not rendered by default.

42 changes: 21 additions & 21 deletions tutorials/Ordinary and Simple Kriging (Basic).ipynb

Large diffs are not rendered by default.

Loading

0 comments on commit 4b9bd65

Please sign in to comment.