Merge branch 'geo-master' into percent_cover

Author: sgoodm

CHANGELOG.txt

@@ -1,5 +1,6 @@
-0.12dev
+0.12.0
 - zone_func argument to apply a function to the masked array before computing stats
+- support shapely 1.6 exceptions
 
 0.11
 - change `band_num` parameter to `band`

src/rasterstats/_version.py

@@ -1 +1 @@
-__version__ = "0.12a1"
+__version__ = "0.12.0"

src/rasterstats/io.py

@@ -10,7 +10,10 @@
 from rasterio.transform import guard_transform
 from affine import Affine
 import numpy as np
-from shapely.geos import ReadingError
+try:
+    from shapely.errors import ReadingError
+except:
+    from shapely.geos import ReadingError
 from shapely import wkt, wkb
 from collections import Iterable, Mapping
 

src/rasterstats/main.py

@@ -161,8 +161,8 @@ def gen_zonal_stats(
         warnings.warn("Use `geojson_out` to preserve feature properties",
                       DeprecationWarning)
 
-    bn = kwargs.get('band_num')
-    if bn:
+    band_num = kwargs.get('band_num')
+    if band_num:
         warnings.warn("Use `band` to specify band number", DeprecationWarning)
         band = band_num
 
@@ -238,8 +238,9 @@ def gen_zonal_stats(
             isnodata = (fsrc.array == fsrc.nodata)
 
             # add nan mask (if necessary)
-            if np.issubdtype(fsrc.array.dtype, float) and \
-               np.isnan(fsrc.array.min()):
+            has_nan = (np.issubdtype(fsrc.array.dtype, float)
+                and np.isnan(fsrc.array.min()))
+            if has_nan:
                 isnodata = (isnodata | np.isnan(fsrc.array))
 
             # Mask the source data array
@@ -325,9 +326,13 @@ def gen_zonal_stats(
                     pctarr = masked.compressed()
                     feature_stats[pctile] = np.percentile(pctarr, q)
 
-            if 'nodata' in stats:
-                featmasked = np.ma.MaskedArray(fsrc.array, mask=np.logical_not(rv_array))
-                feature_stats['nodata'] = float((featmasked == fsrc.nodata).sum())
+            if 'nodata' in stats or 'nan' in stats:
+                featmasked = np.ma.MaskedArray(fsrc.array, mask=(~rv_array))
+
+                if 'nodata' in stats:
+                    feature_stats['nodata'] = float((featmasked == fsrc.nodata).sum())
+                if 'nan' in stats:
+                    feature_stats['nan'] = float(np.isnan(featmasked).sum()) if has_nan else 0
 
             if add_stats is not None:
                 for stat_name, stat_func in add_stats.items():

src/rasterstats/point.py

@@ -53,7 +53,7 @@ def bilinear(arr, x, y):
     if hasattr(arr, 'count') and arr.count() != 4:
         # a masked array with at least one nodata
         # fall back to nearest neighbor
-        val = arr[round(1 - y), round(x)]
+        val = arr[int(round(1 - y)), int(round(x))]
         if val is masked:
             return None
         else:
@@ -172,7 +172,7 @@ def gen_point_query(
             for x, y in geom_xys(geom):
                 if interpolate == 'nearest':
                     r, c = rast.index(x, y)
-                    window = ((r, r+1), (c, c+1))
+                    window = ((int(r), int(r+1)), (int(c), int(c+1)))
                     src_array = rast.read(window=window, masked=True).array
                     val = src_array[0, 0]
                     if val is masked:

src/rasterstats/utils.py

@@ -11,7 +11,7 @@
 
 DEFAULT_STATS = ['count', 'min', 'max', 'mean']
 VALID_STATS = DEFAULT_STATS + \
-    ['sum', 'std', 'median', 'majority', 'minority', 'unique', 'range', 'nodata']
+    ['sum', 'std', 'median', 'majority', 'minority', 'unique', 'range', 'nodata', 'nan']
 #  also percentile_{q} but that is handled as special case
 
 

tests/test_cli.py

@@ -34,14 +34,15 @@ def test_cli_feature_stdin():
     warnings.simplefilter('ignore')
     result = runner.invoke(zonalstats,
                            ['--raster', raster,
-                            '--stats', 'mean',
+                            '--stats', 'all',
                             '--prefix', 'test_'],
                            input=open(vector, 'r').read())
     assert result.exit_code == 0
     outdata = json.loads(result.output)
     assert len(outdata['features']) == 1
     feature = outdata['features'][0]
     assert 'test_mean' in feature['properties']
+    assert 'test_std' in feature['properties']
 
 
 def test_cli_features_sequence():

tests/test_io.py

@@ -268,6 +268,11 @@ def test_Raster():
 
     r2 = Raster(arr, affine, nodata, band=1).read(bounds)
 
+    with pytest.raises(ValueError):
+        r3 = Raster(arr, affine, nodata, band=1).read()
+    with pytest.raises(ValueError):
+        r4 = Raster(arr, affine, nodata, band=1).read(bounds=1, window=1)
+
     # If the abstraction is correct, the arrays are equal
     assert np.array_equal(r1.array, r2.array)
 

tests/test_zonal.py

@@ -29,6 +29,16 @@ def test_main():
     assert round(stats[0]['mean'], 2) == 14.66
 
 
+# remove after band_num alias is removed
+def test_band_alias():
+    polygons = os.path.join(DATA, 'polygons.shp')
+    stats_a = zonal_stats(polygons, raster)
+    stats_b = zonal_stats(polygons, raster, band=1)
+    with pytest.deprecated_call():
+        stats_c = zonal_stats(polygons, raster, band_num=1)
+    assert stats_a[0]['count'] == stats_b[0]['count'] == stats_c[0]['count']
+
+
 def test_zonal_global_extent():
     polygons = os.path.join(DATA, 'polygons.shp')
     stats = zonal_stats(polygons, raster)
@@ -376,6 +386,27 @@ def test_some_nodata():
     assert stats[1]['nodata'] == 19
     assert stats[1]['count'] == 31
 
+
+# update this if nan end up being incorporated into nodata
+def test_nan_nodata():
+    polygon = Polygon([[0, 0], [2, 0], [2, 2], [0, 2]])
+    arr = np.array([
+        [np.nan, 12.25],
+        [-999, 12.75]
+    ])
+    affine = Affine(1, 0, 0,
+                    0, -1, 2)
+
+    stats = zonal_stats(polygon, arr, affine=affine, nodata=-999,
+                        stats='nodata count sum mean min max')
+
+    assert stats[0]['nodata'] == 1
+    assert stats[0]['count'] == 2
+    assert stats[0]['mean'] == 12.5
+    assert stats[0]['min'] == 12.25
+    assert stats[0]['max'] == 12.75
+
+
 def test_some_nodata_ndarray():
     polygons = os.path.join(DATA, 'polygons.shp')
     raster = os.path.join(DATA, 'slope_nodata.tif')
@@ -427,6 +458,67 @@ def test_geojson_out():
         assert 'count' in feature['properties']  # from zonal stats
 
 
+# do not think this is actually testing the line i wanted it to
+# since the read_features func for this data type is generating
+# the properties field
+def test_geojson_out_with_no_properties():
+    polygon = Polygon([[0, 0], [0, 0,5], [1, 1.5], [1.5, 2], [2, 2], [2, 0]])
+    arr = np.array([
+        [100, 1],
+        [100, 1]
+    ])
+    affine = Affine(1, 0, 0,
+                    0, -1, 2)
+
+    stats = zonal_stats(polygon, arr, affine=affine, geojson_out=True)
+    assert 'properties' in stats[0]
+    for key in ['count', 'min', 'max', 'mean']:
+        assert key in stats[0]['properties']
+
+    assert stats[0]['properties']['mean'] == 34
+
+
+# remove when copy_properties alias is removed
+def test_copy_properties_warn():
+    polygons = os.path.join(DATA, 'polygons.shp')
+    # run once to trigger any other unrelated deprecation warnings
+    # so the test does not catch them instead
+    stats_a = zonal_stats(polygons, raster)
+    with pytest.deprecated_call():
+        stats_b = zonal_stats(polygons, raster, copy_properties=True)
+    assert stats_a == stats_b
+    
+
+def test_nan_counts():
+    from affine import Affine
+    transform = Affine(1, 0, 1, 0, -1, 3)
+
+    data = np.array([
+        [np.nan, np.nan, np.nan],
+        [0, 0, 0],
+        [1, 4, 5]
+    ])
+
+    # geom extends an additional row to left
+    geom = 'POLYGON ((1 0, 4 0, 4 3, 1 3, 1 0))'
+
+    # nan stat is requested
+    stats = zonal_stats(geom, data, affine=transform, nodata=0.0, stats="*")
+
+    for res in stats:
+        assert res['count'] == 3  # 3 pixels of valid data
+        assert res['nodata'] == 3  # 3 pixels of nodata
+        assert res['nan'] == 3  # 3 pixels of nans
+
+    # nan are ignored if nan stat is not requested
+    stats = zonal_stats(geom, data, affine=transform, nodata=0.0, stats="count nodata")
+
+    for res in stats:
+        assert res['count'] == 3  # 3 pixels of valid data
+        assert res['nodata'] == 3  # 3 pixels of nodata
+        assert 'nan' not in res
+
+
 def test_percent_cover_zonal_stats():
     polygon = Polygon([[0, 0], [0, 0,5], [1, 1.5], [1.5, 2], [2, 2], [2, 0]])
     arr = np.array([
@@ -480,4 +572,3 @@ def test_geodataframe_zonal():
     expected = zonal_stats(polygons, raster)
     assert zonal_stats(df, raster) == expected
 
-