{"id":1195,"date":"2026-10-19T07:05:17","date_gmt":"2026-10-19T07:05:17","guid":{"rendered":"https:\/\/helpwithmygisassignment.com\/blog\/?p=1195"},"modified":"2026-10-19T07:05:34","modified_gmt":"2026-10-19T07:05:34","slug":"spatial-interpolation-in-python","status":"publish","type":"post","link":"https:\/\/helpwithmygisassignment.com\/blog\/spatial-interpolation-in-python\/","title":{"rendered":"Spatial Interpolation in Python: Kriging and IDW using PyKrige and SciPy"},"content":{"rendered":"<div><span style=\"font-family: inherit; font-size: 42px; font-style: inherit;\">Spatial Interpolation in Python: Kriging and IDW using PyKrige and SciPy<\/span><\/div>\n<div>\n<div>\n<div id=\"ez-toc-container\" class=\"ez-toc-v2_0_89 counter-hierarchy ez-toc-counter ez-toc-grey ez-toc-container-direction\">\n<div class=\"ez-toc-title-container\">\n<p class=\"ez-toc-title\" style=\"cursor:inherit\">Table of Contents<\/p>\n<span class=\"ez-toc-title-toggle\"><a href=\"#\" class=\"ez-toc-pull-right ez-toc-btn ez-toc-btn-xs ez-toc-btn-default ez-toc-toggle\" aria-label=\"Toggle Table of Content\"><span class=\"ez-toc-js-icon-con\"><span class=\"\"><span class=\"eztoc-hide\" style=\"display:none;\">Toggle<\/span><span class=\"ez-toc-icon-toggle-span\"><svg style=\"fill: #999;color:#999\" xmlns=\"http:\/\/www.w3.org\/2000\/svg\" class=\"list-377408\" width=\"20px\" height=\"20px\" viewBox=\"0 0 24 24\" fill=\"none\"><path d=\"M6 6H4v2h2V6zm14 0H8v2h12V6zM4 11h2v2H4v-2zm16 0H8v2h12v-2zM4 16h2v2H4v-2zm16 0H8v2h12v-2z\" fill=\"currentColor\"><\/path><\/svg><svg style=\"fill: #999;color:#999\" class=\"arrow-unsorted-368013\" xmlns=\"http:\/\/www.w3.org\/2000\/svg\" width=\"10px\" height=\"10px\" viewBox=\"0 0 24 24\" version=\"1.2\" baseProfile=\"tiny\"><path d=\"M18.2 9.3l-6.2-6.3-6.2 6.3c-.2.2-.3.4-.3.7s.1.5.3.7c.2.2.4.3.7.3h11c.3 0 .5-.1.7-.3.2-.2.3-.5.3-.7s-.1-.5-.3-.7zM5.8 14.7l6.2 6.3 6.2-6.3c.2-.2.3-.5.3-.7s-.1-.5-.3-.7c-.2-.2-.4-.3-.7-.3h-11c-.3 0-.5.1-.7.3-.2.2-.3.5-.3.7s.1.5.3.7z\"\/><\/svg><\/span><\/span><\/span><\/a><\/span><\/div>\n<nav><ul class='ez-toc-list ez-toc-list-level-1 ' ><li class='ez-toc-page-1 ez-toc-heading-level-2'><a class=\"ez-toc-link ez-toc-heading-1\" href=\"https:\/\/helpwithmygisassignment.com\/blog\/spatial-interpolation-in-python\/#Introduction\" >Introduction<\/a><\/li><li class='ez-toc-page-1 ez-toc-heading-level-2'><a class=\"ez-toc-link ez-toc-heading-2\" href=\"https:\/\/helpwithmygisassignment.com\/blog\/spatial-interpolation-in-python\/#1_Set_up_the_Python_environment_and_data\" >1. Set up the Python environment and data<\/a><ul class='ez-toc-list-level-3' ><li class='ez-toc-heading-level-3'><a class=\"ez-toc-link ez-toc-heading-3\" href=\"https:\/\/helpwithmygisassignment.com\/blog\/spatial-interpolation-in-python\/#11_Install_the_required_libraries\" >1.1 Install the required libraries<\/a><\/li><li class='ez-toc-page-1 ez-toc-heading-level-3'><a class=\"ez-toc-link ez-toc-heading-4\" href=\"https:\/\/helpwithmygisassignment.com\/blog\/spatial-interpolation-in-python\/#12_Prepare_the_input_data\" >1.2 Prepare the input data<\/a><\/li><li class='ez-toc-page-1 ez-toc-heading-level-3'><a class=\"ez-toc-link ez-toc-heading-5\" href=\"https:\/\/helpwithmygisassignment.com\/blog\/spatial-interpolation-in-python\/#13_Load_and_validate_the_observations\" >1.3 Load and validate the observations<\/a><\/li><\/ul><\/li><li class='ez-toc-page-1 ez-toc-heading-level-2'><a class=\"ez-toc-link ez-toc-heading-6\" href=\"https:\/\/helpwithmygisassignment.com\/blog\/spatial-interpolation-in-python\/#2_Inverse_Distance_Weighting_IDW_with_SciPy\" >2. Inverse Distance Weighting (IDW) with SciPy<\/a><ul class='ez-toc-list-level-3' ><li class='ez-toc-heading-level-3'><a class=\"ez-toc-link ez-toc-heading-7\" href=\"https:\/\/helpwithmygisassignment.com\/blog\/spatial-interpolation-in-python\/#21_Understand_the_IDW_equation\" >2.1 Understand the IDW equation<\/a><\/li><li class='ez-toc-page-1 ez-toc-heading-level-3'><a class=\"ez-toc-link ez-toc-heading-8\" href=\"https:\/\/helpwithmygisassignment.com\/blog\/spatial-interpolation-in-python\/#22_Implement_IDW\" >2.2 Implement IDW<\/a><\/li><li class='ez-toc-page-1 ez-toc-heading-level-3'><a class=\"ez-toc-link ez-toc-heading-9\" href=\"https:\/\/helpwithmygisassignment.com\/blog\/spatial-interpolation-in-python\/#23_Generate_the_prediction_grid\" >2.3 Generate the prediction grid<\/a><\/li><li class='ez-toc-page-1 ez-toc-heading-level-3'><a class=\"ez-toc-link ez-toc-heading-10\" href=\"https:\/\/helpwithmygisassignment.com\/blog\/spatial-interpolation-in-python\/#24_Calculate_the_IDW_surface\" >2.4 Calculate the IDW surface<\/a><\/li><\/ul><\/li><li class='ez-toc-page-1 ez-toc-heading-level-2'><a class=\"ez-toc-link ez-toc-heading-11\" href=\"https:\/\/helpwithmygisassignment.com\/blog\/spatial-interpolation-in-python\/#3_Ordinary_Kriging_with_PyKrige\" >3. Ordinary Kriging with PyKrige<\/a><ul class='ez-toc-list-level-3' ><li class='ez-toc-heading-level-3'><a class=\"ez-toc-link ez-toc-heading-12\" href=\"https:\/\/helpwithmygisassignment.com\/blog\/spatial-interpolation-in-python\/#31_Understand_the_variogram\" >3.1 Understand the variogram<\/a><\/li><li class='ez-toc-page-1 ez-toc-heading-level-3'><a class=\"ez-toc-link ez-toc-heading-13\" href=\"https:\/\/helpwithmygisassignment.com\/blog\/spatial-interpolation-in-python\/#32_Fit_and_execute_Ordinary_Kriging\" >3.2 Fit and execute Ordinary Kriging<\/a><\/li><li class='ez-toc-page-1 ez-toc-heading-level-3'><a class=\"ez-toc-link ez-toc-heading-14\" href=\"https:\/\/helpwithmygisassignment.com\/blog\/spatial-interpolation-in-python\/#33_Inspect_the_fitted_variogram\" >3.3 Inspect the fitted variogram<\/a><\/li><li class='ez-toc-page-1 ez-toc-heading-level-3'><a class=\"ez-toc-link ez-toc-heading-15\" href=\"https:\/\/helpwithmygisassignment.com\/blog\/spatial-interpolation-in-python\/#34_Interpret_kriging_variance_correctly\" >3.4 Interpret kriging variance correctly<\/a><\/li><\/ul><\/li><li class='ez-toc-page-1 ez-toc-heading-level-2'><a class=\"ez-toc-link ez-toc-heading-16\" href=\"https:\/\/helpwithmygisassignment.com\/blog\/spatial-interpolation-in-python\/#4_Visualize_the_prediction_surfaces\" >4. Visualize the prediction surfaces<\/a><\/li><li class='ez-toc-page-1 ez-toc-heading-level-2'><a class=\"ez-toc-link ez-toc-heading-17\" href=\"https:\/\/helpwithmygisassignment.com\/blog\/spatial-interpolation-in-python\/#5_Export_the_results_as_GeoTIFF\" >5. Export the results as GeoTIFF<\/a><ul class='ez-toc-list-level-3' ><li class='ez-toc-heading-level-3'><a class=\"ez-toc-link ez-toc-heading-18\" href=\"https:\/\/helpwithmygisassignment.com\/blog\/spatial-interpolation-in-python\/#51_Set_the_coordinate_reference_system\" >5.1 Set the coordinate reference system<\/a><\/li><li class='ez-toc-page-1 ez-toc-heading-level-3'><a class=\"ez-toc-link ez-toc-heading-19\" href=\"https:\/\/helpwithmygisassignment.com\/blog\/spatial-interpolation-in-python\/#52_Write_a_reusable_GeoTIFF_function\" >5.2 Write a reusable GeoTIFF function<\/a><\/li><li class='ez-toc-page-1 ez-toc-heading-level-3'><a class=\"ez-toc-link ez-toc-heading-20\" href=\"https:\/\/helpwithmygisassignment.com\/blog\/spatial-interpolation-in-python\/#53_Export_all_three_rasters\" >5.3 Export all three rasters<\/a><\/li><\/ul><\/li><li class='ez-toc-page-1 ez-toc-heading-level-2'><a class=\"ez-toc-link ez-toc-heading-21\" href=\"https:\/\/helpwithmygisassignment.com\/blog\/spatial-interpolation-in-python\/#6_Evaluate_model_accuracy_using_cross-validation\" >6. Evaluate model accuracy using cross-validation<\/a><ul class='ez-toc-list-level-3' ><li class='ez-toc-heading-level-3'><a class=\"ez-toc-link ez-toc-heading-22\" href=\"https:\/\/helpwithmygisassignment.com\/blog\/spatial-interpolation-in-python\/#61_Define_the_validation_functions\" >6.1 Define the validation functions<\/a><\/li><li class='ez-toc-page-1 ez-toc-heading-level-3'><a class=\"ez-toc-link ez-toc-heading-23\" href=\"https:\/\/helpwithmygisassignment.com\/blog\/spatial-interpolation-in-python\/#62_Evaluate_IDW_configurations\" >6.2 Evaluate IDW configurations<\/a><\/li><li class='ez-toc-page-1 ez-toc-heading-level-3'><a class=\"ez-toc-link ez-toc-heading-24\" href=\"https:\/\/helpwithmygisassignment.com\/blog\/spatial-interpolation-in-python\/#63_Evaluate_the_variogram_models\" >6.3 Evaluate the variogram models<\/a><\/li><\/ul><\/li><li class='ez-toc-page-1 ez-toc-heading-level-2'><a class=\"ez-toc-link ez-toc-heading-25\" href=\"https:\/\/helpwithmygisassignment.com\/blog\/spatial-interpolation-in-python\/#7_Statistical_assumptions_and_common_pitfalls\" >7. Statistical assumptions and common pitfalls<\/a><ul class='ez-toc-list-level-3' ><li class='ez-toc-heading-level-3'><a class=\"ez-toc-link ez-toc-heading-26\" href=\"https:\/\/helpwithmygisassignment.com\/blog\/spatial-interpolation-in-python\/#IDW_and_radial_basis_functions_are_not_the_same_method\" >IDW and radial basis functions are not the same method<\/a><\/li><li class='ez-toc-page-1 ez-toc-heading-level-3'><a class=\"ez-toc-link ez-toc-heading-27\" href=\"https:\/\/helpwithmygisassignment.com\/blog\/spatial-interpolation-in-python\/#Ordinary_Kriging_assumes_an_appropriate_spatial_model\" >Ordinary Kriging assumes an appropriate spatial model<\/a><\/li><li class='ez-toc-page-1 ez-toc-heading-level-3'><a class=\"ez-toc-link ez-toc-heading-28\" href=\"https:\/\/helpwithmygisassignment.com\/blog\/spatial-interpolation-in-python\/#Raster_resolution_is_not_predictive_accuracy\" >Raster resolution is not predictive accuracy<\/a><\/li><li class='ez-toc-page-1 ez-toc-heading-level-3'><a class=\"ez-toc-link ez-toc-heading-29\" href=\"https:\/\/helpwithmygisassignment.com\/blog\/spatial-interpolation-in-python\/#Extrapolation_can_be_misleading\" >Extrapolation can be misleading<\/a><\/li><\/ul><\/li><li class='ez-toc-page-1 ez-toc-heading-level-2'><a class=\"ez-toc-link ez-toc-heading-30\" href=\"https:\/\/helpwithmygisassignment.com\/blog\/spatial-interpolation-in-python\/#Conclusion\" >Conclusion<\/a><\/li><li class='ez-toc-page-1 ez-toc-heading-level-2'><a class=\"ez-toc-link ez-toc-heading-31\" href=\"https:\/\/helpwithmygisassignment.com\/blog\/spatial-interpolation-in-python\/#Official_documentation\" >Official documentation<\/a><\/li><\/ul><\/nav><\/div>\n<h2><span class=\"ez-toc-section\" id=\"Introduction\"><\/span>Introduction<span class=\"ez-toc-section-end\"><\/span><\/h2>\n<p><img loading=\"lazy\" decoding=\"async\" class=\"aligncenter size-large wp-image-1288\" src=\"https:\/\/helpwithmygisassignment.com\/blog\/wp-content\/uploads\/2026\/10\/Gemini_Generated_Image_4t3eax4t3eax4t3e-1024x572.jpg\" alt=\"Spatial Interpolation in Python transforming scatter points into continuous prediction surfaces.\" width=\"1024\" height=\"572\" srcset=\"https:\/\/helpwithmygisassignment.com\/blog\/wp-content\/uploads\/2026\/10\/Gemini_Generated_Image_4t3eax4t3eax4t3e-1024x572.jpg 1024w, https:\/\/helpwithmygisassignment.com\/blog\/wp-content\/uploads\/2026\/10\/Gemini_Generated_Image_4t3eax4t3eax4t3e-300x167.jpg 300w, https:\/\/helpwithmygisassignment.com\/blog\/wp-content\/uploads\/2026\/10\/Gemini_Generated_Image_4t3eax4t3eax4t3e-768x429.jpg 768w, https:\/\/helpwithmygisassignment.com\/blog\/wp-content\/uploads\/2026\/10\/Gemini_Generated_Image_4t3eax4t3eax4t3e.jpg 1376w\" sizes=\"auto, (max-width: 1024px) 100vw, 1024px\" \/><\/p>\n<p>Spatial interpolation transforms scattered point observations into a continuous surface by estimating values at locations where measurements are unavailable. It is widely used to map rainfall, groundwater levels, soil properties, air pollution, temperature, and other environmental variables.<\/p>\n<p>Performing <strong>Spatial Interpolation in Python using Kriging and IDW with PyKrige and SciPy<\/strong> makes these analyses repeatable and easier to integrate into automated geospatial workflows. Instead of manually running an interpolation tool for each dataset, a Python script can clean observations, generate a prediction grid, test different models, calculate validation statistics, and export georeferenced raster files.<\/p>\n<p>Two common methods are <a href=\"https:\/\/helpwithmygisassignment.com\/blog\/idw-vs-kriging-gis\/\">Inverse Distance Weighting (IDW) and Ordinary Kriging<\/a>.<\/p>\n<ul>\n<li><strong>IDW<\/strong> estimates unknown values using a distance-weighted average of observed measurements.<\/li>\n<li><strong>Ordinary Kriging<\/strong> uses a fitted variogram model to calculate spatial prediction weights and estimate uncertainty through kriging variance.<\/li>\n<\/ul>\n<p>The key difference is that IDW assumes that distance determines the influence of observations, whereas Kriging explicitly models spatial dependence. Neither method is universally superior; their performance depends on the data, the underlying spatial process, and the prediction task.<\/p>\n<p>This tutorial builds both models, compares their predictions, exports GeoTIFF surfaces, and evaluates accuracy using cross-validation.<\/p>\n<h2><span class=\"ez-toc-section\" id=\"1_Set_up_the_Python_environment_and_data\"><\/span>1. Set up the Python environment and data<span class=\"ez-toc-section-end\"><\/span><\/h2>\n<h3><span class=\"ez-toc-section\" id=\"11_Install_the_required_libraries\"><\/span>1.1 Install the required libraries<span class=\"ez-toc-section-end\"><\/span><\/h3>\n<p>Install the following packages in your Python environment:<\/p>\n<pre><code>pip install numpy pandas scipy pykrige matplotlib rasterio scikit-learn<\/code><\/pre>\n<p>The workflow uses NumPy and pandas for data preparation, SciPy for nearest-neighbour searches and IDW, PyKrige for Kriging, Matplotlib for visualization, Rasterio for GeoTIFF creation, and scikit-learn for cross-validation.<\/p>\n<h3><span class=\"ez-toc-section\" id=\"12_Prepare_the_input_data\"><\/span>1.2 Prepare the input data<span class=\"ez-toc-section-end\"><\/span><\/h3>\n<p>Create a CSV file named <code>observations.csv<\/code> containing three columns:<\/p>\n<table>\n<tbody>\n<tr>\n<th>Column<\/th>\n<th>Description<\/th>\n<\/tr>\n<tr>\n<td><code>x<\/code><\/td>\n<td>Projected easting or X coordinate<\/td>\n<\/tr>\n<tr>\n<td><code>y<\/code><\/td>\n<td>Projected northing or Y coordinate<\/td>\n<\/tr>\n<tr>\n<td><code>value<\/code><\/td>\n<td>Observed measurement<\/td>\n<\/tr>\n<\/tbody>\n<\/table>\n<p>The <code>value<\/code> field might represent rainfall in millimetres, groundwater depth in metres, or a soil measurement.<\/p>\n<p><strong>Coordinate systems matter.<\/strong> The methods in this tutorial use Euclidean distances, so the coordinates should normally be in an appropriate projected coordinate reference system with linear units. Do not apply planar distance calculations directly to raw longitude and latitude degrees without an appropriate geographic-distance approach.<\/p>\n<p>A CSV file does not inherently identify the CRS of its coordinates. Record the correct CRS from the source dataset before exporting results.<\/p>\n<h3><span class=\"ez-toc-section\" id=\"13_Load_and_validate_the_observations\"><\/span>1.3 Load and validate the observations<span class=\"ez-toc-section-end\"><\/span><\/h3>\n<pre><code>from pathlib import Path\r\n\r\nimport numpy as np\r\nimport pandas as pd\r\n\r\n# Define the input CSV and output directory.\r\ninput_csv = Path(\"observations.csv\")\r\noutput_dir = Path(\"interpolation_outputs\")\r\noutput_dir.mkdir(parents=True, exist_ok=True)\r\n\r\n# Load the measurement table.\r\ndf = pd.read_csv(input_csv)\r\n\r\n# Verify that the expected columns exist.\r\nrequired = [\"x\", \"y\", \"value\"]\r\nmissing = set(required) - set(df.columns)\r\n\r\nif missing:\r\n    raise ValueError(f\"Missing columns: {sorted(missing)}\")\r\n\r\n# Convert coordinate and measurement fields to numeric values.\r\ndf[required] = df[required].apply(pd.to_numeric, errors=\"coerce\")\r\n\r\n# Remove missing and infinite values.\r\ndf = df.replace([np.inf, -np.inf], np.nan)\r\ndf = df.dropna(subset=required)\r\n\r\n# Aggregate observations recorded at identical coordinates.\r\n# Review replicate measurements before using this step in production.\r\ndf = df.groupby([\"x\", \"y\"], as_index=False)[\"value\"].mean()\r\n\r\n# Extract the cleaned data as NumPy arrays.\r\nx = df[\"x\"].to_numpy(dtype=float)\r\ny = df[\"y\"].to_numpy(dtype=float)\r\nz = df[\"value\"].to_numpy(dtype=float)\r\n\r\n# Combine X and Y into a coordinate matrix.\r\nxy = np.column_stack((x, y))\r\n\r\n# Check the number and distribution of unique sample locations.\r\nif len(df) &lt; 10:\r\n    raise ValueError(\r\n        \"Too few unique locations for this example workflow. \"\r\n        \"Use more observations or reassess the validation design.\"\r\n    )\r\n\r\nif np.ptp(x) == 0 or np.ptp(y) == 0:\r\n    raise ValueError(\r\n        \"The sample locations must span both coordinate directions.\"\r\n    )\r\n\r\nprint(\"Unique sample locations:\", len(df))\r\nprint(\"Observed minimum:\", z.min())\r\nprint(\"Observed maximum:\", z.max())<\/code><\/pre>\n<p>Averaging duplicate coordinates is a practical way to prevent ambiguous zero-distance calculations in IDW. However, if repeated measurements represent meaningful temporal variation or measurement uncertainty, they should be assessed rather than automatically averaged.<\/p>\n<p>The minimum sample-count check is a basic safeguard, not a guarantee that the dataset is large or representative enough for reliable Kriging. Inspect the spatial distribution, measurement units, outliers, and sampling gaps before proceeding.<\/p>\n<h2><span class=\"ez-toc-section\" id=\"2_Inverse_Distance_Weighting_IDW_with_SciPy\"><\/span>2. Inverse Distance Weighting (IDW) with SciPy<span class=\"ez-toc-section-end\"><\/span><\/h2>\n<h3><span class=\"ez-toc-section\" id=\"21_Understand_the_IDW_equation\"><\/span>2.1 Understand the IDW equation<span class=\"ez-toc-section-end\"><\/span><\/h3>\n<p>IDW estimates an unknown value by assigning higher weights to closer observations:<\/p>\n<p>[<br \/>\n\\hat{z}(s_0)=<br \/>\n\\frac{\\sum_{i=1}^{n}d_i^{-p}z(s_i)}<br \/>\n{\\sum_{i=1}^{n}d_i^{-p}}<br \/>\n]<\/p>\n<p>Where:<\/p>\n<ul>\n<li>(\\hat{z}(s_0)) is the predicted value at the target location.<\/li>\n<li>(z(s_i)) is the observed value at sample (i).<\/li>\n<li>(d_i) is the distance between the sample and the target.<\/li>\n<li>(p) is the distance-decay power.<\/li>\n<li>(n) is the number of observations used in the prediction.<\/li>\n<\/ul>\n<p>When (p=2), the weight is proportional to inverse squared distance. Increasing the power makes nearby observations more influential.<\/p>\n<p>SciPy does not provide a dedicated standard IDW estimator. Instead, the following implementation uses <code>scipy.spatial.cKDTree<\/code> for efficient nearest-neighbour searches and calculates the IDW weights directly. SciPy documents <code>cKDTree<\/code> as a spatial indexing structure that supports efficient nearest-neighbour queries (<a href=\"https:\/\/docs.scipy.org\/doc\/scipy\/reference\/generated\/scipy.spatial.cKDTree.query.html\" target=\"_blank\" rel=\"noopener\">official documentation<\/a>).<\/p>\n<h3><span class=\"ez-toc-section\" id=\"22_Implement_IDW\"><\/span>2.2 Implement IDW<span class=\"ez-toc-section-end\"><\/span><\/h3>\n<pre><code>from scipy.spatial import cKDTree\r\n\r\n\r\ndef idw_predict(xy, values, query_xy, power=2.0, k=12, tree=None):\r\n    \"\"\"\r\n    Estimate values at query locations using inverse distance weighting.\r\n\r\n    xy:       Known sample coordinates, shape (n_samples, 2).\r\n    values:   Measurements at known coordinates.\r\n    query_xy: Target coordinates, shape (n_targets, 2).\r\n    power:    Positive inverse-distance power.\r\n    k:        Number of nearest observations to use.\r\n    tree:     Optional pre-built cKDTree for repeated predictions.\r\n    \"\"\"\r\n\r\n    # Validate the interpolation parameters.\r\n    if power &lt;= 0:\r\n        raise ValueError(\"Power must be greater than zero.\")\r\n\r\n    if k &lt; 1:\r\n        raise ValueError(\"Neighbour count must be at least one.\")\r\n\r\n    # Limit the neighbour count to the available samples.\r\n    k_eff = min(int(k), len(values))\r\n\r\n    if k_eff == 0:\r\n        raise ValueError(\"No valid sample values were supplied.\")\r\n\r\n    # Build the spatial index only when one was not supplied.\r\n    if tree is None:\r\n        tree = cKDTree(xy)\r\n\r\n    # Find the nearest sample points for each target location.\r\n    distances, indices = tree.query(query_xy, k=k_eff)\r\n\r\n    # Preserve a consistent array shape for k=1.\r\n    if k_eff == 1:\r\n        distances = distances[:, None]\r\n        indices = indices[:, None]\r\n\r\n    # Calculate inverse-distance weights.\r\n    weights = 1.0 \/ np.maximum(distances, 1e-12) ** power\r\n\r\n    # Retrieve the observed values at the neighbouring indices.\r\n    neighbour_values = values[indices]\r\n\r\n    # Compute the weighted mean for every target.\r\n    predictions = (\r\n        np.sum(weights * neighbour_values, axis=1)\r\n        \/ np.sum(weights, axis=1)\r\n    )\r\n\r\n    # At an existing sample location, return its measured value.\r\n    exact = distances == 0.0\r\n    rows = np.flatnonzero(exact.any(axis=1))\r\n\r\n    if rows.size:\r\n        first = exact[rows].argmax(axis=1)\r\n        predictions[rows] = values[indices[rows, first]]\r\n\r\n    return predictions<\/code><\/pre>\n<p>The function accepts an optional spatial tree so that a large prediction grid can be processed in batches without rebuilding the index each time.<\/p>\n<p>The small distance guard prevents division by zero, while the explicit exact-match rule ensures that a prediction at a known sample location returns the measured value.<\/p>\n<h3><span class=\"ez-toc-section\" id=\"23_Generate_the_prediction_grid\"><\/span>2.3 Generate the prediction grid<span class=\"ez-toc-section-end\"><\/span><\/h3>\n<p>The next step is to create a regular grid of cell-centre coordinates. The example uses a 100-unit cell size; choose a resolution that suits the data density and intended output scale.<\/p>\n<pre><code># Define the raster cell size in the units of the projected CRS.\r\nCELL_SIZE = 100.0\r\n\r\n# Snap the raster boundaries to multiples of the cell size.\r\nxmin = np.floor(x.min() \/ CELL_SIZE) * CELL_SIZE\r\nxmax = np.ceil(x.max() \/ CELL_SIZE) * CELL_SIZE\r\nymin = np.floor(y.min() \/ CELL_SIZE) * CELL_SIZE\r\nymax = np.ceil(y.max() \/ CELL_SIZE) * CELL_SIZE\r\n\r\n# Ensure non-zero width and height.\r\nif xmax &lt;= xmin:\r\n    xmax = xmin + CELL_SIZE\r\n\r\nif ymax &lt;= ymin:\r\n    ymax = ymin + CELL_SIZE\r\n\r\n# Calculate raster dimensions.\r\nwidth = int(round((xmax - xmin) \/ CELL_SIZE))\r\nheight = int(round((ymax - ymin) \/ CELL_SIZE))\r\n\r\n# Create coordinates at the centre of each raster cell.\r\ngrid_x = xmin + (np.arange(width) + 0.5) * CELL_SIZE\r\ngrid_y = ymin + (np.arange(height) + 0.5) * CELL_SIZE\r\n\r\n# Create the two-dimensional coordinate arrays.\r\nX, Y = np.meshgrid(grid_x, grid_y)\r\n\r\n# Flatten the grid into an array of target coordinate pairs.\r\nquery_xy = np.column_stack((X.ravel(), Y.ravel()))\r\n\r\n# Store the extent for plotting.\r\nplot_extent = (xmin, xmax, ymin, ymax)\r\n\r\nprint(f\"Raster dimensions: {width} columns \u00d7 {height} rows\")\r\nprint(f\"Grid cell size: {CELL_SIZE}\")<\/code><\/pre>\n<p>This grid covers the rectangular extent of the observations. That does not mean every cell is equally supported by the data. For final analysis, apply a study-area boundary mask and consider masking or flagging predictions far from sampled locations.<\/p>\n<h3><span class=\"ez-toc-section\" id=\"24_Calculate_the_IDW_surface\"><\/span>2.4 Calculate the IDW surface<span class=\"ez-toc-section-end\"><\/span><\/h3>\n<pre><code># Build the nearest-neighbour index once.\r\ntree = cKDTree(xy)\r\n\r\n# Allocate memory for all predicted cell values.\r\nidw_flat = np.empty(len(query_xy), dtype=float)\r\n\r\n# Predict in chunks to limit temporary memory use.\r\nchunk_size = 50_000\r\n\r\nfor start in range(0, len(query_xy), chunk_size):\r\n    stop = min(start + chunk_size, len(query_xy))\r\n\r\n    idw_flat[start:stop] = idw_predict(\r\n        xy=xy,\r\n        values=z,\r\n        query_xy=query_xy[start:stop],\r\n        power=2.0,\r\n        k=12,\r\n        tree=tree,\r\n    )\r\n\r\n# Convert the predictions into a two-dimensional grid.\r\nidw_grid = idw_flat.reshape(height, width)\r\n\r\nprint(\"IDW prediction range:\",\r\n      np.nanmin(idw_grid),\r\n      np.nanmax(idw_grid))<\/code><\/pre>\n<p>The output is a continuous surface estimated from the point observations.<\/p>\n<p><strong>Statistical limitation:<\/strong> IDW does not estimate a variogram or provide a kriging variance. Its performance depends heavily on the chosen power, neighbour count, sample distribution, and spatial structure of the variable. Select these settings using validation rather than visual smoothness alone.<\/p>\n<h2><span class=\"ez-toc-section\" id=\"3_Ordinary_Kriging_with_PyKrige\"><\/span>3. Ordinary Kriging with PyKrige<span class=\"ez-toc-section-end\"><\/span><\/h2>\n<h3><span class=\"ez-toc-section\" id=\"31_Understand_the_variogram\"><\/span>3.1 Understand the variogram<span class=\"ez-toc-section-end\"><\/span><\/h3>\n<p>Ordinary Kriging accounts for spatial dependence through a variogram, which describes how dissimilar observations tend to become as the distance between them increases.<\/p>\n<p>The experimental semivariogram is commonly written as:<\/p>\n<p>[<br \/>\n\\gamma(h)=<br \/>\n\\frac{1}{2N(h)}<br \/>\n\\sum_{(i,j)\\in N(h)}<br \/>\n[z(s_i)-z(s_j)]^2<br \/>\n]<\/p>\n<p>Here, (h) is a separation-distance bin and (N(h)) represents the observation pairs assigned to that bin.<\/p>\n<p>A theoretical variogram model is fitted to the experimental semivariogram. PyKrige can estimate model parameters automatically when they are not supplied explicitly (<a href=\"https:\/\/geostat-framework.readthedocs.io\/projects\/pykrige\/en\/stable\/examples\/00_ordinary.html\" target=\"_blank\" rel=\"noopener\">PyKrige Ordinary Kriging example<\/a>).<\/p>\n<p>Three commonly used models are:<\/p>\n<table>\n<tbody>\n<tr>\n<th>Model<\/th>\n<th>General behaviour<\/th>\n<th>Typical consideration<\/th>\n<\/tr>\n<tr>\n<td>Spherical<\/td>\n<td>Reaches the sill at a finite range<\/td>\n<td>Useful when dependence appears to level off at a distinct distance<\/td>\n<\/tr>\n<tr>\n<td>Exponential<\/td>\n<td>Approaches the sill asymptotically<\/td>\n<td>Represents spatial dependence that diminishes gradually<\/td>\n<\/tr>\n<tr>\n<td>Gaussian<\/td>\n<td>Smooth increase near the origin<\/td>\n<td>Can represent a relatively smooth spatial field<\/td>\n<\/tr>\n<\/tbody>\n<\/table>\n<p>The <strong>nugget<\/strong> reflects measurement error and\/or unresolved short-range variability. The <strong>sill<\/strong> represents the level that semivariance approaches, while the <strong>range<\/strong> describes the characteristic scale of spatial dependence. The exact parameter interpretation depends on the model and its parameterization.<\/p>\n<h3><span class=\"ez-toc-section\" id=\"32_Fit_and_execute_Ordinary_Kriging\"><\/span>3.2 Fit and execute Ordinary Kriging<span class=\"ez-toc-section-end\"><\/span><\/h3>\n<pre><code>from pykrige.ok import OrdinaryKriging\r\n\r\n# Select a variogram model for the initial run.\r\nVARIOGRAM_MODEL = \"spherical\"\r\n\r\n# Build the Ordinary Kriging model.\r\nok = OrdinaryKriging(\r\n    x,\r\n    y,\r\n    z,\r\n    variogram_model=VARIOGRAM_MODEL,\r\n    nlags=6,\r\n    weight=False,\r\n    verbose=False,\r\n    enable_plotting=False,\r\n    exact_values=True,\r\n    pseudo_inv=True,\r\n)\r\n\r\n# Predict the value at each grid cell.\r\n# PyKrige also returns the associated kriging variance.\r\nkriging_result, variance_result = ok.execute(\r\n    \"grid\",\r\n    grid_x,\r\n    grid_y,\r\n    backend=\"vectorized\",\r\n)\r\n\r\n# Convert masked outputs into arrays containing NaN where masked.\r\nkriging_grid = np.asarray(\r\n    np.ma.filled(kriging_result, np.nan),\r\n    dtype=float,\r\n)\r\n\r\nkriging_variance = np.asarray(\r\n    np.ma.filled(variance_result, np.nan),\r\n    dtype=float,\r\n)\r\n\r\n# Check that the result matches the grid dimensions.\r\nif kriging_grid.shape != (height, width):\r\n    raise RuntimeError(\"Unexpected Kriging raster dimensions.\")\r\n\r\nif kriging_variance.shape != (height, width):\r\n    raise RuntimeError(\"Unexpected Kriging variance dimensions.\")\r\n\r\nprint(\"Kriging prediction range:\",\r\n      np.nanmin(kriging_grid),\r\n      np.nanmax(kriging_grid))\r\n\r\nprint(\"Kriging variance range:\",\r\n      np.nanmin(kriging_variance),\r\n      np.nanmax(kriging_variance))<\/code><\/pre>\n<p>When <code>variogram_parameters<\/code> is omitted, PyKrige estimates the parameters from the experimental semivariogram. Here, <code>nlags=6<\/code> defines the number of distance bins, while <code>weight=False<\/code> leaves the automatic fitting routine&#8217;s optional short-lag weighting disabled.<\/p>\n<p>The <code>exact_values=True<\/code> setting treats observations as exact measurements for prediction purposes. If measurements contain substantial error, the assumptions should be reconsidered rather than blindly relying on exact interpolation.<\/p>\n<p>The vectorized backend can be fast, but its memory requirements may become large for fine grids. For a large raster, test tiled prediction or use a lower-memory execution strategy supported by your installed PyKrige version.<\/p>\n<h3><span class=\"ez-toc-section\" id=\"33_Inspect_the_fitted_variogram\"><\/span>3.3 Inspect the fitted variogram<span class=\"ez-toc-section-end\"><\/span><\/h3>\n<pre><code>import matplotlib.pyplot as plt\r\n\r\n# Display the experimental variogram and fitted theoretical model.\r\nok.display_variogram_model()\r\n\r\nplt.show()<\/code><\/pre>\n<p>Check whether the fitted curve follows the general pattern of the empirical semivariogram. A poor fit, an implausible range, or a weakly defined sill may indicate that the selected model or assumptions are inappropriate.<\/p>\n<p>Repeat the model construction and prediction for the alternative models:<\/p>\n<pre><code># Alternative variogram names:\r\n# \"exponential\"\r\n# \"gaussian\"<\/code><\/pre>\n<p>Do not choose the final variogram from its appearance alone. Compare candidate models using the same validation design.<\/p>\n<h3><span class=\"ez-toc-section\" id=\"34_Interpret_kriging_variance_correctly\"><\/span>3.4 Interpret kriging variance correctly<span class=\"ez-toc-section-end\"><\/span><\/h3>\n<p>Ordinary Kriging returns a predicted value and a kriging variance at each output location.<\/p>\n<p>The variance describes uncertainty conditional on the fitted model and its assumptions. It is not the same as the actual error of the prediction. It also does not automatically capture all uncertainty associated with measurement quality, variogram selection, or an incorrectly specified spatial process.<\/p>\n<p>For a well-supported model, variance often reflects the geometry of the sampling network and the assumed spatial dependence. Use it alongside validation statistics and data-coverage information, not as a guarantee of predictive accuracy.<\/p>\n<h2><span class=\"ez-toc-section\" id=\"4_Visualize_the_prediction_surfaces\"><\/span>4. Visualize the prediction surfaces<span class=\"ez-toc-section-end\"><\/span><\/h2>\n<p><img loading=\"lazy\" decoding=\"async\" class=\"aligncenter size-large wp-image-1289\" src=\"https:\/\/helpwithmygisassignment.com\/blog\/wp-content\/uploads\/2026\/10\/Gemini_Generated_Image_2x8si2x8si2x8si2-1024x572.jpg\" alt=\"Spatial Interpolation in Python comparing Inverse Distance Weighting and Ordinary Kriging prediction maps.\" width=\"1024\" height=\"572\" srcset=\"https:\/\/helpwithmygisassignment.com\/blog\/wp-content\/uploads\/2026\/10\/Gemini_Generated_Image_2x8si2x8si2x8si2-1024x572.jpg 1024w, https:\/\/helpwithmygisassignment.com\/blog\/wp-content\/uploads\/2026\/10\/Gemini_Generated_Image_2x8si2x8si2x8si2-300x167.jpg 300w, https:\/\/helpwithmygisassignment.com\/blog\/wp-content\/uploads\/2026\/10\/Gemini_Generated_Image_2x8si2x8si2x8si2-768x429.jpg 768w, https:\/\/helpwithmygisassignment.com\/blog\/wp-content\/uploads\/2026\/10\/Gemini_Generated_Image_2x8si2x8si2x8si2.jpg 1376w\" sizes=\"auto, (max-width: 1024px) 100vw, 1024px\" \/><\/p>\n<p>A side-by-side map helps reveal differences between IDW and Kriging. Use the same value scale for both prediction maps so that colour differences represent the same numerical differences.<\/p>\n<pre><code>import matplotlib.pyplot as plt\r\n\r\n# Derive common colour limits for the two prediction surfaces.\r\nvmin = min(np.nanmin(idw_grid), np.nanmin(kriging_grid))\r\nvmax = max(np.nanmax(idw_grid), np.nanmax(kriging_grid))\r\n\r\n# Create two prediction panels and one variance panel.\r\nfig, axes = plt.subplots(\r\n    1, 3,\r\n    figsize=(16, 5),\r\n    constrained_layout=True,\r\n)\r\n\r\n# Plot IDW and Kriging using identical colour limits.\r\nfor ax, surface, title in [\r\n    (axes[0], idw_grid, \"IDW Prediction\"),\r\n    (axes[1], kriging_grid, \"Ordinary Kriging\"),\r\n]:\r\n    image = ax.imshow(\r\n        surface,\r\n        origin=\"lower\",\r\n        extent=plot_extent,\r\n        aspect=\"equal\",\r\n        cmap=\"viridis\",\r\n        vmin=vmin,\r\n        vmax=vmax,\r\n    )\r\n\r\n    # Mark observed locations without confusing them\r\n    # with the raster's numerical colour scale.\r\n    ax.scatter(\r\n        x, y,\r\n        marker=\"o\",\r\n        facecolors=\"none\",\r\n        edgecolors=\"black\",\r\n        s=24,\r\n    )\r\n\r\n    ax.set_title(title)\r\n    ax.set_xlabel(\"Projected X\")\r\n    ax.set_ylabel(\"Projected Y\")\r\n    fig.colorbar(image, ax=ax, label=\"Predicted value\")\r\n\r\n# Plot kriging variance with its own numerical scale.\r\nvariance_image = axes[2].imshow(\r\n    kriging_variance,\r\n    origin=\"lower\",\r\n    extent=plot_extent,\r\n    aspect=\"equal\",\r\n    cmap=\"magma\",\r\n)\r\n\r\naxes[2].scatter(\r\n    x, y,\r\n    marker=\"o\",\r\n    facecolors=\"none\",\r\n    edgecolors=\"black\",\r\n    s=24,\r\n)\r\n\r\naxes[2].set_title(\"Kriging Variance\")\r\naxes[2].set_xlabel(\"Projected X\")\r\naxes[2].set_ylabel(\"Projected Y\")\r\n\r\nfig.colorbar(\r\n    variance_image,\r\n    ax=axes[2],\r\n    label=\"Kriging variance\",\r\n)\r\n\r\n# Save the comparison figure.\r\nfig.savefig(\r\n    output_dir \/ \"interpolation_comparison.png\",\r\n    dpi=200,\r\n    bbox_inches=\"tight\",\r\n)\r\n\r\nplt.show()<\/code><\/pre>\n<p>Interpret the maps in context. IDW may display localized peaks around measured points, while Kriging can produce a different pattern depending on the fitted variogram. Neither pattern alone establishes which model is better.<\/p>\n<h2><span class=\"ez-toc-section\" id=\"5_Export_the_results_as_GeoTIFF\"><\/span>5. Export the results as GeoTIFF<span class=\"ez-toc-section-end\"><\/span><\/h2>\n<p>The NumPy arrays contain predictions but not their geographic reference information. Rasterio writes the numerical values alongside the raster dimensions, affine transform, CRS, and NoData definition. See the <a href=\"https:\/\/rasterio.readthedocs.io\/en\/stable\/topics\/writing.html\" target=\"_blank\" rel=\"noopener\">Rasterio writing documentation<\/a>.<\/p>\n<h3><span class=\"ez-toc-section\" id=\"51_Set_the_coordinate_reference_system\"><\/span>5.1 Set the coordinate reference system<span class=\"ez-toc-section-end\"><\/span><\/h3>\n<p>For illustration, the following setting uses a projected CRS identifier. <strong>Replace it with the actual CRS of the input coordinates.<\/strong> A CRS identifier must not be guessed from the coordinate magnitudes alone.<\/p>\n<pre><code>import rasterio\r\nfrom rasterio.transform import from_origin\r\n\r\n# Example only: change to the CRS used by your input X and Y coordinates.\r\nDATA_CRS = \"EPSG:32737\"\r\n\r\n# Define the transform from the upper-left corner.\r\ntransform = from_origin(\r\n    xmin,\r\n    ymax,\r\n    CELL_SIZE,\r\n    CELL_SIZE,\r\n)<\/code><\/pre>\n<p>Because <code>grid_y<\/code> was created from south to north, the arrays must be flipped vertically when written to a north-up GeoTIFF. The transform describes the upper-left pixel corner and pixel dimensions; see the <a href=\"https:\/\/rasterio.readthedocs.io\/en\/stable\/api\/rasterio.transform.html\" target=\"_blank\" rel=\"noopener\">Rasterio transform documentation<\/a>.<\/p>\n<h3><span class=\"ez-toc-section\" id=\"52_Write_a_reusable_GeoTIFF_function\"><\/span>5.2 Write a reusable GeoTIFF function<span class=\"ez-toc-section-end\"><\/span><\/h3>\n<pre><code>def write_geotiff(output_path, array):\r\n    \"\"\"Write a 2D float array to a single-band GeoTIFF.\"\"\"\r\n\r\n    # Convert the prediction array to floating-point values.\r\n    data = np.asarray(array, dtype=np.float32)\r\n\r\n    # Confirm that the array matches the expected grid.\r\n    if data.shape != (height, width):\r\n        raise ValueError(\"Array dimensions do not match the output grid.\")\r\n\r\n    # Convert from bottom-up array order to north-up raster order.\r\n    data = np.flipud(data)\r\n\r\n    # Replace non-finite values with an explicit NoData value.\r\n    nodata = -9999.0\r\n\r\n    data = np.where(\r\n        np.isfinite(data),\r\n        data,\r\n        nodata,\r\n    ).astype(np.float32)\r\n\r\n    # Define the GeoTIFF metadata.\r\n    profile = {\r\n        \"driver\": \"GTiff\",\r\n        \"height\": height,\r\n        \"width\": width,\r\n        \"count\": 1,\r\n        \"dtype\": \"float32\",\r\n        \"crs\": DATA_CRS,\r\n        \"transform\": transform,\r\n        \"nodata\": nodata,\r\n        \"compress\": \"deflate\",\r\n    }\r\n\r\n    # Write the raster using a context manager.\r\n    with rasterio.open(output_path, \"w\", **profile) as dst:\r\n        dst.write(data, 1)\r\n\r\n    print(f\"Saved: {output_path}\")<\/code><\/pre>\n<h3><span class=\"ez-toc-section\" id=\"53_Export_all_three_rasters\"><\/span>5.3 Export all three rasters<span class=\"ez-toc-section-end\"><\/span><\/h3>\n<pre><code># Export the IDW prediction surface.\r\nwrite_geotiff(\r\n    output_dir \/ \"idw_surface.tif\",\r\n    idw_grid,\r\n)\r\n\r\n# Export the Ordinary Kriging prediction surface.\r\nwrite_geotiff(\r\n    output_dir \/ \"ordinary_kriging.tif\",\r\n    kriging_grid,\r\n)\r\n\r\n# Export the associated kriging variance surface.\r\nwrite_geotiff(\r\n    output_dir \/ \"kriging_variance.tif\",\r\n    kriging_variance,\r\n)<\/code><\/pre>\n<p>Open the files in QGIS or ArcGIS Pro and verify their CRS, cell size, grid alignment, extent, and NoData values.<\/p>\n<p>This workflow creates a rectangular raster around the sample network. For a production map, apply a study-area mask if one is available and assess whether areas outside the supported sampling footprint should be excluded or flagged.<\/p>\n<h2><span class=\"ez-toc-section\" id=\"6_Evaluate_model_accuracy_using_cross-validation\"><\/span>6. Evaluate model accuracy using cross-validation<span class=\"ez-toc-section-end\"><\/span><\/h2>\n<p>An interpolation surface may look convincing while still producing large prediction errors. Cross-validation tests the model by withholding known observations, predicting their values from the remaining observations, and comparing predictions with measurements.<\/p>\n<p>The Root Mean Square Error is:<\/p>\n<p>[<br \/>\nRMSE=\\sqrt{\\frac{1}{n}<br \/>\n\\sum_{i=1}^{n}(z_i-\\hat{z}_i)^2}<br \/>\n]<\/p>\n<p>Lower RMSE indicates smaller squared prediction errors on the validation observations, but it does not by itself establish that the model is unbiased or reliable throughout the study area.<\/p>\n<h3><span class=\"ez-toc-section\" id=\"61_Define_the_validation_functions\"><\/span>6.1 Define the validation functions<span class=\"ez-toc-section-end\"><\/span><\/h3>\n<p>The following workflow compares several IDW settings and three Kriging variogram models using the same five-fold partition.<\/p>\n<pre><code>from sklearn.model_selection import KFold\r\n\r\n\r\ndef rmse(observed, predicted):\r\n    \"\"\"Calculate RMSE while ignoring invalid prediction pairs.\"\"\"\r\n\r\n    observed = np.asarray(observed, dtype=float)\r\n    predicted = np.asarray(predicted, dtype=float)\r\n\r\n    valid = np.isfinite(observed) &amp; np.isfinite(predicted)\r\n\r\n    if not valid.any():\r\n        return np.nan\r\n\r\n    return float(\r\n        np.sqrt(np.mean((observed[valid] - predicted[valid]) ** 2))\r\n    )\r\n\r\n\r\n# Create reproducible training and validation folds.\r\ncv = KFold(\r\n    n_splits=5,\r\n    shuffle=True,\r\n    random_state=42,\r\n)\r\n\r\n# Store the splits so each model uses the same partitions.\r\nsplits = list(cv.split(xy))<\/code><\/pre>\n<h3><span class=\"ez-toc-section\" id=\"62_Evaluate_IDW_configurations\"><\/span>6.2 Evaluate IDW configurations<span class=\"ez-toc-section-end\"><\/span><\/h3>\n<pre><code>def validate_idw(power, k):\r\n    \"\"\"Return cross-validated RMSE for one IDW configuration.\"\"\"\r\n\r\n    predictions = np.full(len(z), np.nan, dtype=float)\r\n\r\n    for train_idx, test_idx in splits:\r\n        # Fit IDW using only the training observations.\r\n        train_xy = xy[train_idx]\r\n        train_z = z[train_idx]\r\n\r\n        # Predict at the withheld coordinates.\r\n        predictions[test_idx] = idw_predict(\r\n            xy=train_xy,\r\n            values=train_z,\r\n            query_xy=xy[test_idx],\r\n            power=power,\r\n            k=k,\r\n        )\r\n\r\n    # Compare all withheld predictions with their observations.\r\n    return rmse(z, predictions)\r\n\r\n\r\n# Compare different powers and neighbourhood sizes.\r\nidw_results = []\r\n\r\nfor power in [1.0, 2.0, 3.0]:\r\n    for k in [8, 12, 16]:\r\n        score = validate_idw(power, k)\r\n\r\n        idw_results.append({\r\n            \"model\": \"IDW\",\r\n            \"power\": power,\r\n            \"neighbors\": k,\r\n            \"variogram\": None,\r\n            \"RMSE\": score,\r\n        })<\/code><\/pre>\n<p>This parameter search helps determine whether stronger distance decay or a different neighbourhood size improves prediction accuracy.<\/p>\n<h3><span class=\"ez-toc-section\" id=\"63_Evaluate_the_variogram_models\"><\/span>6.3 Evaluate the variogram models<span class=\"ez-toc-section-end\"><\/span><\/h3>\n<pre><code>def validate_kriging(model_name):\r\n    \"\"\"Return cross-validated RMSE for a variogram model.\"\"\"\r\n\r\n    predictions = np.full(len(z), np.nan, dtype=float)\r\n\r\n    for train_idx, test_idx in splits:\r\n        # Use only training locations for variogram fitting.\r\n        train_xy = xy[train_idx]\r\n        train_z = z[train_idx]\r\n\r\n        test_xy = xy[test_idx]\r\n\r\n        # Construct the fold-specific Kriging model.\r\n        model = OrdinaryKriging(\r\n            train_xy[:, 0],\r\n            train_xy[:, 1],\r\n            train_z,\r\n            variogram_model=model_name,\r\n            nlags=6,\r\n            weight=False,\r\n            verbose=False,\r\n            enable_plotting=False,\r\n            exact_values=True,\r\n            pseudo_inv=True,\r\n        )\r\n\r\n        # Predict only at the withheld coordinates.\r\n        predicted, _ = model.execute(\r\n            \"points\",\r\n            test_xy[:, 0],\r\n            test_xy[:, 1],\r\n        )\r\n\r\n        predictions[test_idx] = np.asarray(\r\n            np.ma.filled(predicted, np.nan),\r\n            dtype=float,\r\n        ).ravel()\r\n\r\n    return rmse(z, predictions)\r\n\r\n\r\n# Evaluate spherical, exponential, and Gaussian variograms.\r\nkriging_results = []\r\n\r\nfor model_name in [\"spherical\", \"exponential\", \"gaussian\"]:\r\n    score = validate_kriging(model_name)\r\n\r\n    kriging_results.append({\r\n        \"model\": \"Ordinary Kriging\",\r\n        \"power\": None,\r\n        \"neighbors\": None,\r\n        \"variogram\": model_name,\r\n        \"RMSE\": score,\r\n    })\r\n\r\n# Combine and rank all evaluated configurations.\r\nresults = pd.DataFrame(idw_results + kriging_results)\r\n\r\nresults = results.sort_values(\"RMSE\", na_position=\"last\")\r\n\r\nprint(results.to_string(index=False))<\/code><\/pre>\n<p>The code fits the variogram inside each training fold instead of fitting it once using all observations. This avoids letting the withheld measurements influence the variogram fitted for that fold.<\/p>\n<p>The comparison reports RMSE for the evaluated configurations. The lowest score identifies the best-performing tested configuration under this particular validation design; it does not prove that the selected model is best for every location or future dataset.<\/p>\n<p><strong>Important spatial-validation limitation:<\/strong> Randomly shuffled folds can overestimate performance when neighbouring locations appear in both training and validation sets. For spatial prediction, consider spatial block cross-validation, where geographically separated blocks are withheld together. Choose a block size appropriate to the study extent and spatial dependence. See the <a href=\"https:\/\/scikit-learn.org\/stable\/modules\/cross_validation.html\" target=\"_blank\" rel=\"noopener\">scikit-learn cross-validation guide<\/a> for the general validation framework.<\/p>\n<h2><span class=\"ez-toc-section\" id=\"7_Statistical_assumptions_and_common_pitfalls\"><\/span>7. Statistical assumptions and common pitfalls<span class=\"ez-toc-section-end\"><\/span><\/h2>\n<h3><span class=\"ez-toc-section\" id=\"IDW_and_radial_basis_functions_are_not_the_same_method\"><\/span>IDW and radial basis functions are not the same method<span class=\"ez-toc-section-end\"><\/span><\/h3>\n<p>Although both methods are available in SciPy-related interpolation workflows, IDW and radial basis function interpolation use different mathematical formulations.<\/p>\n<p><code>scipy.interpolate.Rbf<\/code> is a legacy interface. For new implementations, SciPy recommends <code>RBFInterpolator<\/code>. RBF methods combine radial kernels and can include a polynomial component; they are not simply inverse-distance weighted averages. See the <a href=\"https:\/\/docs.scipy.org\/doc\/scipy\/reference\/generated\/scipy.interpolate.RBFInterpolator.html\" target=\"_blank\" rel=\"noopener\">SciPy RBFInterpolator reference<\/a>.<\/p>\n<h3><span class=\"ez-toc-section\" id=\"Ordinary_Kriging_assumes_an_appropriate_spatial_model\"><\/span>Ordinary Kriging assumes an appropriate spatial model<span class=\"ez-toc-section-end\"><\/span><\/h3>\n<p>Ordinary Kriging assumes a locally constant but unknown mean and a suitable variogram describing spatial dependence. A strong regional trend, anisotropy, abrupt geological boundaries, or spatial barriers may violate the assumptions of a simple isotropic model.<\/p>\n<p>Investigate the empirical variogram and consider alternative modelling approaches where the data do not support ordinary Kriging.<\/p>\n<h3><span class=\"ez-toc-section\" id=\"Raster_resolution_is_not_predictive_accuracy\"><\/span>Raster resolution is not predictive accuracy<span class=\"ez-toc-section-end\"><\/span><\/h3>\n<p>A 10-metre raster cell does not mean that predictions are accurate to 10 metres. Resolution controls the representation of the surface, while accuracy depends on sampling density, measurement quality, spatial structure, model assumptions, and validation.<\/p>\n<h3><span class=\"ez-toc-section\" id=\"Extrapolation_can_be_misleading\"><\/span>Extrapolation can be misleading<span class=\"ez-toc-section-end\"><\/span><\/h3>\n<p>The algorithms can produce predictions away from sample points, but these estimates may be poorly supported. Mask the final raster to the study area where appropriate, inspect the distance to the nearest observation, and avoid overinterpreting large unsampled regions.<\/p>\n<h2><span class=\"ez-toc-section\" id=\"Conclusion\"><\/span>Conclusion<span class=\"ez-toc-section-end\"><\/span><\/h2>\n<p>Python provides a reproducible approach to spatial interpolation that can be integrated into environmental monitoring, hydrological analysis, and automated geospatial data pipelines.<\/p>\n<p>The workflow presented here prepares point data, constructs a regular prediction grid, runs IDW using SciPy, fits Ordinary Kriging models using PyKrige, visualizes both prediction surfaces, exports GeoTIFFs with spatial metadata, and compares model accuracy through cross-validation.<\/p>\n<p>The most important distinction is that <strong>IDW uses a distance-decay rule, whereas Kriging models spatial dependence through a variogram<\/strong>. Their suitability must be evaluated using the characteristics of the observations and an appropriate validation strategy.<\/p>\n<p>A reliable interpolation workflow is not simply one that produces a smooth map. It is one that preserves coordinate reference information, makes its assumptions explicit, evaluates predictive performance, and communicates where its estimates are\u2014and are not\u2014well supported.<\/p>\n<h2><span class=\"ez-toc-section\" id=\"Official_documentation\"><\/span>Official documentation<span class=\"ez-toc-section-end\"><\/span><\/h2>\n<ul>\n<li><a href=\"https:\/\/geostat-framework.readthedocs.io\/projects\/pykrige\/en\/stable\/examples\/00_ordinary.html\" target=\"_blank\" rel=\"noopener\">PyKrige Ordinary Kriging example<\/a><\/li>\n<li><a href=\"https:\/\/geostat-framework.readthedocs.io\/projects\/pykrige\/en\/latest\/api.html\" target=\"_blank\" rel=\"noopener\">PyKrige API documentation<\/a><\/li>\n<li><a href=\"https:\/\/docs.scipy.org\/doc\/scipy\/reference\/generated\/scipy.spatial.cKDTree.query.html\" target=\"_blank\" rel=\"noopener\">SciPy cKDTree nearest-neighbour queries<\/a><\/li>\n<li><a href=\"https:\/\/docs.scipy.org\/doc\/scipy\/reference\/generated\/scipy.interpolate.RBFInterpolator.html\" target=\"_blank\" rel=\"noopener\">SciPy RBFInterpolator<\/a><\/li>\n<li><a href=\"https:\/\/rasterio.readthedocs.io\/en\/stable\/topics\/writing.html\" target=\"_blank\" rel=\"noopener\">Rasterio writing datasets<\/a><\/li>\n<li><a href=\"https:\/\/rasterio.readthedocs.io\/en\/stable\/api\/rasterio.transform.html\" target=\"_blank\" rel=\"noopener\">Rasterio coordinate transforms<\/a><\/li>\n<li><a href=\"https:\/\/scikit-learn.org\/stable\/modules\/cross_validation.html\" target=\"_blank\" rel=\"noopener\">scikit-learn cross-validation<\/a><\/li>\n<\/ul>\n<\/div>\n<p>The biggest improvement is that the article now has a consistent path from raw point data to validated GeoTIFF outputs, with explicit warnings where the statistical assumptions or coordinate metadata matter. I also removed the implication that RBF interpolation is another form of IDW and made the model comparison more meaningful by using the same validation folds.<\/p>\n<\/div>\n","protected":false},"excerpt":{"rendered":"<p>Spatial Interpolation in Python: Kriging and IDW using PyKrige and SciPy Introduction Spatial interpolation transforms scattered point observations into a continuous surface by estimating values at locations where measurements are unavailable. It is widely used to map rainfall, groundwater levels, soil properties, air pollution, temperature, and other environmental variables. Performing Spatial Interpolation in Python using &#8230; <a title=\"Spatial Interpolation in Python: Kriging and IDW using PyKrige and SciPy\" class=\"read-more\" href=\"https:\/\/helpwithmygisassignment.com\/blog\/spatial-interpolation-in-python\/\" aria-label=\"Read more about Spatial Interpolation in Python: Kriging and IDW using PyKrige and SciPy\">Read more<\/a><\/p>\n","protected":false},"author":1,"featured_media":1288,"comment_status":"closed","ping_status":"closed","sticky":false,"template":"","format":"standard","meta":{"footnotes":""},"categories":[17,4],"tags":[],"class_list":["post-1195","post","type-post","status-publish","format-standard","has-post-thumbnail","hentry","category-civil-engineering","category-gis-articles"],"_links":{"self":[{"href":"https:\/\/helpwithmygisassignment.com\/blog\/wp-json\/wp\/v2\/posts\/1195","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/helpwithmygisassignment.com\/blog\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/helpwithmygisassignment.com\/blog\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/helpwithmygisassignment.com\/blog\/wp-json\/wp\/v2\/users\/1"}],"replies":[{"embeddable":true,"href":"https:\/\/helpwithmygisassignment.com\/blog\/wp-json\/wp\/v2\/comments?post=1195"}],"version-history":[{"count":6,"href":"https:\/\/helpwithmygisassignment.com\/blog\/wp-json\/wp\/v2\/posts\/1195\/revisions"}],"predecessor-version":[{"id":1291,"href":"https:\/\/helpwithmygisassignment.com\/blog\/wp-json\/wp\/v2\/posts\/1195\/revisions\/1291"}],"wp:featuredmedia":[{"embeddable":true,"href":"https:\/\/helpwithmygisassignment.com\/blog\/wp-json\/wp\/v2\/media\/1288"}],"wp:attachment":[{"href":"https:\/\/helpwithmygisassignment.com\/blog\/wp-json\/wp\/v2\/media?parent=1195"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/helpwithmygisassignment.com\/blog\/wp-json\/wp\/v2\/categories?post=1195"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/helpwithmygisassignment.com\/blog\/wp-json\/wp\/v2\/tags?post=1195"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}