From 41d4cb2d760bae8fde70d4aea88a98780b05a6e1 Mon Sep 17 00:00:00 2001 From: userw Date: Tue, 28 Mar 2017 13:36:53 -0400 Subject: [PATCH 1/7] Add proof of concept for no_overlap stat (only works for np array rasters) along with test case/data --- src/rasterstats/main.py | 58 ++++++++++++++++-- src/rasterstats/utils.py | 2 +- tests/data/single_polygon_partial_overlap.cpg | 1 + tests/data/single_polygon_partial_overlap.dbf | Bin 0 -> 76 bytes tests/data/single_polygon_partial_overlap.prj | 1 + tests/data/single_polygon_partial_overlap.qpj | 1 + tests/data/single_polygon_partial_overlap.shp | Bin 0 -> 268 bytes tests/data/single_polygon_partial_overlap.shx | Bin 0 -> 108 bytes tests/test_zonal.py | 40 ++++++++++++ 9 files changed, 97 insertions(+), 6 deletions(-) create mode 100644 tests/data/single_polygon_partial_overlap.cpg create mode 100644 tests/data/single_polygon_partial_overlap.dbf create mode 100644 tests/data/single_polygon_partial_overlap.prj create mode 100644 tests/data/single_polygon_partial_overlap.qpj create mode 100644 tests/data/single_polygon_partial_overlap.shp create mode 100644 tests/data/single_polygon_partial_overlap.shx diff --git a/src/rasterstats/main.py b/src/rasterstats/main.py index 8b2b450..fd90b87 100644 --- a/src/rasterstats/main.py +++ b/src/rasterstats/main.py @@ -33,6 +33,7 @@ def gen_zonal_stats( layer=0, band=1, nodata=None, + no_overlap=None, affine=None, stats=None, all_touched=False, @@ -139,7 +140,26 @@ def gen_zonal_stats( warnings.warn("Use `band` to specify band number", DeprecationWarning) band = band_num - with Raster(raster, affine, nodata, band) as rast: + + if 'no_overlap' in stats: + if no_overlap is None: + no_overlap = -998 + if no_overlap == nodata: + no_overlap = -997 + + warnings.warn("Setting no_overlap to {0}; specify no_overlap " + "explicitly when request the `no_overlap` stat".format( + no_overlap)) + elif no_overlap == nodata: + raise Exception("`no_overlap` value is equal to `nodata` value. " + "Values must be distinct to calculate `no_overlap`") + + tmp_nodata = no_overlap if 'no_overlap' in stats else nodata + + + + with Raster(raster, affine=affine, nodata=tmp_nodata, band=band) as rast: + features_iter = read_features(vectors, layer) for _, feat in enumerate(features_iter): geom = shape(feat['geometry']) @@ -149,7 +169,14 @@ def gen_zonal_stats( geom_bounds = tuple(geom.bounds) - fsrc = rast.read(bounds=geom_bounds) + fsrc = rast.read(bounds=geom_bounds, masked=False) + + print stats + print "no_overlap: {0}".format(no_overlap) + print "nodata: {0}".format(nodata) + print "tmp_nodata: {0}".format(tmp_nodata) + print "fsrc.nodata: {0}".format(fsrc.nodata) + # rasterized geometry rv_array = rasterize_geom(geom, like=fsrc, all_touched=all_touched) @@ -157,17 +184,32 @@ def gen_zonal_stats( # nodata mask isnodata = (fsrc.array == fsrc.nodata) + # include actual nodata val when no_overlap is used + if nodata is not None: + isnodata = (isnodata | (fsrc.array == nodata)) + # add nan mask (if necessary) if np.issubdtype(fsrc.array.dtype, float) and \ np.isnan(fsrc.array.min()): isnodata = (isnodata | np.isnan(fsrc.array)) + # Mask the source data array # mask everything that is not a valid value or not within our geom masked = np.ma.MaskedArray( fsrc.array, mask=(isnodata | ~rv_array)) + + print "fsrc.array" + print fsrc.array + print "isnodata" + print isnodata + print "rv_array" + print rv_array + print "masked" + print masked + # execute zone_func on masked zone ndarray if zone_func is not None: if not callable(zone_func): @@ -232,9 +274,15 @@ 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 'no_overlap' in stats: + featmasked = np.ma.MaskedArray(fsrc.array, mask=(~rv_array)) + + if 'nodata' in stats: + feature_stats['nodata'] = float((featmasked == fsrc.nodata).sum()) + if 'no_overlap' in stats: + feature_stats['no_overlap'] = float((featmasked == no_overlap).sum()) + if add_stats is not None: for stat_name, stat_func in add_stats.items(): diff --git a/src/rasterstats/utils.py b/src/rasterstats/utils.py index 9987789..129b67a 100644 --- a/src/rasterstats/utils.py +++ b/src/rasterstats/utils.py @@ -9,7 +9,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', 'no_overlap'] # also percentile_{q} but that is handled as special case diff --git a/tests/data/single_polygon_partial_overlap.cpg b/tests/data/single_polygon_partial_overlap.cpg new file mode 100644 index 0000000..3ad133c --- /dev/null +++ b/tests/data/single_polygon_partial_overlap.cpg @@ -0,0 +1 @@ +UTF-8 \ No newline at end of file diff --git a/tests/data/single_polygon_partial_overlap.dbf b/tests/data/single_polygon_partial_overlap.dbf new file mode 100644 index 0000000000000000000000000000000000000000..063d6f7e8d5edc0030a8a21f21205408eeaabc54 GIT binary patch literal 76 lcmZRMXP07RU|?`$;0BVIATtFn<_BVN!MP9yuL2wx0|0Tv19kua literal 0 HcmV?d00001 diff --git a/tests/data/single_polygon_partial_overlap.prj b/tests/data/single_polygon_partial_overlap.prj new file mode 100644 index 0000000..f2e7a2c --- /dev/null +++ b/tests/data/single_polygon_partial_overlap.prj @@ -0,0 +1 @@ +PROJCS["Albers",GEOGCS["GCS_GRS 1980(IUGG, 1980)",DATUM["D_unknown",SPHEROID["GRS80",6378137,298.257222101]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Albers"],PARAMETER["standard_parallel_1",43],PARAMETER["standard_parallel_2",48],PARAMETER["latitude_of_origin",34],PARAMETER["central_meridian",-120],PARAMETER["false_easting",600000],PARAMETER["false_northing",0],UNIT["Meter",1]] \ No newline at end of file diff --git a/tests/data/single_polygon_partial_overlap.qpj b/tests/data/single_polygon_partial_overlap.qpj new file mode 100644 index 0000000..a656ae9 --- /dev/null +++ b/tests/data/single_polygon_partial_overlap.qpj @@ -0,0 +1 @@ +PROJCS["unnamed",GEOGCS["GRS 1980(IUGG, 1980)",DATUM["unknown",SPHEROID["GRS80",6378137,298.257222101]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]],PROJECTION["Albers_Conic_Equal_Area"],PARAMETER["standard_parallel_1",43],PARAMETER["standard_parallel_2",48],PARAMETER["latitude_of_center",34],PARAMETER["longitude_of_center",-120],PARAMETER["false_easting",600000],PARAMETER["false_northing",0],UNIT["Meter",1]] diff --git a/tests/data/single_polygon_partial_overlap.shp b/tests/data/single_polygon_partial_overlap.shp new file mode 100644 index 0000000000000000000000000000000000000000..7b7268f160c6b320f62d4b037b32b5575dcd23fe GIT binary patch literal 268 zcmZQzQ0HR64%%KYGcd3M<;-L&gO(rTbxb|svc$Ya&+$M-!)~8DypGN7jBkY7^c=Bk zM-gQN@*F^B^-a3x{rn=Y;{@xlrA@7Rj^XJAT@$bJI*JIDY~}@;1yTr72Ldp)*XwmZ z@c`w})y|E)70?UR3p2|pRBiVcpxNO)nJ&#hv(fd#{E@7_{+}RFo?GbE29LA6j$5{> J-xLIz4*=sIOS=F7 literal 0 HcmV?d00001 diff --git a/tests/data/single_polygon_partial_overlap.shx b/tests/data/single_polygon_partial_overlap.shx new file mode 100644 index 0000000000000000000000000000000000000000..412eebf61bba29e0836267187e09fe915e82baac GIT binary patch literal 108 zcmZQzQ0HR64$NLKGcd3M<;-L&gO(rTbxb|svc$Ya&+$M-!)~8DypGN7jBkY7^c=Bk KM-epw@*DtBI}XYK literal 0 HcmV?d00001 diff --git a/tests/test_zonal.py b/tests/test_zonal.py index c91b207..1556e6d 100644 --- a/tests/test_zonal.py +++ b/tests/test_zonal.py @@ -424,6 +424,46 @@ def test_geojson_out(): assert 'id' in feature['properties'] # from orig assert 'count' in feature['properties'] # from zonal stats + +def test_array_overlap_counts(): + from affine import Affine + + nodata = -9999 + no_overlap = -8888 + + transform = Affine(1, 0, 1, 0, -1, 3) + + data = np.array([ + [-9999, 0.0, 516.2840576171875, 524.4825439453125], + [-9999, 178.74169921875, 573.80126953125, 415.345947265625], + [-9999, 397.3150939941406, 568.3016357421875, 185.3491973876953] + ]) + + geom = 'POLYGON ((0 0, 5 0, 5 3, 0 3, 0 0))' + + stats = zonal_stats(geom, data, affine=transform, stats="*", nodata=nodata, no_overlap=no_overlap) + + for res in stats: + assert res['count'] == 9 # Two pixels of valid data + assert res['nodata'] == 3 # Two pixels of nodata + assert res['no_overlap'] == 3 # Three pixels of no overlap + + +def test_raster_overlap_counts() + nodata = -9999 + no_overlap = -8888 + + # same shape/overlap/nodata-pixel as test_array_overlap_counts + polygons = os.path.join(DATA, 'single_polygon_partial_overlap.shp') + + stats = zonal_stats(polygons, raster, stats="*", nodata=nodata, no_overlap=no_overlap) + + for res in stats: + assert res['count'] == 9 # Two pixels of valid data + assert res['nodata'] == 3 # Two pixels of nodata + assert res['no_overlap'] == 3 # Three pixels of no overlap + + # Optional tests def test_geodataframe_zonal(): polygons = os.path.join(DATA, 'polygons.shp') From 0fdc6eb23802d667c753452c6489be41b519021f Mon Sep 17 00:00:00 2001 From: userw Date: Tue, 28 Mar 2017 13:47:08 -0400 Subject: [PATCH 2/7] Fix syntax --- tests/test_zonal.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_zonal.py b/tests/test_zonal.py index d758842..6a7e918 100644 --- a/tests/test_zonal.py +++ b/tests/test_zonal.py @@ -482,7 +482,7 @@ def test_array_overlap_counts(): assert res['no_overlap'] == 3 # Three pixels of no overlap -def test_raster_overlap_counts() +def test_raster_overlap_counts(): nodata = -9999 no_overlap = -8888 From ac66944a9de41f3d15521f3f4a29f38a4cb21d28 Mon Sep 17 00:00:00 2001 From: userw Date: Tue, 28 Mar 2017 14:35:07 -0400 Subject: [PATCH 3/7] Clean up test order --- tests/test_zonal.py | 78 +++++++++++++++++++++------------------------ 1 file changed, 37 insertions(+), 41 deletions(-) diff --git a/tests/test_zonal.py b/tests/test_zonal.py index 6a7e918..9c5442a 100644 --- a/tests/test_zonal.py +++ b/tests/test_zonal.py @@ -457,46 +457,6 @@ def test_geojson_out(): assert 'id' in feature['properties'] # from orig assert 'count' in feature['properties'] # from zonal stats - -def test_array_overlap_counts(): - from affine import Affine - - nodata = -9999 - no_overlap = -8888 - - transform = Affine(1, 0, 1, 0, -1, 3) - - data = np.array([ - [-9999, 0.0, 516.2840576171875, 524.4825439453125], - [-9999, 178.74169921875, 573.80126953125, 415.345947265625], - [-9999, 397.3150939941406, 568.3016357421875, 185.3491973876953] - ]) - - geom = 'POLYGON ((0 0, 5 0, 5 3, 0 3, 0 0))' - - stats = zonal_stats(geom, data, affine=transform, stats="*", nodata=nodata, no_overlap=no_overlap) - - for res in stats: - assert res['count'] == 9 # Two pixels of valid data - assert res['nodata'] == 3 # Two pixels of nodata - assert res['no_overlap'] == 3 # Three pixels of no overlap - - -def test_raster_overlap_counts(): - nodata = -9999 - no_overlap = -8888 - - # same shape/overlap/nodata-pixel as test_array_overlap_counts - polygons = os.path.join(DATA, 'single_polygon_partial_overlap.shp') - - stats = zonal_stats(polygons, raster, stats="*", nodata=nodata, no_overlap=no_overlap) - - for res in stats: - assert res['count'] == 9 # Two pixels of valid data - assert res['nodata'] == 3 # Two pixels of nodata - assert res['no_overlap'] == 3 # Three pixels of no overlap - - # 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 @@ -529,7 +489,6 @@ def test_copy_properties_warn(): def test_nan_counts(): - from affine import Affine transform = Affine(1, 0, 1, 0, -1, 3) data = np.array([ @@ -558,6 +517,43 @@ def test_nan_counts(): assert 'nan' not in res +def test_array_overlap_counts(): + nodata = -9999 + no_overlap = -8888 + + transform = Affine(1, 0, 1, 0, -1, 3) + + data = np.array([ + [-9999, 0.0, 516.2840576171875, 524.4825439453125], + [-9999, 178.74169921875, 573.80126953125, 415.345947265625], + [-9999, 397.3150939941406, 568.3016357421875, 185.3491973876953] + ]) + + geom = 'POLYGON ((0 0, 5 0, 5 3, 0 3, 0 0))' + + stats = zonal_stats(geom, data, affine=transform, stats="*", nodata=nodata, no_overlap=no_overlap) + + for res in stats: + assert res['count'] == 9 # Two pixels of valid data + assert res['nodata'] == 3 # Two pixels of nodata + assert res['no_overlap'] == 3 # Three pixels of no overlap + + +def test_raster_overlap_counts(): + nodata = -9999 + no_overlap = -8888 + + # same shape/overlap/nodata-pixel as test_array_overlap_counts + polygons = os.path.join(DATA, 'single_polygon_partial_overlap.shp') + + stats = zonal_stats(polygons, raster, stats="*", nodata=nodata, no_overlap=no_overlap) + + for res in stats: + assert res['count'] == 9 # Two pixels of valid data + assert res['nodata'] == 3 # Two pixels of nodata + assert res['no_overlap'] == 3 # Three pixels of no overlap + + # Optional tests def test_geodataframe_zonal(): polygons = os.path.join(DATA, 'polygons.shp') From 13cdc84741a0735aab53b8a115638a27ec08c3b0 Mon Sep 17 00:00:00 2001 From: userw Date: Tue, 28 Mar 2017 15:18:54 -0400 Subject: [PATCH 4/7] Add warning to force nodata value when using no_overlap (and code to manage related pieces) and add/update tests --- src/rasterstats/main.py | 44 +++++++++++++++++++++---------------- tests/test_zonal.py | 48 +++++++++++++++++++++++++++++++++++------ 2 files changed, 67 insertions(+), 25 deletions(-) diff --git a/src/rasterstats/main.py b/src/rasterstats/main.py index b0e472a..adb0442 100644 --- a/src/rasterstats/main.py +++ b/src/rasterstats/main.py @@ -142,6 +142,12 @@ def gen_zonal_stats( if 'no_overlap' in stats: + + if 'nodata' in stats and nodata is None: + nodata = -999 + warnings.warn("Setting nodata to -999; specify nodata explicitly " + "when requesting no_overlap stat") + if no_overlap is None: no_overlap = -998 if no_overlap == nodata: @@ -157,7 +163,6 @@ def gen_zonal_stats( tmp_nodata = no_overlap if 'no_overlap' in stats else nodata - with Raster(raster, affine=affine, nodata=tmp_nodata, band=band) as rast: features_iter = read_features(vectors, layer) @@ -171,21 +176,17 @@ def gen_zonal_stats( fsrc = rast.read(bounds=geom_bounds, masked=False) - print stats - print "no_overlap: {0}".format(no_overlap) - print "nodata: {0}".format(nodata) - print "tmp_nodata: {0}".format(tmp_nodata) - print "fsrc.nodata: {0}".format(fsrc.nodata) - - # rasterized geometry rv_array = rasterize_geom(geom, like=fsrc, all_touched=all_touched) + if nodata is None and no_overlap is None: + nodata = fsrc.nodata + # nodata mask isnodata = (fsrc.array == fsrc.nodata) # include actual nodata val when no_overlap is used - if nodata is not None: + if nodata is not None and no_overlap is not None: isnodata = (isnodata | (fsrc.array == nodata)) # add nan mask (if necessary) @@ -201,15 +202,20 @@ def gen_zonal_stats( fsrc.array, mask=(isnodata | ~rv_array)) - - print "fsrc.array" - print fsrc.array - print "isnodata" - print isnodata - print "rv_array" - print rv_array - print "masked" - print masked + # print stats + # print "no_overlap: {0}".format(no_overlap) + # print "nodata: {0}".format(nodata) + # print "tmp_nodata: {0}".format(tmp_nodata) + # print "fsrc.nodata: {0}".format(fsrc.nodata) + + # print "fsrc.array" + # print fsrc.array + # print "isnodata" + # print isnodata + # print "rv_array" + # print rv_array + # print "masked" + # print masked # execute zone_func on masked zone ndarray if zone_func is not None: @@ -280,7 +286,7 @@ def gen_zonal_stats( featmasked = np.ma.MaskedArray(fsrc.array, mask=(~rv_array)) if 'nodata' in stats: - feature_stats['nodata'] = float((featmasked == fsrc.nodata).sum()) + feature_stats['nodata'] = float((featmasked == nodata).sum()) if 'nan' in stats: feature_stats['nan'] = float(np.isnan(featmasked).sum()) if has_nan else 0 if 'no_overlap' in stats: diff --git a/tests/test_zonal.py b/tests/test_zonal.py index 9c5442a..fcf5266 100644 --- a/tests/test_zonal.py +++ b/tests/test_zonal.py @@ -501,7 +501,7 @@ def test_nan_counts(): 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="*") + stats = zonal_stats(geom, data, affine=transform, nodata=0.0, stats="count nodata nan") for res in stats: assert res['count'] == 3 # 3 pixels of valid data @@ -516,6 +516,15 @@ def test_nan_counts(): assert res['nodata'] == 3 # 3 pixels of nodata assert 'nan' not in res + # nan stat is not impacted by no_overlap + 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 + assert res['no_overlap'] == 0 # 3 pixels of nans + def test_array_overlap_counts(): nodata = -9999 @@ -539,19 +548,46 @@ def test_array_overlap_counts(): assert res['no_overlap'] == 3 # Three pixels of no overlap -def test_raster_overlap_counts(): +def test_array_overlap_counts_with_nan(): nodata = -9999 no_overlap = -8888 - # same shape/overlap/nodata-pixel as test_array_overlap_counts - polygons = os.path.join(DATA, 'single_polygon_partial_overlap.shp') + transform = Affine(1, 0, 1, 0, -1, 3) - stats = zonal_stats(polygons, raster, stats="*", nodata=nodata, no_overlap=no_overlap) + data = np.array([ + [-9999, 0.0, 516.2840576171875, np.nan], + [-9999, 178.74169921875, 573.80126953125, 415.345947265625], + [-9999, 397.3150939941406, 568.3016357421875, 185.3491973876953] + ]) + + geom = 'POLYGON ((0 0, 5 0, 5 3, 0 3, 0 0))' + + stats = zonal_stats(geom, data, affine=transform, stats="*", nodata=nodata, no_overlap=no_overlap) for res in stats: - assert res['count'] == 9 # Two pixels of valid data + assert res['count'] == 8 # Two pixels of valid data assert res['nodata'] == 3 # Two pixels of nodata assert res['no_overlap'] == 3 # Three pixels of no overlap + assert res['nan'] == 1 # One pixel of nan + + +# def test_raster_overlap_counts(): +# nodata = -9999 +# no_overlap = -8888 + +# # same shape/overlap/nodata-pixel as test_array_overlap_counts +# polygons = os.path.join(DATA, 'single_polygon_partial_overlap.shp') + +# stats = zonal_stats(polygons, raster, stats="*", nodata=nodata, no_overlap=no_overlap) + +# for res in stats: +# assert res['count'] == 9 # Two pixels of valid data +# assert res['nodata'] == 3 # Two pixels of nodata +# assert res['no_overlap'] == 3 # Three pixels of no overlap + + +# def test_raster_overlap_counts_with_nan(): +# pass # Optional tests From 9e3d2ef4783dcf61d650763cdb88ce887b3c1fba Mon Sep 17 00:00:00 2001 From: userw Date: Tue, 28 Mar 2017 15:54:47 -0400 Subject: [PATCH 5/7] Update no_overlap tests --- tests/test_zonal.py | 48 +++++++++++++-------------------------------- 1 file changed, 14 insertions(+), 34 deletions(-) diff --git a/tests/test_zonal.py b/tests/test_zonal.py index fcf5266..35506e3 100644 --- a/tests/test_zonal.py +++ b/tests/test_zonal.py @@ -532,8 +532,14 @@ def test_array_overlap_counts(): transform = Affine(1, 0, 1, 0, -1, 3) + # data = np.array([ + # [-9999, 0.0, 516.2840576171875, 524.4825439453125], + # [-9999, 178.74169921875, 573.80126953125, 415.345947265625], + # [-9999, 397.3150939941406, 568.3016357421875, 185.3491973876953] + # ]) + data = np.array([ - [-9999, 0.0, 516.2840576171875, 524.4825439453125], + [-9999, 0.0, 516.2840576171875, np.nan], [-9999, 178.74169921875, 573.80126953125, 415.345947265625], [-9999, 397.3150939941406, 568.3016357421875, 185.3491973876953] ]) @@ -543,51 +549,25 @@ def test_array_overlap_counts(): stats = zonal_stats(geom, data, affine=transform, stats="*", nodata=nodata, no_overlap=no_overlap) for res in stats: - assert res['count'] == 9 # Two pixels of valid data + assert res['count'] == 8 # Two pixels of valid data assert res['nodata'] == 3 # Two pixels of nodata assert res['no_overlap'] == 3 # Three pixels of no overlap + assert res['nan'] == 1 # One pixel of nan -def test_array_overlap_counts_with_nan(): +def test_raster_overlap_counts(): nodata = -9999 no_overlap = -8888 - transform = Affine(1, 0, 1, 0, -1, 3) + # same shape/overlap/nodata-pixel as test_array_overlap_counts + polygons = os.path.join(DATA, 'single_polygon_partial_overlap.shp') - data = np.array([ - [-9999, 0.0, 516.2840576171875, np.nan], - [-9999, 178.74169921875, 573.80126953125, 415.345947265625], - [-9999, 397.3150939941406, 568.3016357421875, 185.3491973876953] - ]) - - geom = 'POLYGON ((0 0, 5 0, 5 3, 0 3, 0 0))' - - stats = zonal_stats(geom, data, affine=transform, stats="*", nodata=nodata, no_overlap=no_overlap) + stats = zonal_stats(polygons, raster, stats="*", nodata=nodata, no_overlap=no_overlap) for res in stats: - assert res['count'] == 8 # Two pixels of valid data + assert res['count'] == 9 # Two pixels of valid data assert res['nodata'] == 3 # Two pixels of nodata assert res['no_overlap'] == 3 # Three pixels of no overlap - assert res['nan'] == 1 # One pixel of nan - - -# def test_raster_overlap_counts(): -# nodata = -9999 -# no_overlap = -8888 - -# # same shape/overlap/nodata-pixel as test_array_overlap_counts -# polygons = os.path.join(DATA, 'single_polygon_partial_overlap.shp') - -# stats = zonal_stats(polygons, raster, stats="*", nodata=nodata, no_overlap=no_overlap) - -# for res in stats: -# assert res['count'] == 9 # Two pixels of valid data -# assert res['nodata'] == 3 # Two pixels of nodata -# assert res['no_overlap'] == 3 # Three pixels of no overlap - - -# def test_raster_overlap_counts_with_nan(): -# pass # Optional tests From 117390dbc76ff9580fc43f74393386d3a2ab00ac Mon Sep 17 00:00:00 2001 From: userw Date: Tue, 28 Mar 2017 16:05:37 -0400 Subject: [PATCH 6/7] Add additional no_overlap test --- tests/test_zonal.py | 23 +++++++++++++++++++++++ 1 file changed, 23 insertions(+) diff --git a/tests/test_zonal.py b/tests/test_zonal.py index 35506e3..69c366e 100644 --- a/tests/test_zonal.py +++ b/tests/test_zonal.py @@ -555,6 +555,29 @@ def test_array_overlap_counts(): assert res['nan'] == 1 # One pixel of nan +def test_missing_no_overlap_logic(): + nodata = -998 + no_overlap = None + + transform = Affine(1, 0, 1, 0, -1, 3) + + data = np.array([ + [-998, 0.0, 516.2840576171875, np.nan], + [-998, 178.74169921875, 573.80126953125, 415.345947265625], + [-998, 397.3150939941406, 568.3016357421875, 185.3491973876953] + ]) + + geom = 'POLYGON ((0 0, 5 0, 5 3, 0 3, 0 0))' + + stats = zonal_stats(geom, data, affine=transform, stats="*", nodata=nodata, no_overlap=no_overlap) + + for res in stats: + assert res['count'] == 8 # Two pixels of valid data + assert res['nodata'] == 3 # Two pixels of nodata + assert res['no_overlap'] == 3 # Three pixels of no overlap + assert res['nan'] == 1 # One pixel of nan + + def test_raster_overlap_counts(): nodata = -9999 no_overlap = -8888 From 8f35c27465c824b4eb71e8fc1f6d6390e0ffbada Mon Sep 17 00:00:00 2001 From: userw Date: Tue, 28 Mar 2017 16:13:02 -0400 Subject: [PATCH 7/7] Add additional no_overlap test --- tests/test_zonal.py | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) diff --git a/tests/test_zonal.py b/tests/test_zonal.py index 69c366e..5c52556 100644 --- a/tests/test_zonal.py +++ b/tests/test_zonal.py @@ -578,6 +578,24 @@ def test_missing_no_overlap_logic(): assert res['nan'] == 1 # One pixel of nan +def test_nodata_and_no_overlap_check(): + nodata = -9999 + no_overlap = -9999 + + transform = Affine(1, 0, 1, 0, -1, 3) + + data = np.array([ + [-998, 0.0, 516.2840576171875, np.nan], + [-998, 178.74169921875, 573.80126953125, 415.345947265625], + [-998, 397.3150939941406, 568.3016357421875, 185.3491973876953] + ]) + + geom = 'POLYGON ((0 0, 5 0, 5 3, 0 3, 0 0))' + + with pytest.raises(Exception): + stats = zonal_stats(geom, data, affine=transform, stats="*", nodata=nodata, no_overlap=no_overlap) + + def test_raster_overlap_counts(): nodata = -9999 no_overlap = -8888