Merge pull request #430 from DataverseLabs/atp-variance-error

Updated variance error calculations in area-to-point poisson kriging
DataverseLabs · Oct 26, 2024 · e4a5ff3 · e4a5ff3
2 parents fd8c757 + f08e893
commit e4a5ff3
Show file tree

Hide file tree

Showing 8 changed files with 78 additions and 24 deletions.
diff --git a/README.md b/README.md
@@ -4,7 +4,7 @@
 
 # Pyinterpolate
 
-**version 0.5.2** - *Mykolaiv*
+**version 0.5.3** - *Mykolaiv*
 
 ---
 

diff --git a/changelog.rst b/changelog.rst
@@ -6,6 +6,12 @@ Pyinterpolate is the Python library for **geostatistics** and **spatial statisti
 Changes by date
 ===============
 
+2024-10-26
+----------
+**version 0.5.3**
+
+* (logic) debugged variance error calculations for Area-to-Point Poisson Kriging
+
 2024-06-26
 ----------
 **version 0.5.2**

diff --git a/docs/source/conf.py b/docs/source/conf.py
@@ -15,7 +15,7 @@
 project = 'Pyinterpolate'
 copyright = '2024, Szymon Moliński'
 author = 'Szymon Moliński'
-release = '0.5.2'
+release = '0.5.3'
 
 # -- General configuration ---------------------------------------------------
 # https://www.sphinx-doc.org/en/master/usage/configuration.html#general-configuration

diff --git a/docs/source/index.rst b/docs/source/index.rst
@@ -6,7 +6,7 @@
 Pyinterpolate
 =============
 
-**version 0.5.2** - *Mykolaiv*
+**version 0.5.3** - *Mykolaiv*
 ------------------------------------
 
 .. image:: imgs/pyinterpolate-banner.png

diff --git a/pyinterpolate/__init__.py b/pyinterpolate/__init__.py
@@ -46,4 +46,4 @@
 from pyinterpolate.viz import interpolate_raster
 
 
-__version__ = "0.5.1"
+__version__ = "0.5.3"
diff --git a/pyinterpolate/kriging/models/block/area_to_point_poisson_kriging.py b/pyinterpolate/kriging/models/block/area_to_point_poisson_kriging.py
@@ -18,7 +18,7 @@
 from pyinterpolate.kriging.utils.kwarnings import ExperimentalFeatureWarning
 from pyinterpolate.processing.preprocessing.blocks import Blocks, PointSupport
 from pyinterpolate.processing.select_values import select_poisson_kriging_data, prepare_pk_known_areas, \
-    get_aggregated_point_support_values
+    get_aggregated_point_support_values, get_distances_within_unknown
 from pyinterpolate.processing.transform.transform import transform_ps_to_dict, get_areal_values_from_agg, sem_to_cov
 from pyinterpolate.variogram import TheoreticalVariogram
 
@@ -194,6 +194,18 @@ def area_to_point_pk(semivariogram_model: TheoreticalVariogram,
     transposed = avg_b2p_covariances.T
     logging.info("POISSON KRIGING: AREA-TO-POINT | Kriging system weights has been prepared")
 
+    if isinstance(unknown_block[0], np.ndarray):
+        u_idx = unknown_block[0][0]
+    else:
+        u_idx = unknown_block[0]
+
+    distances_within_unknown_block = get_distances_within_unknown(unknown_block_point_support)
+    semivariance_within_unknown = b2b_semivariance.calculate_average_semivariance({
+        u_idx: distances_within_unknown_block
+    })[u_idx]
+
+    covariance_within_unknown = sem_to_cov([semivariance_within_unknown], sill)[0]
+
     predicted_points = []
 
     for idx, point in enumerate(transposed):
@@ -225,6 +237,8 @@ def area_to_point_pk(semivariogram_model: TheoreticalVariogram,
 
         # Calculate error
         sigmasq = np.matmul(w.T, point)
+        sigmasq = covariance_within_unknown - sigmasq
+
         if sigmasq < 0:
             if raise_when_negative_error:
                 logging.info(f"POISSON KRIGING: AREA-TO-POINT | Negative variance error for point {point}")

diff --git a/setup.cfg b/setup.cfg
@@ -3,7 +3,7 @@ name = pyinterpolate
 description = Spatial Interpolation in Python
 long_description = file: README.md
 long_description_content_type = text/markdown; charset=UTF-8
-version = 0.5.2
+version = 0.5.3
 url = https://github.com/DataverseLabs/pyinterpolate
 download_url = https://github.com/DataverseLabs/pyinterpolate/archive/
 author = Szymon Moliński

diff --git a/tutorials/C15-Poisson-Kriging-Area-to-Point.ipynb b/tutorials/C15-Poisson-Kriging-Area-to-Point.ipynb