Filter out missing data when interpolating from point data (#10)

* Filter out missing data when interpolating from point data * use compression to reduce filesizes for test data
corteva · Dec 12, 2019 · 4cc1705 · 4cc1705
1 parent 744d49e
commit 4cc1705
Show file tree

Hide file tree

Showing 21 changed files with 81 additions and 4 deletions.
diff --git a/geocube/rasterize.py b/geocube/rasterize.py
@@ -3,6 +3,7 @@
 This module contains tools for rasterizing vector data.
 """
 import numpy
+import pandas
 import rasterio.features
 import rasterio.transform
 import rasterio.warp
@@ -65,6 +66,19 @@ def rasterize_image(
         raise
 
 
+def _remove_missing_data(data_values, geometry_array):
+    """
+    Missing data causes issues with interpolation of point data
+    https://github.com/corteva/geocube/issues/9
+
+    This filters the data so those issues don't cause problems.
+    """
+    not_missing_data = ~pandas.isnull(data_values)
+    geometry_array = geometry_array[not_missing_data]
+    data_values = data_values[not_missing_data]
+    return data_values, geometry_array
+
+
 def rasterize_points_griddata(
     geometry_array,
     data_values,
@@ -104,6 +118,7 @@ def rasterize_points_griddata(
     if data_values.dtype == object:
         return None
     try:
+        data_values, geometry_array = _remove_missing_data(data_values, geometry_array)
         return griddata(
             points=(geometry_array.x, geometry_array.y),
             values=data_values,
@@ -152,6 +167,7 @@ def rasterize_points_radial(
     logger = get_logger()
 
     try:
+        data_values, geometry_array = _remove_missing_data(data_values, geometry_array)
         interp = Rbf(geometry_array.x, geometry_array.y, data_values, function=method)
         return interp(*numpy.meshgrid(grid_coords["x"], grid_coords["y"]))
     except ValueError as ter:

diff --git a/geocube/vector_to_cube.py b/geocube/vector_to_cube.py
@@ -5,17 +5,17 @@
 import numpy
 import pandas
 import xarray
-
-from geocube.geo_utils.geobox import load_vector_data
-from geocube.logger import get_logger
-from geocube.rasterize import rasterize_image
 from rioxarray.rioxarray import (
     DEFAULT_GRID_MAP,
     add_spatial_ref,
     add_xy_grid_meta,
     affine_to_coords,
 )
 
+from geocube.geo_utils.geobox import load_vector_data
+from geocube.logger import get_logger
+from geocube.rasterize import rasterize_image
+
 
 def _format_series_data(data_series):
     """

diff --git a/sphinx/history.rst b/sphinx/history.rst
@@ -1,6 +1,10 @@
 History
 =======
 
+0.0.10
+------
+- Filter out missing data when interpolating from point data (issue #9)
+
 0.0.9
 -----
 - Added `rescale` kwarg to `geocube.rasterize.rasterize_points_griddata`. (pull #8)

diff --git a/test/integration/api/test_core_integration.py b/test/integration/api/test_core_integration.py
@@ -711,3 +711,34 @@ def test_make_geocube__custom_rasterize_function(function, compare_name, tmpdir)
         os.path.join(TEST_COMPARE_DATA_DIR, compare_name), mask_and_scale=False
     ) as xdc:
         xarray.testing.assert_allclose(out_grid, xdc, rtol=0.1, atol=0.1)
+
+
+@pytest.mark.parametrize(
+    "function,compare_name",
+    [
+        (rasterize_points_griddata, "rasterize_griddata_nearest_nodata.nc"),
+        (
+            partial(rasterize_points_griddata, method="cubic"),
+            "rasterize_griddata_cubic_nodata.nc",
+        ),
+        (rasterize_points_radial, "rasterize_radial_linear_nodata.nc"),
+    ],
+)
+def test_make_geocube__custom_rasterize_function__filter_null(
+    function, compare_name, tmpdir
+):
+    input_geodata = os.path.join(TEST_INPUT_DATA_DIR, "point_with_null.geojson")
+    out_grid = make_geocube(
+        vector_data=input_geodata,
+        resolution=(-0.00001, 0.00001),
+        rasterize_function=function,
+    )
+
+    # test writing to netCDF
+    out_grid.to_netcdf(str(tmpdir.mkdir("geocube_custom").join(compare_name)))
+
+    # test output data
+    with xarray.open_dataset(
+        os.path.join(TEST_COMPARE_DATA_DIR, compare_name), mask_and_scale=False
+    ) as xdc:
+        xarray.testing.assert_allclose(out_grid, xdc, rtol=0.1, atol=0.1)
diff --git a/test/test_data/compare/rasterize_griddata_cubic.nc b/test/test_data/compare/rasterize_griddata_cubic.nc
diff --git a/test/test_data/compare/rasterize_griddata_cubic_nodata.nc b/test/test_data/compare/rasterize_griddata_cubic_nodata.nc
diff --git a/test/test_data/compare/rasterize_griddata_nearest.nc b/test/test_data/compare/rasterize_griddata_nearest.nc
diff --git a/test/test_data/compare/rasterize_griddata_nearest_nodata.nc b/test/test_data/compare/rasterize_griddata_nearest_nodata.nc
diff --git a/test/test_data/compare/rasterize_griddata_rescale.nc b/test/test_data/compare/rasterize_griddata_rescale.nc
diff --git a/test/test_data/compare/rasterize_image_sum.nc b/test/test_data/compare/rasterize_image_sum.nc
diff --git a/test/test_data/compare/rasterize_radial_linear.nc b/test/test_data/compare/rasterize_radial_linear.nc
diff --git a/test/test_data/compare/rasterize_radial_linear_nodata.nc b/test/test_data/compare/rasterize_radial_linear_nodata.nc
diff --git a/test/test_data/compare/soil_grid_flat_categorical.nc b/test/test_data/compare/soil_grid_flat_categorical.nc
diff --git a/test/test_data/compare/soil_grid_flat_interpolate_na.nc b/test/test_data/compare/soil_grid_flat_interpolate_na.nc
diff --git a/test/test_data/compare/soil_grid_group.nc b/test/test_data/compare/soil_grid_group.nc
diff --git a/test/test_data/compare/soil_grid_group_categorical.nc b/test/test_data/compare/soil_grid_group_categorical.nc
diff --git a/test/test_data/compare/soil_grid_group_no_geom.nc b/test/test_data/compare/soil_grid_group_no_geom.nc
diff --git a/test/test_data/compare/time_vector_data.nc b/test/test_data/compare/time_vector_data.nc
diff --git a/test/test_data/compare/vector_data_group.nc b/test/test_data/compare/vector_data_group.nc
diff --git a/test/test_data/compare/vector_time_data_group.nc b/test/test_data/compare/vector_time_data_group.nc
diff --git a/test/test_data/input/point_with_null.geojson b/test/test_data/input/point_with_null.geojson
@@ -0,0 +1,26 @@
+{
+"type": "FeatureCollection",
+"crs": { "type": "name", "properties": { "name": "urn:ogc:def:crs:OGC:1.3:CRS84" } },
+"features": [
+{ "type": "Feature", "properties": { "test_attr": 1.3}, "geometry": { "type": "Point", "coordinates": [ -47.266811, 44.219318 ] } },
+{ "type": "Feature", "properties": { "test_attr": 1.3}, "geometry": { "type": "Point", "coordinates": [ -47.266798, 44.219318 ] } },
+{ "type": "Feature", "properties": { "test_attr": null}, "geometry": { "type": "Point", "coordinates": [ -47.266807, 44.219318 ] } },
+{ "type": "Feature", "properties": { "test_attr": 1.3}, "geometry": { "type": "Point", "coordinates": [ -47.266803, 44.219318 ] } },
+{ "type": "Feature", "properties": { "test_attr": 1.3}, "geometry": { "type": "Point", "coordinates": [ -47.266794, 44.219318 ] } },
+{ "type": "Feature", "properties": { "test_attr": 1.3}, "geometry": { "type": "Point", "coordinates": [ -47.266789, 45.219318 ] } },
+{ "type": "Feature", "properties": { "test_attr": 5.3}, "geometry": { "type": "Point", "coordinates": [ -47.266816, 45.219318 ] } },
+{ "type": "Feature", "properties": { "test_attr": 1.3}, "geometry": { "type": "Point", "coordinates": [ -47.266851, 45.219318 ] } },
+{ "type": "Feature", "properties": { "test_attr": 1.3}, "geometry": { "type": "Point", "coordinates": [ -47.266842, 45.219318 ] } },
+{ "type": "Feature", "properties": { "test_attr": 1.3}, "geometry": { "type": "Point", "coordinates": [ -47.266825, 45.219318 ] } },
+{ "type": "Feature", "properties": { "test_attr": null}, "geometry": { "type": "Point", "coordinates": [ -47.266834, 44.219318 ] } },
+{ "type": "Feature", "properties": { "test_attr": 1.3}, "geometry": { "type": "Point", "coordinates": [ -47.26682, 44.219318 ] } },
+{ "type": "Feature", "properties": { "test_attr": 2.3}, "geometry": { "type": "Point", "coordinates": [ -47.266829, 44.219318 ] } },
+{ "type": "Feature", "properties": { "test_attr": 1.3}, "geometry": { "type": "Point", "coordinates": [ -47.266838, 44.219318 ] } },
+{ "type": "Feature", "properties": { "test_attr": 1.3}, "geometry": { "type": "Point", "coordinates": [ -47.266856, 44.219318 ] } },
+{ "type": "Feature", "properties": { "test_attr": 1.3}, "geometry": { "type": "Point", "coordinates": [ -47.266847, 45.219318 ] } },
+{ "type": "Feature", "properties": { "test_attr": 1.3}, "geometry": { "type": "Point", "coordinates": [ -47.266891, 45.219318 ] } },
+{ "type": "Feature", "properties": { "test_attr": null}, "geometry": { "type": "Point", "coordinates": [ -47.266865, 45.219318 ] } },
+{ "type": "Feature", "properties": { "test_attr": 1.3}, "geometry": { "type": "Point", "coordinates": [ -47.266873, 45.219318 ] } },
+{ "type": "Feature", "properties": { "test_attr": 1.7}, "geometry": { "type": "Point", "coordinates": [ -47.266869, 45.219318 ] } }
+]
+}