Add spatial interpolation with xr_interpolate notebook (#1233)

* Add ensemble tide modelling functionality to model_tides

* Update test_coastal.py

* Remove test

* Updates to IDW, xr_interpolate and ensemble tide modelling code"

* Doco updates

* Switch ensemble rankings from high to low = good

* Update docstring

* Fix doco

* Add interpolation notebook

* Remove coastal files from branch

* Add points data

* Review feedback;

* Add p param to IDW

* Fix test

How_to_guides/README.rst

@@ -21,6 +21,7 @@ A recipe book of simple code examples demonstrating how to perform common geospa
    Generating_geomedian_composites.ipynb
    Image_segmentation.ipynb
    Imagery_on_web_map.ipynb
+   Interpolation.ipynb
    Land_cover_animated_plots.ipynb
    Land_cover_change_mapping.ipynb
    Land_cover_pixel_drill.ipynb

Supplementary_data/Interpolation/interpolation_points.geojson

@@ -0,0 +1,24 @@
+{
+"type": "FeatureCollection",
+"name": "interpolation_points",
+"crs": { "type": "name", "properties": { "name": "urn:ogc:def:crs:OGC:1.3:CRS84" } },
+"features": [
+{ "type": "Feature", "properties": { "z": 1 }, "geometry": { "type": "Point", "coordinates": [ 149.041358543960143, -35.381110059879738 ] } },
+{ "type": "Feature", "properties": { "z": 2 }, "geometry": { "type": "Point", "coordinates": [ 149.050043769856671, -35.376996005507706 ] } },
+{ "type": "Feature", "properties": { "z": 3 }, "geometry": { "type": "Point", "coordinates": [ 149.069699807411894, -35.371053482525888 ] } },
+{ "type": "Feature", "properties": { "z": 4 }, "geometry": { "type": "Point", "coordinates": [ 149.077470799003521, -35.363739608086725 ] } },
+{ "type": "Feature", "properties": { "z": 5 }, "geometry": { "type": "Point", "coordinates": [ 149.081584853375546, -35.345454921988825 ] } },
+{ "type": "Feature", "properties": { "z": 6 }, "geometry": { "type": "Point", "coordinates": [ 149.087070259204921, -35.3125424870126 ] } },
+{ "type": "Feature", "properties": { "z": 7 }, "geometry": { "type": "Point", "coordinates": [ 149.092098547881847, -35.286029692170644 ] } },
+{ "type": "Feature", "properties": { "z": 8 }, "geometry": { "type": "Point", "coordinates": [ 149.098041070863644, -35.268659240377637 ] } },
+{ "type": "Feature", "properties": { "z": 9 }, "geometry": { "type": "Point", "coordinates": [ 149.108554765369945, -35.256317077261549 ] } },
+{ "type": "Feature", "properties": { "z": 10 }, "geometry": { "type": "Point", "coordinates": [ 149.129582154382547, -35.251288788584631 ] } },
+{ "type": "Feature", "properties": { "z": 11 }, "geometry": { "type": "Point", "coordinates": [ 149.148781074785319, -35.252660140041975 ] } },
+{ "type": "Feature", "properties": { "z": 12 }, "geometry": { "type": "Point", "coordinates": [ 149.162494589358772, -35.259516897328687 ] } },
+{ "type": "Feature", "properties": { "z": 13 }, "geometry": { "type": "Point", "coordinates": [ 149.180779275456672, -35.257688428718893 ] } },
+{ "type": "Feature", "properties": { "z": 14 }, "geometry": { "type": "Point", "coordinates": [ 149.190378735658072, -35.247174734212599 ] } },
+{ "type": "Feature", "properties": { "z": 15 }, "geometry": { "type": "Point", "coordinates": [ 149.198149727249671, -35.226147345200012 ] } },
+{ "type": "Feature", "properties": { "z": 16 }, "geometry": { "type": "Point", "coordinates": [ 149.213691710432897, -35.211976713474137 ] } },
+{ "type": "Feature", "properties": { "z": 17 }, "geometry": { "type": "Point", "coordinates": [ 149.229233693616095, -35.20374860473008 ] } }
+]
+}

Tests/dea_tools/test_spatial.py

@@ -13,6 +13,7 @@
     xr_vectorize,
     xr_rasterize,
     xr_interpolate,
+    idw,
 )
 from dea_tools.validation import eval_metrics
 
@@ -335,13 +336,84 @@ def test_subpixel_contours_dim(satellite_da):
     subpixel_contours(satellite_da_date, z_values=600, dim="date")
 
 
-# def test_subpixel_contours_dem_crs(dem_da):
-#     # Verify that an error is raised if data passed with no spatial ref/geobox
-#     with pytest.raises(ValueError):
-#         subpixel_contours(dem_da.drop_vars("spatial_ref"), z_values=700)
-
-#     # Verify that no error is raised if we provide the correct CRS
-#     subpixel_contours(dem_da.drop_vars("spatial_ref"), z_values=700, crs="EPSG:4326")
+# Test Inverse Distance Weighted function
+def test_idw():
+    # Basic psuedo-1D example
+    input_z = [1, 2, 3, 4, 5]
+    input_x = [0, 1, 2, 3, 4]
+    input_y = [0, 0, 0, 0, 0]
+    output_x = [0.5, 1.5, 2.5, 3.5]
+    output_y = [0.0, 0.0, 0.0, 0.0]
+    out = idw(input_z, input_x, input_y, output_x, output_y, k=2)
+    assert np.allclose(out, [1.5, 2.5, 3.5, 4.5])
+
+    # Verify that k > input points gives error
+    with pytest.raises(ValueError):
+        idw(input_z, input_x, input_y, output_x, output_y, k=6)
+
+    # 2D nearest neighbour case
+    input_z = [1, 2, 3, 4]
+    input_x = [0, 4, 0, 4]
+    input_y = [0, 0, 4, 4]
+    output_x = [1, 4, 0, 3]
+    output_y = [0, 1, 3, 4]
+    out = idw(input_z, input_x, input_y, output_x, output_y, k=1)
+    assert np.allclose(out, [1, 2, 3, 4])
+
+    # Two neighbours
+    input_z = [1, 2, 3, 4]
+    input_x = [0, 4, 0, 4]
+    input_y = [0, 0, 4, 4]
+    output_x = [2, 0, 4, 2]
+    output_y = [0, 2, 2, 4]
+    out = idw(input_z, input_x, input_y, output_x, output_y, k=2)
+    assert np.allclose(out, [1.5, 2, 3, 3.5])
+
+    # Four neighbours
+    out = idw(input_z, input_x, input_y, output_x, output_y, k=4)
+    assert np.allclose(out, [2.11, 2.30, 2.69, 2.88], rtol=0.01)
+
+    # Four neighbours; max distance of 2
+    out = idw(input_z, input_x, input_y, output_x, output_y, k=4, max_dist=2)
+    assert np.allclose(out, [1.5, 2, 3, 3.5])
+
+    # Four neighbours; max distance of 2, k_min of 4 (should return NaN)
+    out = idw(input_z, input_x, input_y, output_x, output_y, k=4, max_dist=2, k_min=4)
+    assert np.isnan(out).all()
+
+    # Four neighbours; power function p=0
+    out = idw(input_z, input_x, input_y, output_x, output_y, k=4, p=0)
+    assert np.allclose(out, [2.5, 2.5, 2.5, 2.5])
+
+    # Four neighbours; power function p=2
+    out = idw(input_z, input_x, input_y, output_x, output_y, k=4, p=2)
+    assert np.allclose(out, [1.83, 2.17, 2.83, 3.17], rtol=0.01)
+
+    # Different units, nearest neighbour case
+    input_z = [10, 20, 30, 40]
+    input_x = [1125296, 1155530, 1125296, 1155530]
+    input_y = [-4169722, -4169722, -4214782, -4214782]
+    output_x = [1124952, 1159593, 1120439, 1155284]
+    output_y = [-4169749, -4172892, -4211108, -4214332]
+    out = idw(input_z, input_x, input_y, output_x, output_y, k=1)
+    assert np.allclose(out, [10, 20, 30, 40])
+
+    # Verify distance works on different units
+    output_x = [1142134, 1138930]
+    output_y = [-4171232, -4213451]
+    out = idw(input_z, input_x, input_y, output_x, output_y, k=4, max_dist=20000)
+    assert np.allclose(out, [15, 35], rtol=0.1)
+
+    # Test multidimensional input
+    input_z = np.column_stack(([1, 2, 3, 4], [10, 20, 30, 40]))
+    input_x = [0, 4, 0, 4]
+    input_y = [0, 0, 4, 4]
+    output_x = [1, 4, 0, 3]
+    output_y = [0, 1, 3, 4]
+    out = idw(input_z, input_x, input_y, output_x, output_y, k=1)
+    assert input_z.shape == out.shape
+    assert np.allclose(out[:, 0], [1, 2, 3, 4])
+    assert np.allclose(out[:, 1], [10, 20, 30, 40])
 
 
 @pytest.mark.parametrize(

Tests/test_notebooks.sh

@@ -11,6 +11,6 @@ pip3 install ./Tools
 pytest Tests/dea_tools
 
 # Test Juputer Notebooks
-pytest --durations=10 --nbval-lax Beginners_guide DEA_products --ignore DEA_products/DEA_Wetlands_Insight_Tool.ipynb How_to_guides/Animated_timeseries.ipynb How_to_guides/Contour_extraction.ipynb How_to_guides/Calculating_band_indices.ipynb How_to_guides/Downloading_data_with_STAC.ipynb How_to_guides/Exporting_GeoTIFFs.ipynb How_to_guides/Generating_composites.ipynb How_to_guides/Image_segmentation.ipynb How_to_guides/Opening_GeoTIFFs_NetCDFs.ipynb How_to_guides/Pansharpening.ipynb How_to_guides/Planetary_computer.ipynb How_to_guides/Polygon_drill.ipynb How_to_guides/Principal_component_analysis.ipynb How_to_guides/Rasterize_vectorize.ipynb How_to_guides/Tidal_modelling.ipynb How_to_guides/Using_load_ard.ipynb How_to_guides/Virtual_products.ipynb Real_world_examples/Coastal_erosion.ipynb Real_world_examples/Intertidal_elevation.ipynb
+pytest --durations=10 --nbval-lax Beginners_guide DEA_products --ignore DEA_products/DEA_Wetlands_Insight_Tool.ipynb How_to_guides/Animated_timeseries.ipynb How_to_guides/Contour_extraction.ipynb How_to_guides/Calculating_band_indices.ipynb How_to_guides/Downloading_data_with_STAC.ipynb How_to_guides/Exporting_GeoTIFFs.ipynb How_to_guides/Generating_composites.ipynb How_to_guides/Image_segmentation.ipynb How_to_guides/Interpolation.ipynb How_to_guides/Opening_GeoTIFFs_NetCDFs.ipynb How_to_guides/Pansharpening.ipynb How_to_guides/Planetary_computer.ipynb How_to_guides/Polygon_drill.ipynb How_to_guides/Principal_component_analysis.ipynb How_to_guides/Rasterize_vectorize.ipynb How_to_guides/Tidal_modelling.ipynb How_to_guides/Using_load_ard.ipynb How_to_guides/Virtual_products.ipynb Real_world_examples/Coastal_erosion.ipynb Real_world_examples/Intertidal_elevation.ipynb