From 95609c3770bcd0b43a000f2a1293b86748755e73 Mon Sep 17 00:00:00 2001
From: vagrant <vagrant@precise64.(none)>
Date: Wed, 18 Sep 2013 05:44:06 +0000
Subject: [PATCH] initial commit

---
 .gitignore                   |  12 ++++
 README.md                    |  63 ++++++++++++++++++++
 scripts/coverage.sh          |   1 +
 scripts/provision.sh         |  16 +++++
 scripts/rasterstats          |  47 +++++++++++++++
 setup.cfg                    |   5 ++
 setup.py                     |  49 +++++++++++++++
 src/rasterstats/__init__.py  |   5 ++
 src/rasterstats/main.py      | 112 +++++++++++++++++++++++++++++++++++
 src/rasterstats/utils.py     |  16 +++++
 tests/__init__.py            |   0
 tests/conftest.py            |   3 +
 tests/data/lines.dbf         | Bin 0 -> 87 bytes
 tests/data/lines.prj         |  16 +++++
 tests/data/lines.qpj         |   1 +
 tests/data/lines.shp         | Bin 0 -> 436 bytes
 tests/data/lines.shx         | Bin 0 -> 116 bytes
 tests/data/points.dbf        | Bin 0 -> 98 bytes
 tests/data/points.prj        |  16 +++++
 tests/data/points.qpj        |   1 +
 tests/data/points.shp        | Bin 0 -> 184 bytes
 tests/data/points.shx        | Bin 0 -> 124 bytes
 tests/data/points_noproj.dbf | Bin 0 -> 98 bytes
 tests/data/points_noproj.shp | Bin 0 -> 184 bytes
 tests/data/points_noproj.shx | Bin 0 -> 124 bytes
 tests/data/polygons.dbf      | Bin 0 -> 87 bytes
 tests/data/polygons.prj      |  16 +++++
 tests/data/polygons.qpj      |   1 +
 tests/data/polygons.shp      | Bin 0 -> 436 bytes
 tests/data/polygons.shx      | Bin 0 -> 116 bytes
 tests/data/slope.tif         | Bin 0 -> 27542 bytes
 tests/test_zonal.py          |  30 ++++++++++
 32 files changed, 410 insertions(+)
 create mode 100644 .gitignore
 create mode 100644 README.md
 create mode 100644 scripts/coverage.sh
 create mode 100644 scripts/provision.sh
 create mode 100644 scripts/rasterstats
 create mode 100644 setup.cfg
 create mode 100644 setup.py
 create mode 100644 src/rasterstats/__init__.py
 create mode 100644 src/rasterstats/main.py
 create mode 100644 src/rasterstats/utils.py
 create mode 100644 tests/__init__.py
 create mode 100644 tests/conftest.py
 create mode 100644 tests/data/lines.dbf
 create mode 100644 tests/data/lines.prj
 create mode 100644 tests/data/lines.qpj
 create mode 100644 tests/data/lines.shp
 create mode 100644 tests/data/lines.shx
 create mode 100644 tests/data/points.dbf
 create mode 100644 tests/data/points.prj
 create mode 100644 tests/data/points.qpj
 create mode 100644 tests/data/points.shp
 create mode 100644 tests/data/points.shx
 create mode 100644 tests/data/points_noproj.dbf
 create mode 100644 tests/data/points_noproj.shp
 create mode 100644 tests/data/points_noproj.shx
 create mode 100644 tests/data/polygons.dbf
 create mode 100644 tests/data/polygons.prj
 create mode 100644 tests/data/polygons.qpj
 create mode 100644 tests/data/polygons.shp
 create mode 100644 tests/data/polygons.shx
 create mode 100644 tests/data/slope.tif
 create mode 100644 tests/test_zonal.py

diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000..4223611
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,12 @@
+*.pyc
+.coverage
+TEST*.xml
+build/*
+dist/*
+*.swp
+.DS_Store
+*.orig
+*.egg-info
+MANIFEST
+Vagrantfile
+.vagrant/*
diff --git a/README.md b/README.md
new file mode 100644
index 0000000..86df251
--- /dev/null
+++ b/README.md
@@ -0,0 +1,63 @@
+# python-raster-stats
+
+Summary statistics of raster dataset values based on vector geometries.
+
+Docs (API, Topics)
+Examples
+Build status
+Test coverage
+
+## Quickstart
+
+You've got a raster dataset representing elevation and a vector dataset representing county boundaries. 
+Show the average, min and max elevation for each county
+Show percentage of land cover by county
+Output to csv (via CLI)
+Pass geometries via wkt or wkb (single & list)
+Integrate with other python packages via __geo_interface__
+
+## Features
+
+* Raster data support: 
+  * Any continuous raster band supported by GDAL
+* Vector data support:
+  * OGR layer
+* Python module (returns built in python data structures - list of dicts)
+* Depends on GDAL, GEOS, shapely and numpy
+* Full coverage unit testing
+
+## Issues
+To report a bug via github issues: provide smallest possible raster, vector and code to reproduce it
+
+## Docs
+
+## TODO 
+* respects null/no-data values
+* covering edge cases for input datasets
+* command line interface which returns csv data and optionally copies over original vector attrs
+* supports categorical
+* Vector data:
+  * Points, Lines, Polygon and Multi-* geometries
+  * Can specify:
+    * OGR layer
+    * single geoms or list of geometries represented by wkts, wkbs or any object that supports the __geo_interface__
+* projection support
+* pip installable
+* python 2 & 3 support
+* buildthedocs OR use some sort of literate programming
+* travis-ci and https://coveralls.io/info/features
+* Examples for PyShp, GeoDjango, Fiona, Path to OGR resource, ArcPy (as Ipython notebooks?)
+* reproject on the fly using projection metadata
+* heavily profiled using a wide range of input data. The resulting heuristics used to automatically configure for optimal performance. Optimzation heuristic for determining global_src_extent - number of features - extent of features - available system memory vs raster_extent
+* CLI: pivots for categorical
+* support parallel processing on multiple CPUs via the `multiprocessing` approach
+* zonal majority (http://stackoverflow.com/questions/6252280/find-the-most-frequent-number-in-a-numpy-vector)
+
+## Alternatives
+
+Grass r.stats
+R spatialdataframe
+starspan
+zonal statistics arcpy
+QGIS
+
diff --git a/scripts/coverage.sh b/scripts/coverage.sh
new file mode 100644
index 0000000..9d15f67
--- /dev/null
+++ b/scripts/coverage.sh
@@ -0,0 +1 @@
+coverage run setup.py test && coverage report -m --include=*rasterstats*
diff --git a/scripts/provision.sh b/scripts/provision.sh
new file mode 100644
index 0000000..da8db18
--- /dev/null
+++ b/scripts/provision.sh
@@ -0,0 +1,16 @@
+#!/bin/bash
+# Provision base software required for running raster_stats
+
+#apt-get install -y python-software-properties
+#add-apt-repository -y ppa:ubuntugis/ubuntugis-unstable
+
+apt-get update
+
+apt-get install -y libgdal-dev gdal-bin \
+                   python-gdal python-pip python-numpy \
+                   libspatialindex-dev libspatialindex1 \
+                   build-essential git atop python-dev python-dateutil
+
+cd /usr/local/src/python-raster-stats
+sudo pip install shapely pytest coverage
+python setup.py develop
diff --git a/scripts/rasterstats b/scripts/rasterstats
new file mode 100644
index 0000000..ac941a1
--- /dev/null
+++ b/scripts/rasterstats
@@ -0,0 +1,47 @@
+#!/bin/env python
+"""
+rasterstats
+
+Usage:
+  rasterstats COMMAND DATASET [ARGS ...]
+  rasterstats (-h | --help)
+  rasterstats --version
+
+Options:
+  -h --help     Show this screen.
+  --version     Show version.
+"""
+import sys
+from docopt import docopt
+from rasterstats import raster_stats
+
+
+def guess_type(thing):
+    try:
+        return float(thing)
+    except:
+        pass
+    return str(thing)
+
+
+if __name__ == '__main__':
+    opt = docopt(__doc__, version='geofu 0.0.1')
+    if opt['DATASET'] == "-":
+        # use stdin
+        opt['DATASET'] = sys.stdin.readlines()[-1].strip()
+    lyr = load(opt['DATASET'])
+    args = []
+    kwargs = {}
+    for arg in opt['ARGS']:
+        if "=" in arg:
+            k, v = arg.split("=")
+            v = guess_type(v)
+            kwargs[k] = v
+        else:
+            arg = guess_type(arg)
+            args.append(arg)
+    res = getattr(lyr, opt['COMMAND'])(*args, **kwargs)
+    if isinstance(res, lyr.__class__):
+        print res.path
+    else:
+        print res
diff --git a/setup.cfg b/setup.cfg
new file mode 100644
index 0000000..875240f
--- /dev/null
+++ b/setup.cfg
@@ -0,0 +1,5 @@
+# content of setup.cfg
+[pytest]
+norecursedirs = examples* src* scripts*
+addopts = -rf -v
+#addopts = --maxfail=1 --pdb -rf -v
diff --git a/setup.py b/setup.py
new file mode 100644
index 0000000..de722db
--- /dev/null
+++ b/setup.py
@@ -0,0 +1,49 @@
+import os
+import sys
+from setuptools import setup
+from setuptools.command.test import test as TestCommand
+
+
+def read(fname):
+    return open(os.path.join(os.path.dirname(__file__), fname)).read()
+
+
+class PyTest(TestCommand):
+    def finalize_options(self):
+        TestCommand.finalize_options(self)
+        self.test_args = []
+        self.test_suite = True
+
+    def run_tests(self):
+        import pytest
+        errno = pytest.main(self.test_args)
+        sys.exit(errno)
+
+
+setup(
+    name="rasterstats",
+    version="0.1a",
+    author="Matthew Perry",
+    author_email="perrygeo@gmail.com",
+    description=("Summarize geospatial raster datasets based on vector geometries"),
+    license="BSD",
+    keywords="gis geospatial geographic raster vector zonal statistics",
+    url="https://github.com/perrygeo/python-raster-stats",
+    package_dir={'': 'src'},
+    packages=['rasterstats'],
+    long_description=read('README.md'),
+    install_requires=[
+        'shapely',
+        'numpy',
+        'GDAL',
+        # pandas, pyproj?
+    ],
+    scripts=['scripts/rasterstats'],
+    tests_require=['pytest'],
+    cmdclass = {'test': PyTest},
+    classifiers=[
+        "Development Status :: 1 - Planning",
+        "Topic :: Utilities",
+        "License :: OSI Approved :: BSD License",
+    ],
+)
diff --git a/src/rasterstats/__init__.py b/src/rasterstats/__init__.py
new file mode 100644
index 0000000..4abe42d
--- /dev/null
+++ b/src/rasterstats/__init__.py
@@ -0,0 +1,5 @@
+# -*- coding: utf-8 -*-
+from .main import raster_stats
+
+__all__ = ['raster_stats']
+
diff --git a/src/rasterstats/main.py b/src/rasterstats/main.py
new file mode 100644
index 0000000..37045c7
--- /dev/null
+++ b/src/rasterstats/main.py
@@ -0,0 +1,112 @@
+# -*- coding: utf-8 -*-
+from shapely.geometry import mapping, shape
+import numpy as np
+from osgeo import gdal, ogr
+from osgeo.gdalconst import GA_ReadOnly
+from .utils import bbox_to_pixel_offsets
+gdal.PushErrorHandler('CPLQuietErrorHandler')
+
+
+def raster_stats(vector, raster, layer_num=0, band_num=1, nodata_value=None, 
+                 global_src_extent=False):
+    rds = gdal.Open(raster, GA_ReadOnly)
+    assert(rds)
+    rb = rds.GetRasterBand(band_num)
+    rgt = rds.GetGeoTransform()
+
+    if nodata_value is not None:
+        nodata_value = float(nodata_value)
+        rb.SetNoDataValue(nodata_value)
+
+    vds = ogr.Open(vector, GA_ReadOnly)
+    assert(vds)
+    vlyr = vds.GetLayer(layer_num)
+
+    # create an in-memory numpy array of the source raster data
+    # covering the whole extent of the vector layer
+    if global_src_extent:
+        # use global source extent
+        # useful only when disk IO or raster scanning inefficiencies are your limiting factor
+        # advantage: reads raster data in one pass
+        # disadvantage: large vector extents may have big memory requirements
+        src_offset = bbox_to_pixel_offsets(rgt, vlyr.GetExtent())
+        src_array = rb.ReadAsArray(*src_offset)
+
+        # calculate new geotransform of the layer subset
+        new_gt = (
+            (rgt[0] + (src_offset[0] * rgt[1])),
+            rgt[1],
+            0.0,
+            (rgt[3] + (src_offset[1] * rgt[5])),
+            0.0,
+            rgt[5]
+        )
+
+    mem_drv = ogr.GetDriverByName('Memory')
+    driver = gdal.GetDriverByName('MEM')
+
+    #geometries = iter(vlyr, None)
+    #for feat in geometries:
+    stats = []
+
+    # Loop through geometries
+    feat = vlyr.GetNextFeature()
+    while feat:
+        if not global_src_extent:
+            # use local source extent
+            # fastest option when you have fast disks and well indexed raster (ie tiled Geotiff)
+            # advantage: each feature uses the smallest raster chunk
+            # disadvantage: lots of reads on the source raster
+            src_offset = bbox_to_pixel_offsets(rgt, feat.geometry().GetEnvelope())
+            src_array = rb.ReadAsArray(*src_offset)
+
+            # calculate new geotransform of the feature subset
+            new_gt = (
+                (rgt[0] + (src_offset[0] * rgt[1])),
+                rgt[1],
+                0.0,
+                (rgt[3] + (src_offset[1] * rgt[5])),
+                0.0,
+                rgt[5]
+            )
+
+        # Create a temporary vector layer in memory
+        mem_ds = mem_drv.CreateDataSource('out')
+        mem_layer = mem_ds.CreateLayer('poly', None, ogr.wkbPolygon)
+        mem_layer.CreateFeature(feat.Clone())
+
+        # Rasterize it
+        rvds = driver.Create('', src_offset[2], src_offset[3], 1, gdal.GDT_Byte)
+        rvds.SetGeoTransform(new_gt)
+        gdal.RasterizeLayer(rvds, [1], mem_layer, burn_values=[1])
+        rv_array = rvds.ReadAsArray()
+
+        # Mask the source data array with our current feature
+        # we take the logical_not to flip 0<->1 to get the correct mask effect
+        # we also mask out nodata values explictly
+        masked = np.ma.MaskedArray(
+            src_array,
+            mask=np.logical_or(
+                src_array == nodata_value,
+                np.logical_not(rv_array)
+            )
+        )
+
+        feature_stats = {
+            'min': float(masked.min()),
+            'mean': float(masked.mean()),
+            'max': float(masked.max()),
+            'std': float(masked.std()),
+            'sum': float(masked.sum()),
+            'count': int(masked.count()),
+            'fid': int(feat.GetFID())}
+
+        stats.append(feature_stats)
+
+        rvds = None
+        mem_ds = None
+        feat = vlyr.GetNextFeature()
+
+    vds = None
+    rds = None
+    return stats
diff --git a/src/rasterstats/utils.py b/src/rasterstats/utils.py
new file mode 100644
index 0000000..82478c6
--- /dev/null
+++ b/src/rasterstats/utils.py
@@ -0,0 +1,16 @@
+# -*- coding: utf-8 -*-
+
+def bbox_to_pixel_offsets(gt, bbox):
+    originX = gt[0]
+    originY = gt[3]
+    pixel_width = gt[1]
+    pixel_height = gt[5]
+    x1 = int((bbox[0] - originX) / pixel_width)
+    x2 = int((bbox[1] - originX) / pixel_width) + 1
+
+    y1 = int((bbox[3] - originY) / pixel_height)
+    y2 = int((bbox[2] - originY) / pixel_height) + 1
+
+    xsize = x2 - x1
+    ysize = y2 - y1
+    return (x1, y1, xsize, ysize)
diff --git a/tests/__init__.py b/tests/__init__.py
new file mode 100644
index 0000000..e69de29
diff --git a/tests/conftest.py b/tests/conftest.py
new file mode 100644
index 0000000..b554816
--- /dev/null
+++ b/tests/conftest.py
@@ -0,0 +1,3 @@
+import logging
+import sys
+logging.basicConfig(stream=sys.stderr, level=logging.INFO)
diff --git a/tests/data/lines.dbf b/tests/data/lines.dbf
new file mode 100644
index 0000000000000000000000000000000000000000..ac8457820bae7f137af807b7b921c5f600d35bb3
GIT binary patch
literal 87
rcmZRMXP07PU|?`$;0BUtAe@0AGX*Z@2V!x-xex}g0vs5^Sqeq~v>63H

literal 0
HcmV?d00001

diff --git a/tests/data/lines.prj b/tests/data/lines.prj
new file mode 100644
index 0000000..932909e
--- /dev/null
+++ b/tests/data/lines.prj
@@ -0,0 +1,16 @@
+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["metre",1,
+        AUTHORITY["EPSG","9001"]]]
+
diff --git a/tests/data/lines.qpj b/tests/data/lines.qpj
new file mode 100644
index 0000000..5fbc831
--- /dev/null
+++ b/tests/data/lines.qpj
@@ -0,0 +1 @@
+GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]]
diff --git a/tests/data/lines.shp b/tests/data/lines.shp
new file mode 100644
index 0000000000000000000000000000000000000000..dfe374a476e4bb809f2826cc55322afc6ae44eb2
GIT binary patch
literal 436
zcmZQzQ0HR64sN|*W?%p!hOo)3B1hlxI%;f|_xxO^=P0SuadDpjpQA~w?z)^VJxA=?
zQA8PmyabS0@iVsMl(X>RGz+AV6RH(v54ycx%p&u=8uT3Nl)bjH1Kq*Xm)NM^sOLC;
zv$ET6Ha<s&AN<;DTJ#*R9+_XK$-(D%=<4nq)pk9{w8O7{asm0vH~d|a*s13Ta|_HL
zg!n*bF#)jya+pkR+8J#PbPKxOAb+qz{Q(O<nA&D`Z^7IjypB<m8u+_`?ym8?`0~YH
gUdL0<CV#SR)^ogAp~iZHk<YR1O*-oiARpaq05d6kGXMYp

literal 0
HcmV?d00001

diff --git a/tests/data/lines.shx b/tests/data/lines.shx
new file mode 100644
index 0000000000000000000000000000000000000000..a52806c68fb868c5f70c8e47232bb91b79dc8afa
GIT binary patch
literal 116
zcmZQzQ0HR64y;}<GcW)VL)c_ik)!W;9W^$~dw#CdbClHSxVTS%&(Wk-cU?}Go+Ebc
PD56F{UIGwL17ZgNg6j^H

literal 0
HcmV?d00001

diff --git a/tests/data/points.dbf b/tests/data/points.dbf
new file mode 100644
index 0000000000000000000000000000000000000000..3aee2835477f4a0fd80d129f2d3096912182e6a7
GIT binary patch
literal 98
ucmZRMXP07TU|?`$;0BUtAe@0AGX*Z@2V!x-xex}g0vs5^Sqer7rZE8abOy2j

literal 0
HcmV?d00001

diff --git a/tests/data/points.prj b/tests/data/points.prj
new file mode 100644
index 0000000..932909e
--- /dev/null
+++ b/tests/data/points.prj
@@ -0,0 +1,16 @@
+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["metre",1,
+        AUTHORITY["EPSG","9001"]]]
+
diff --git a/tests/data/points.qpj b/tests/data/points.qpj
new file mode 100644
index 0000000..5fbc831
--- /dev/null
+++ b/tests/data/points.qpj
@@ -0,0 +1 @@
+GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]]
diff --git a/tests/data/points.shp b/tests/data/points.shp
new file mode 100644
index 0000000000000000000000000000000000000000..3fc12bd04da047f9310631626bc61d565b09f320
GIT binary patch
literal 184
zcmZQzQ0HR64q{#~GcYg$<)%6AX#e?~*HPDEres{bo+H1~RD+lNe2zxf3qCrx>p5cA
qjv@*&h6~v&B)v>9xtcp7EjPa)>;Z}cT>#Msvy-PUu~ENK&k+EZogq~K

literal 0
HcmV?d00001

diff --git a/tests/data/points.shx b/tests/data/points.shx
new file mode 100644
index 0000000000000000000000000000000000000000..740f6b2e37ca30b67dbfc24acf9537eef68a6220
GIT binary patch
literal 124
zcmZQzQ0HR64(whqGcYg$<)%6AX#e?~*HPDEres{bo+H1~RD+lNe2zxf3qCrx>p5cA
Rjv{IV<bmvQfYN>-8UQHf4#5Bb

literal 0
HcmV?d00001

diff --git a/tests/data/points_noproj.dbf b/tests/data/points_noproj.dbf
new file mode 100644
index 0000000000000000000000000000000000000000..3aee2835477f4a0fd80d129f2d3096912182e6a7
GIT binary patch
literal 98
ucmZRMXP07TU|?`$;0BUtAe@0AGX*Z@2V!x-xex}g0vs5^Sqer7rZE8abOy2j

literal 0
HcmV?d00001

diff --git a/tests/data/points_noproj.shp b/tests/data/points_noproj.shp
new file mode 100644
index 0000000000000000000000000000000000000000..3fc12bd04da047f9310631626bc61d565b09f320
GIT binary patch
literal 184
zcmZQzQ0HR64q{#~GcYg$<)%6AX#e?~*HPDEres{bo+H1~RD+lNe2zxf3qCrx>p5cA
qjv@*&h6~v&B)v>9xtcp7EjPa)>;Z}cT>#Msvy-PUu~ENK&k+EZogq~K

literal 0
HcmV?d00001

diff --git a/tests/data/points_noproj.shx b/tests/data/points_noproj.shx
new file mode 100644
index 0000000000000000000000000000000000000000..740f6b2e37ca30b67dbfc24acf9537eef68a6220
GIT binary patch
literal 124
zcmZQzQ0HR64(whqGcYg$<)%6AX#e?~*HPDEres{bo+H1~RD+lNe2zxf3qCrx>p5cA
Rjv{IV<bmvQfYN>-8UQHf4#5Bb

literal 0
HcmV?d00001

diff --git a/tests/data/polygons.dbf b/tests/data/polygons.dbf
new file mode 100644
index 0000000000000000000000000000000000000000..ac8457820bae7f137af807b7b921c5f600d35bb3
GIT binary patch
literal 87
rcmZRMXP07PU|?`$;0BUtAe@0AGX*Z@2V!x-xex}g0vs5^Sqeq~v>63H

literal 0
HcmV?d00001

diff --git a/tests/data/polygons.prj b/tests/data/polygons.prj
new file mode 100644
index 0000000..932909e
--- /dev/null
+++ b/tests/data/polygons.prj
@@ -0,0 +1,16 @@
+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["metre",1,
+        AUTHORITY["EPSG","9001"]]]
+
diff --git a/tests/data/polygons.qpj b/tests/data/polygons.qpj
new file mode 100644
index 0000000..5fbc831
--- /dev/null
+++ b/tests/data/polygons.qpj
@@ -0,0 +1 @@
+GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]]
diff --git a/tests/data/polygons.shp b/tests/data/polygons.shp
new file mode 100644
index 0000000000000000000000000000000000000000..80d58e8c72ae1f66163096057270420064438672
GIT binary patch
literal 436
zcmZQzQ0HR64sN|*W?*0i%Bfv;n|t~muVb2%L2hTOp5q7RJ-!*de2z8`lP5`b=s9B7
zjv~qk<OLv`b=CRl$K^n?gmh#N>ObRk{3{VGp4G1B2vW!n)e5tB;SJTYB%s_&=0Mw<
zj}U5MdX26Zd~^oNqw^PQW!MM;&4$?pGw%gU?$v!jyJ38o`9N!#VD60B`8~pvoe!s<
zk==RX_WP59Tzrl#37i33!0_<?*frrcHzIss;dNVjuflSmJJI>(cbE7X0QJMd7iPB^
NXHT+Jo1P=O`2aisf+GL`

literal 0
HcmV?d00001

diff --git a/tests/data/polygons.shx b/tests/data/polygons.shx
new file mode 100644
index 0000000000000000000000000000000000000000..a7d2fe5d56ba00a9c6b07cb20b4993820f4a61e4
GIT binary patch
literal 116
zcmZQzQ0HR64y;}<Gcd3M<<u^_%{_gO*D=k>Ah)ws&+!BE9^VXJK1Z8}$&(~I^c=Bk
PM-epw@&bUk4MYO~w?q#s

literal 0
HcmV?d00001

diff --git a/tests/data/slope.tif b/tests/data/slope.tif
new file mode 100644
index 0000000000000000000000000000000000000000..c0be86fef27bf3d72bc2b8bbe7c019be24b4df7a
GIT binary patch
literal 27542
zcmeFZc{G*%|2}G-G7q6r(LfneX)tA9ubU`ElW5R{N~K9@9t@G8G9{S`nMLL?5+P+s
zDwGsbkyN6RIQR4UKHul}eSYVUbJkgBt@F>>Yi;lMzW2S?-uHFChu7=6u9s!YW^)bg
z=Him%;^OAz;^yMwTE{-Q|L$9`eV)JjQ`tW6-~GSN6Zju}DQ-S4KK9sa_9^#2j*tJl
zufQIc|NERewm<2AoU``tJ}-M*++CFGul6(hbn_795@1^a_RG?}c)2#S{e^7*j+ZEx
z@L%V!{R}T2uC0F^=VAZ+-+1Uva_2tismuBA_W$Z5eNB6M?Hk1~bi+JonXeZ9`+V;I
zASK83*G2!;XD{tPEe~IYm(X7)v46O^#?EmIWq7S;U-10>THrtJ<!M5HzZMp{<HgNh
z<+%T}f?Q(%Y1ubL|LgNSf498<X~%K#{imHZO+fUo&v9`H{-@>U68cZe%O(7umWNB^
zKP?~Exc{{L|NXQ7t37YR`LB-I-!1-bx7_>d_<2i}f3A~#$%$WQgT!CG!AiFdbN+h&
zugB5-`?ote*q=6K)=srxds7Xknyy(pb*}DWZM9`qrluPI>CQTBy#2^gC-%*w`}giY
zaBTl+lYRRf4>})kwEwRUEY?1~>ae4e!(XSE9(Hipf6UQt_vzDITn3XW|NrOz%NYFE
zo#VPUR?En>JYp)nm!aZOAinO($JUbvQM+IXj30!c`~7`<w^U+!ujMf=74I3LcU??T
z<r_wPWgD|u=_MmS(86dOX<<$mwlLhATN&}BRwn6AE7Nzal?m5sWsd5%F~4o!GG}hS
zVq!w-m>rd^%qpj9rfEqm6BZTDgxZBMuTMlU1qxA2%Gz+ILpq4r;_Aa}RlC6KYIkL1
z_n%>2+PN@_Qw}gulkA!CMb?aL>uN^BdIl5!Ac+3HjQ{VN{~nuv|AJL1jm(Fqb<CbD
z57<p8z=hA5xWC5_)~B;DMX><}i{2t5whS|}L-4010v01LuxZje*u3w=>&SN0=e1((
z;d*2bM58866yw8VIa0f)B8bl$OApm!$09*G75otn8Yc8v#*KojW-(@ut}^1)X^hZ`
zTg(%mNT%sX1atXRB-3LO!Sr0Y#<(o=XYx7TjHZJ-<M8<uGkWtNV{^-a3D(-oj7+s+
zOkS;K{2LsZeL9Xzv%dvXJKvnK%wm{!=f#XpvKI4OZU(a;OP;wrLx?F!8=|XYA8Gl`
zE^_(ML6(-E>Ckk3W?QK&6BMe<<e6$R!so@At;KHt&+g6thNs@>GlqA56H`-EgHDri
zG_IVVoP-A9=l%sYzVGog>m~MEJwTXPI6lOM!??K=_aq*|dU6ePgsQP^bryD|oWp~)
zvvECyAG5SGI3A&sv9HJtzgFCXhWmKBxnUK#sBI^O`ZXkwD@mjM!gPE;QICo`-AG-B
z?Gj%(mV9HJXSK&ras3NirN82CK@DCCM8jMt0T&7qP*vFpfjl8nn$`%F_FFh6dyHnh
zkzu5Uj2P(xFg8Dpm@AtWF|V4{nVPt%%*ad`W_JJ|lO^(<`WAksuT7uGZT<&3c={b}
zZ~skxYowV^`&5{1vRcdv&N60V_BuxB=L+V8sTK41;08v0cmwl9V+(Wqg(I`E^8n*`
zVkdL#$vVc*+muPuGGg?<FJ^xHSit<VZvH(!McBE?AMlb%Y50!qMgpWej*AX_`HiXT
zd1zM@KS?s&B+%TC(+at$o*sj2$s`;WzYK5bRd8RT3ukqE<m4Yj-odT--MJJ7+T!4u
z>57XvPvEQ}L!RGs$tGtbRrKzqZUZZt6R%47cLb?)mmDc;YLNSvJWS?uL5^V<tmaKZ
zO?eRtr`MqJ(<5}1R^r{)GGy1hV&{z_H8o37vHf!x?Qe!SZ$B*W@{?SdB;AY{M?!%^
z_`uZ*Gr<@5?a+sJxnr>Ic!2JDMeGY^ko0UOe$5ZylvzLK$g|&nDoBGq<&0s|?w=^&
zQ>1fqxtY?H^BFb4MT}1CV&<fq9%JHU%qZ?#%Lu)(VRouoF@_;4m;=HLqquf4bF4;_
z34gBsf5J1)>Mc|MtdkkY7N$+7gsI`&4{ROog}B-`6usf1zQrHG@0^I_<dw*1a6-rU
z05pzVK-!B%cxWhxv2}XTYP*iN-4U2`+!lNqs}ZyFBqSC`Au;4Ba?gB%xv&Dc?PRDX
zd<9J~BdRoAOdVx&>2!cJS?-L7_qg-e)^`mnDpC;>R}PIyS@73R$Lhvvj9qKT&gk*9
zCuuglit2*llT*myx5o6newdh8jGU5w7#y5L8)|gvkEa%8Nr{q_$rG#?DTH)|9iA8o
z;J_XmY+USz`~|w0b@(3VXR14*#!aBRygKA}cOi94ET$m;g|wj&G+o+;GBfPx_`z-D
z(Qt@XSi6y~g$wB(Uq|(gV>qAV2A*jpkXjvpcWHL`-^b_H<`2xM%_k=Ht02X{9LD~D
zN*HGP!QrPX+Je2|e#;FvzBys9?i^SwKMK_czGy4Df;W0@uutC#8Z*H#?+rW~Ov5i_
zcU-ph0e4^^CjE}bw}2+NmGjaT%p%^*DRe%l10G%(D0v);vC|Ro4Ge)_{#888$;QC5
z9PGFwOh&IHXmnx=ynZF(dB8)c4~r20C2d+LJc%@q<zS__H}um_VC^|a%=>x(N4^E)
zK~@DWjB7%wc@buZdqM0;6=(S!E2Li;!O`lev}c<<{rEhQMoR?grb#x6PlRzaE<0hm
z`#b1wk|Ymt6|ydvMtt9=Q1d-CYO2>E&Qo0)Kd_umSQ(OBksLh{AHk$ewQ$&(i!B!K
zkhAJGwkqBHKjGQJ^fSAreq-|2T*c_K6Oh<!h~7vocJ3TQ$hS){jrM_S=rx$PU%_>W
zQ0#sdh#4WSI3VVTqzzU$P=6llrBiUd?-m?ZMPS&V2-ErMutU5RvRsX@JS0Kk+;gaI
z=@1^;dm*4|H?+Atv7<8>bvqL<am#DW+{{IPt(gl4e_{MNVY<Lb(1t`w;$1L_I%h4Q
z=*R_BZ1D#f?7V0gbHcnimvHi|Ii^)A!u_5<ifux0|F}2?m&$PLdJbTM?ic7CTu3(h
z2S{s@D|zx9q&SVWw9{b`Ssj}|UsjYOJ|>bAsU{D~y@{rW-%zG8o^<X`BhOE=G*hY-
zC95*A^db*^78pmZDMQfo8^s<ODcbW@lOB1`qgv%DWa~VE{^$7Y%b3G#v=lZ@HmG1G
z{SLw%#Q@ygbpxwQA7T}|J_&6647HjpoUL$y`>aEFZ0`=0Uv5x+b`V|@wn6^WKFHc$
z$9~f^yh~4omvI)Xlp0ac(1pn@kD-xRgbN~KBoIA`CR)_s)8ruJ%n!wS-CR`p*28?s
zcNi}aCQXq(<mnY*;fq{sjH<yl`CgcO6d^@PMar}hqPF}$*!v+7pWB8w`?oJb*Tz%?
zPq_xiRUGU|x`E^c+aURH89q-R#F()TwRo;3i}Lk!;fFaDEC5-Lx1ep28_3Mwl=7x7
zBx#>1boE^+?(XX1xW1apskc`{)cz}&*^`MJpT{`*z7+nd+34^Pp`ILb3Ym3?CY3wW
zsvZZruWU)y`Ew{??{DN*voXQ{*x#ALO6GN(n6Y*I3Wj??pR8N0sJMPJeXm<XI|COJ
z1o+8kWhC~hdZ2Ut1st%t0EvLJa6EYmGKq)Lz2X}D{A2MTlfB1_3*azTjI}(?$ZGtI
z*{$WcbG{D2u7cDiIs*TvC3rt022+|MP(S?!;v{aPH{la1`M#i8?gnm*a|1QGqDFEP
zdh^YZDR33?KO^uk+8C4MZy>%m1=Z_a(02a>@*@hN_xm00%e=;tNe^Ji34~OIIPu>x
zAmIojno^`glh)5B)eo~sX6__H)Kofr+nl<)t!So=F<FVvrsA(3aVmBNwoFgw^z~|>
zXSpAiX(xjBRxr$?wa|AxoYQqy!bn%Sms47&h*L+0Ib5$karP&=aio=J<IKHk%suq%
zUwCF*8)a@9a2xZ8O4Bn#3$ijjL@#pA&}qLT<Z8T>Orv;cI3W+-&T&|EBOb=(VbJ^K
zfjRzesQem;PX7!@hDPG%%d6NbRe*-ywV06AjmVv!FumpxT+JF$r2P(p+|Q6dQh<J)
zn>ZH`iwi4n<9EqDB&mFWTFX_e6yd^M{_WsaK8+ulN;oB|hT#%_TsOaf@0%B6^-Mox
zde<Yyy$x9-4cPej3F=f{<89d=^mB@E^5Y}8%@!o);W$cEl_1X)B@+0qMho1PN$p(=
zlBo#qd<3boa~8F(Swa@A3#m?8n7Do1F@32Gd}CzjSf3f`=-E;3xnorH#)Eo}UnG&6
z?v(TK1ij%pK|(x7=}etHwVmEfZ&n_m+>_R{c39<K`+H0SmvKQXukqk25i<7Dq1Hn)
z$?1nOoqMT9Q<`VfyVyze;7lV9Mc&5s(|ItsdmAyvSD-5AiQB`)2tQnd*VZA($i9j_
zt1|HIYa^P@b)xIr3oK7;!m0&@xTc?rB)=B4Y_EXX)?0XgBMIU6voY<^M+k~Pfy&is
z$bL2sThq_Ozxpb?U)W=C{t2AbPeV0-KE`V%VRZQ|+*%d`*>BN!-VlMsbAw@P_X>|A
zTd^_k61+mw@O)+wqKkX*IByJDZ@)p>cn~io3!q?g1@4&z2)x~k(sPsP!}$4hmfhD!
z+qEdj^)4bCJF&Qa9UZpWMAc%4>7bhjDXLv2^HrBAslc0_HG2|It~;&hv!@Hz^XSrE
zAxe5JPQAWQ5aqDwUwAsa<1%)1<1<dZp-A=D49NeoD2X=QK*PlCcvRMpn#;;GPh%K5
zM=LS^LM=8eEXK#i5OkdjhvBE!s4K2Ty~hO%KJh}SQ5;gdYGAyh3D&vKaW1kALmi3G
z-WP*`;C7e>)#22_1W3!XaX?Es=1zMDp}}ZG=x%@$OJSI4Uxs;GF6N}wLGnW*ob{fe
z=y(-7=Tq^VFBpr$PQXfiEk1Z1z#+3t%*w9Dl`H<Jjl6+JVO#KU?|yh5e2Ve+KVipS
zZgOeir`p=<I1+giy`9f7$~6SZQ5BlHU7ftbrjzm;En2cC5sV!dRq^kpabFDSZK)1f
zyt1HDZ)@6AqDOwSWGO@GF)q)@gJAnJY%O^YgM=)cT-FNhL;3JlrGMesaEZ(Kl@p(_
zU9mE~nrTkoCn%GGdkHpcUB;{F+@wE2j0!J`l3^-88CdX=t?x@X%H2e9X#u=kx}kWw
z7_*oAqUE(GmIcHhig^Hj(<&Iw>%q;Dmsl~*7th6Fu<lYfOd6g+DLet2PbNV#y$U&s
zX%MQ)0?)U%koA3w3kq#Gva1=dha1o-Q;&r0jri@x?!!&#@bwJA%TzD)PIQC*v`F-;
zwP5<%i^$Um=KPRR$J&bdIC#+n%RiQ3(XT<A3Z74PlT@i=b`E++v(U1)6C0NcP}W^_
z%E_HgAJylP>^v1JayCZFNEh5Xc9FspYkF%bMILGc;FEfWr%T&llEzKDB(m@`JpqeS
zUZZ>ccq*=(N<x*BXt1A)tio3Q3s0|2T*iTu_>4ErRHo!j>&aPgHLYMKl84+dcD89z
z#W+c-OZ<cB6@v6|y)>zo_o0eE4Ra1>qrAKcXE^sEE*F7I?3~&*9D|;@_rNdy99g5k
z;CZhSGb2x+d1@?N>)x~1*bGca#yqpD7}1Z0_>C6y`G^r$s5o8f<{@#@KBVbChr!7l
zXr&~AOF0>fwkF`0cRaHGgks$7NQkE2!GYdl@U!b-Xvrj4)LcYUVHjrJ&B6Y<6dV=E
zhf{?#Wli5egMUQGaAprDkkg16YC<&TJ&r7Nr_qvcvuI(yI)%3hP|xxuxNcj;&W#P!
zB(R;d-s@2B0udV2sl~|s-+1^~oz`2<r+~)UbbZ`(dQ~@(*42yB^<6wP===q{9&h>=
zo(b!?jBOivjh!kd(gjBc;-2D44o+&6{V^2c(HfKxIf1%<$`Es0jUt?8k{ug|Wv;1a
zW10+z=O!an^e$or5)kVWj6&;Z$odsQYsx4_&3{AmWf<Hao<R7gXc#DWBGkPW@%+)4
z_oM*ThkwC&s}dblWOL6eN~E-3h=Pn>A>?cUcmfk(-WUj>Q>$U(xCMJ6uHvy{5-!B&
zK;l*l1R{B9=y@CLQ@f!fFHcnh;xt}<5T_gO!@0Q<vrKj9W`!Afm4xHwA5-*{d_gbq
zQuXP6hzCf}oGv4(O;@8j#a3Ld-w6w!6!>5T(PJwzZCOcGdu?cTts3n;w~%&f84>^Q
zMP%$Vn_4xe(w$p!lrF|aCO<x7K=C^y;vD~lC&IXlb4q!QA9P4iOVS=Xvm=1AH%_JX
z58^R2#G%(^OG)yj1&O~iC!^Pj<oUG;M^D~^L(mOO;R%QH<_!Gox&^DXfw;Li5+8CN
z!B}39zHfPp>5*=j_RR&^JHv3y@-cQ?Z-8(|2Mn48=ux*goi6911!6V$5YPxMnOA5`
zy^VVd{E-%wfbVKoanFH)ko!8EknzHNpB8+IQlg&$^NGfcD9U6G8>eWK#Fz#>Yn@2W
zqGGf@yAOSA{6A)BPRFx;Lw9#EN6+vgl=_>oNM`^|C*EMMhzuzjn$p)s1-htHj*K{G
z@U;eGx`6?uRxKqJn_X1#(VuSE>(kr0zwshsHdUHvQ1$ppRQE-Oy0Rqc!ka-9dY9tP
z{#zJ5<@_%^%eA?Uhh_PUSM1@TWoGMXtyB;lP*)_gib=HZ^>UhAvw|Xbt)qqf8YC3)
z5F70a&>Q{~*H`DjFQFXcf0w}4Bn}5;LlEs2g&Dg)pu+AuMsgCcQO*r&*IY5`9fm(_
z9>^oqiSK&f(c}0EozBVlF=08L9+N_ld;pw|-$S`jDz<;Rjr^DKc(;8ej;~pXrkF5n
z+4BLPmQAPXi>9>g?ONI`W=!!$Din53iDp;IlG<}gD&9Yls&7lu!2}Ju^+1}w>+0i@
z(Kz_~dO_6SE(RsK(4^c3Kebu3{OTg&U+@mmm+#?i-Z^}J?1&@ov*^UsMP%CPMYb9@
zNbaK<Ejr$cYj1@~DL|dLFHNOuy9e=rNr(SBH}Ltnf-5H!32RUN3(wRpZezO)K4Xi>
z&$yGino?)D(HYGg@XVP=F6XsLN@EGBPBEjZ$CA{TSjlo^@9}-#H(ax*XXDE%cx{M5
z*Xv6N6}*HWt+z4Y@C^OsH_^E$6wQBP5uq1|F?oNST@VUxuanqcW`^zcK;qbXlx1(m
z`oqg${w5ed1PfsMG6z3w?<3#m3KUlJV!iAZoD+HjC;6EqBX3Br-fEKQs4*QEF`*SJ
zRmst19Mx;T#;=G}6sO$7&nKzq{XT{->u%w~R8w@=y5b%i%RSwYjPs87aM!aJ`r4ZG
zo^LA2CpJKLs01z(ywLhc5CR|Nso<goJ@39s9qWQAFLW*auAe|3W3EEew+4LLGGx~*
zN}7di92@5isiaVJE`5%9UGD$Fvu;0+vBM&MV~y}O6s+G)eP1@wj?-I_8#S38jyIz6
z1q`LD&ZZkBJs6nv6<c(FL96{fT8cyQJox}bUu}YH+D^D0v%?0*Xy_LQKwZWZWvi~?
zk6t<)ojuUJMj6M<mm%!x4!n<7!p1aSScm9h+xZK~VRPKKyJAq_8H?_Bnb<G#2s>k5
zp;)E`Uz|Qd_QwP&U8hZNgEc92p9%#}Hl#1vGPE_n3hTu#W1-Y}ybQUDHR{=j<t@S4
z-D31cwheMm9wK$%0j4EAg8q?*C~~dA`yg40FjOZ^Sphn~?=^T2m15}D3HXozpeat>
z=g!maZPzFy;}~f_*-CT&=+FZlUsQy8A)ZT&4vrX5a*GxTKM^NG8o{VaA9k(t{TH6r
z**wN89r=w#1K&W2+msrgE+Rd52~tqfC+#K!;`yUM*K0rE`r>A^)(B9QT^GbP!!h%U
z7s8vi!c}q`3@2}fc%}_3R$j!PCNpHU4si6ccVe7G9==d8B<C!|qaCKOn0OJ349?@q
zN^7iLYl_^N?%3{?i-6d8Oqvk}<1^j36(dh`7*)D4e;T>VE7Hy@5_HYsEhfuHV4=r7
zq)gDHAKkJvSYM6f&c68J6pwq<gWBEG=-Ct_`Xi)7yOPsTYMzX3ZVwPGSA!t4Vs!EM
zz+${M33f~)sk`5BFZ3Pw9PYzQ*Bs}buLb|Wcv|uDJiCrurj)~5N!fWTHSiI6l|IG1
z35BrZQXr+6<#ceNCH;(^L$l*2(z8G@+Pf|2UwFFQ<1yA<%x^5S_$^KwO3{QEmQTrf
z0~?c3I5IyFX7m&a8f8cp<)+UKLKHVU6IvWLKPq-bl8-%vC+~x4@-|%da)5iA6$Gt>
zVZE!7^WqH$=`ImyxEc>Pw{RSlyo4|Iu9*LQJxn4A`&jORG%rB^S`ozj>sei-2_p7y
zvDowto(h#i<aaFISBD~!jm5(X!qLz2qoPLHcsQa&-@ohA=?M!+NPh~QFqNeXzeFg7
z-RJLR%2R650|f8B1_$j_eD<gVSK%|L%&bOr*%X=;vY2f5@luOUH&Tq>qM_;@CN(SI
zb7mOK=NZxGAFfn;^9&hZTu%~pmQ*HhN++zP$h<urOCQ%D>Ct53&s<46Elucr;X<+&
zpGgAk5&yz-#~*Iv&BSMXj$Id3DhDClRRWV4ELZ-_9ks@J7!zeVpS8IVn!@JAp#t>E
zHUcx_0<q&M%QNlYf!Eu2fal02jMg|JZ{<!jE;qwjd0lwKEJep9M+C84pWuWbboz!M
z=9)h|FWF$K(Mp_XIDv@WH{qC71@Q~dv2I}u44qe_!_padlzg#r&rawEIN`TZDYjJp
z!dNvg8I38>_Zfz?C2A!FZ&^wi*HtN{fR}oO%Mmhf7oNYF!O_ySgHzRQB$egEQS$}j
z)L!G>zD7J%m_r}-YSA$xHa0%ria!Z2a8IBRv&J`Zk~P{88FqnIXav#onhoTXx{$mi
zG%3tmpArk^QsarQ7!ZGmO)q<4B&$JhTv%RU`v}q<TG8r%>tA^83*s`~H=oz|!{jF@
zkAH|GZ9ee2Hq3c1I1WM{Y#hyX7ja58SZXRn8$$Xq*(Va`M6RRPzzdxC^%%LZ8oq}d
zF#FbStomt#P;C=bpIwb>OV7e1)CX_od*j$r4>ni3f{@-QFbi+7`&1}i=s4lh9!sdG
z@4=M^!6@Gqh_3B6=r3`Hs(w0(6{^vi@B~$zytKbXj;580P|T1p9gR_<%wzM(X6-y`
zh!dh$AsHC^bqHmDf-o{1g*%%(;g%Z)F`HO4ZG43A8E?_^n3q<G8&NYyl}tZ1Blpi^
zIM)|JUFsJ2uf5_7=W5XNlhKr$=1p_!)adyuR);wuOH+1hQpD;-6caa@5-bGC_GtjB
zKjolsaX$`oy~4qLPatWL^e;Sr+#O>E<#~-W+p8e?Iu{!EFC)K+2Q~96IG?vg;Xp$f
zo-__Yyq%l!HsnC?Oe9XR`KA~r6l-Rm!{+YYxaNNxy^_0-Teu2lTQ<Nf>j2I-U4~il
zHFRfRM$?}VH20i_U6&DREcYYzK>%DTS^adDJ@%RJL0x+&YzE#yY_BN!H%O6_KdaUC
zenz)iF*f=?LB~2ds_hXauc_h`y-||l!ez*Mjx6<Ljv{A48@{k}_w2(4+*I{K)TN_1
zzBU+qY%XEBKNIUM9)bV#EP5PnN~19XB)z)^qD|QtoOKf|EH6GcmLK~DcaiJ%II=yl
zn~bCsC?M_|5~@d0+M!B~VGBt0{!&spGMy^Nq9AoC0iEhC=-J3eCLB4ER!aXDp7+`O
zNJF33m=(++dnpyu8?T}=ND9-^uW+ugIZD>O7JQ5orlh5vpu|8#7+=SLeFBE&WFp-+
z4lC~ZqITIC*e=|GpBr7!+3pMPE_WoGT0&@sF0@?;(>Dddo#!GH(l<iic|S}KM?;=z
zfr4Kj;%s{%o%;jp%@k=>*F0KRZA7`3jVP>f7M-~=fZqboG0jkonk7f^bweq*Sp4)Q
ze<0#L9~}-BqI(miseOkSY5Ub7TVoGgv{~N1Bo<l}hN%;xaF6>gcBHBhy62Idfd~!Y
zJKkE<A^v_ID&=D!<>dvF)9z$$5lqf8%P7lLj8<ATz?XLzpXSXY0r}O$CE!ALi;Ss;
zUziHaC$b!?G|dm?r&udtN?norFFduRxs5v}vv_Xv$M5z8WE)4}M};~v=Y2D}IOh;H
z@DH=OkO)O2-N(=PXq+jE#`5?$%vMRk#PO*(u5btApM{`?QOCtJ2`EpR$#QDTU^l!1
zdQO@UybgSS;ffy;cj3R`9d?QiV29%<EK+%>^eGSZL~_$aM={cSqD6z>IK(Hmin6M0
z=*o~Kop75#b90ix5mz8AVE1fWWvmaV!HRW!G-x!A1~WwI-ED0W89$x8%PQfz!vo(G
zSK;}w)i^zm)dZtX;bZS}2nLJLSCh{OX#0$gygvNYdIXiOGL*Nl@f^FBEm-bGuXvp)
zsCXG0*UFRRSrM|Drb)JwZD@|xDca(4n56o2X!}$FdVXO6r8X?1yVK@VZ@dC2kKX+U
z&o_7HF=oyZ#%8JYj5D_n65ii}>H7qv9kj%vggT=hZ9Tjw?Pj&fG31|n0qdXlAQPSf
z(-qN>-g_Ajt=tis8-i<71JJR@0!^WF;j?HS%du@mLy-qAU%iL~$za@6yu;?Ab=b^$
zAKV}JqJQ)oKI#0x-lIQpN^Agc4W(&q%OX<mRHtBFahfy2O(nilX#OSzGTb^9o1gO1
z$dVJ-b2))CyDtN?e5EN+M}ekFE3>@ua$5P<d|p|KkJhfJdvA-F@*S}I>WuLLEccls
zOTt596gw#orw$iDIq@mlHk6~ws1iq`B2c-r4}t3r(o@;>WPj3_bf#-iz1LjI+^0=%
z>X*>W<HodNU<Mtwo=l6vS<d0H9;QtZqaPCuNXly*b6UCOA3UuRD;WzJabt%Ky3CsM
z*U@5>3cnS3*fib^pWnoA1T~fLS|bCSztrJQ&t3G*je^}$KfDzUN5pV2#*Yg`>@z={
zi;ZM)ZN?jyW;ADI;Gw@iK5aP-lJvl1*D%yYWMJw2`*<N&kKi9K5MtK`qrh&=$oK=B
z!|(Cp<pUg57(mGKGAJp0f$k?RS{0{3MX@vKc&jI*pT0zgoD%%~8aY~5laU@KO`edU
z3qi~1)YXkNS3s6{>&hT2;R)l?i!e^{LtvaIk_76Z&tmUu^9Fa;cwyA}I+VOBz`eK_
zHk*r~ShfU9j<><1d==TeG$*~dg|xoKfDV?dB=->x6{%W~dD%)jyUCm$*=!{z13Qv2
zQKR`EhB*6ro?xM1GGpxa_8&ZJm;Yk4Q~xmH_A%gDaTE8u)6r*DjwNg3(6K`t$5I}0
z464uIVto!W&zIprDXaZdgkY+1B!o(1QR0|^ip+f6G<t~Biz;vxS)k4sY`PVRopa*x
zMlA=Ww`=iOqZ{8%e<CS-47NV6vGCG$7z=HMrnM)VLq|g4OEferxyjq<7fvTnBGcpZ
zsD!a3D>i4g+sBRE*Q??7eJ(8AZ(_5&0A0*rdCR4fX}q2pts2!P2c-e*)@6A<w;0^t
z5`{X;OW;z!jD(af=vng8Mu%+F#(1GfJr2DOi{KEIf|J#E(C+_*W3ZpypTy^}YqlOu
z95<gBAEJu~HdD-W2RiHFNGD(HAVqIu;vAVz!cP|vueC8Pwt0tlm&37G;W<-aTK5m0
zL8pco!Nb3qsG3BCq}_m#VjQ0DD8z!AT!dUF?A@ixnQ_Aqa?UX<FHs5S+DFLuu7jX{
z4Q6<jV8P5JEcqD%8J71b53qxMl{Y+9V^K9&4-f9&_*TkKJ$#>`6_JM3)2@SIHHwmx
zaTsd926^u#Q1em6u3-xdZK}t~b28-7lLCpL0EqVJligt_ic6YC+r&>`cxND{$tHp|
zG{B>`3cJTTab}Ac1#?ZLpBsDN-QJ0*Swk2bWcjy!j}R$SjCR{DEEkm_ztka&>&V4s
zgE08;-2e~IEj*~Wjz<;oNc^41nS8Yx51hoQwm_P!S&n#_rZ%m8r$Hh{^T@(infRQh
zQ6aZJ)kBk1c^6a5CD0UME_x^92lJ#xMzH(wKX@*D{GIXY{KX8MN`c_Qc(f!&!z4Ku
zJ#kfd?CcIT>lYmD5(^x*xdE}K`Pd_J4?D8zpuM>ro{p_3zM6}?mmzp}{0tHmeIdK4
z90K7laba&fl0JEYx4<9UPhZ7C_mkKmy#i9t#-ls^AS%@25fXM0(&{>xc|DIauG|c@
z$8%wl%EmxV=SfiKDygfxlFN~aq-c2;oR&(g*_(<%?uAg%J%BBHpCD$tG?m{R!Lvgp
zSa10ZL91VbOZz!w%D+H4LV$j<`OCwK$uxIj0(dHdu*f_BtCB+SLn9n@-!5Wv(ph8(
zzBGDM@dQTJ{b=}AiBXRh)IC<B(Cr#DZv1TG;a4JmOZGlnEu?eWx^zB|<zWI>kWM2v
zU6H)UIapKATnzkQc*g!1U^YwqVl>qg5Ij8^g~B%=?3@R!#$GIZT?EZF`nY<-3H{>f
zcr8_h$paM-JX8hc-yJx%r32Hp-D6|CbQqS{qAO?ucyy28((fdsh-M>A>JrSJx}kky
zIIjLKz=yYSa9rR7CT%H-=@OnM2Sc}1A6J~Lv7XgkH8?sfk9&#^?f0Y=aH7nKGb!R|
z6595z#)?5BbbD{X;b(>@e6tRzUae@>{)IcrdDwU30hY}fM9IvLsA<ZAa!4sA)@oAx
z{57<En>fug&w#K>2=YdQ5Ym4Ti4P^<fA2D8q+J2umNFyrq-QwL@(Db<qOdI@7ExYY
zEDybq=GZKt9UE2YouUFgYG(N;GbIvyH;IgG49V71nRFjm;=z?V=7)RjKjZVj^KVS(
z&<|#9S_14F!mzC;7#;mN=r<UGoq8whUkSx7k6lpxngnsdCcH^}2CbiEShB7ddDS;i
z(0LMvZr{S{`LQh5!)i*Wk}>OKHNMpGQSl!^vRL*5J}i%NR5lH-1l;g@mpRPd*<wFs
zAVWR_hbz6Xd_@$jj=w<cOMVLCNRa$KDe9FFq{jLinC`2B!bgU9t7C?zuHNALoCAZZ
zY-qF=;ZAcN>W-d<l<Rp6FD!>|&>JY9zlX1ZUvQwvh`!g(r!N7l*3D|5u@UJ|W%cd_
zdK#!Wrwsv*Y{Ut1(Hp*D&Y8L8xOaO5J0jze+2M`+pRstiLyGc*)M=fbDt)$}Ln*uF
z(9{qyx;ZRBQxx=Qt(qlunvH_5?Ga;?{lC^vuNPmK4#gqHLo*A4k+C>`B?PKj&u}~X
z8$P(V!g^&hHc738{%uhhb)CnArD6E-r4~iK6>wK|#b>Vb;G9gy<WnPvKPyek--?su
zH(m-*|APm~?0!A+2Cebs2t9fe_BQU&Hg`c`Kq3N9J%ns!GS*23A-5p_I%Da0yN$i(
zPjTS3?0}2QGvrwxL!iPTs2D}yvTri>F1Zi;@>ZM_e1$%id;Z><h$Xi!!BxZumIDRg
zef<ht%btUE`p}#A3usV(1Z5Y8QLO(3N4J!seK-M!^|+w>D;IlAO7Z#WTsr?}058SW
zkgd;6HsikI;puo*|BS(%)2%2I5Ti%$XVK2P29)%hLxK&ubU}baL(1mlW44N3O3kDV
z4pmIo!TNv3r`eRx%#qXqrf|~(B<;z<s_xtPCe{o4_3fCjy#Z0N{TP0B7Jp{;8JPwu
z<J>eqtl9Gj&#Ugk{ZBD=GOyucJc=w;ehRG_gQR*pR&1_Ba7rz5H$MZ1)tB3LRG>z*
z5E=sE2o1Ulb2fg~d-DiWz9c|IF&u7MVaTC;v`@W{o`4dJynMsvn3+%?S_gwrN4QPd
zi^v(9Vb9}(W32xn*}M;!F@i4tpSV{$1f^>{6s_2cvKJ)?Z2F3AF`9JMa6E0;U5)mn
zA=FrWL*9i7+)np~+{;YNE(U$QtVvmRD{0=;HDokv09uZAh@C7$hT~+&zOox-)-f#q
z5(Q&%8CrgFHo2@8r4pH8Y!aSDH#?jtL&us_KI_ne)s>8lYr{WyE?eHq9E$$R#ISm#
zmUcP%7qz2v(-*9lFT(n%^)Pwfk7VV`U?#Y5Bvy*!C+k1dwy(khrMp;A*Mz;XAJEDA
z1|D#};-kcStoXs|ydn9xP<<O3#dnZ4rx>p;RO9(rB^pgqFc4Y_S*-?m9le2L>1nuh
zsszU$XQR!f6jGB)AaOYlx@&IZl+Ss1w{C*r*G+I5-iC5tM{E%~haIf$J>%#X2=;wN
z`q><udT|$CTy^+qkqG(aNic~MqtOFObSvTm-1w?t`R*I09sK}iUOX(+d|)5XP^i*D
zDlVN%Mk*W<vD!qgpR9<-j>SC97`9&PsIJwJc3V3_eE(;T)%KT|n>&%xF7VN;U5yZ2
zzLZj$d}zbaHZmE~C-=S8%tzD4fAIW%^ds~7B0J}!KVbg-E_mD?hTXJR*e#uat%FrK
z$h=1wyH5$&^Ww|)9*#P{7bfN0N0nGJy#2o5SnMCvEaxKb_z`3XcfznC7piAY!M}eN
z9u^v)dQUncUcbQ+zo!VPuEe@k4HybeLh$|q*s$8^mvt7@cwB|LdVV4|>;-nddjN+K
z)+^@Ygmc!;u=-<<IU<g*{%{48bZ)bL=;zq^DIanpei;3E6)|jIX=OG9D_E~V9}hhq
zo=Tq=@K7NCU2F=hg4Be5oXIRhjJG2sE0pMJ?*&rNT1|W!k|gFbix%tYk&ekydgQ4>
zf$}Ht#M6Xk^f*#}ge`<;=%dI#0#|p3VPQu<nntJ7$m^4&H^-I2oL5oNOf%Z=S;-ih
z{;%=5==n!xiR%|eie1Z%>;y>uE34<J=Atk$2^*H*hg923e9O9n<f~gzw|@~%UGYYd
zMI}u2KA^*C49xF7^smi8f07@{zk8$1md&3STb6$kMb<PG7#rTevKj4=6MTmso`YC^
z_%VKzUBI!E+N`d!m5c_i(A2P9B;q-Vjz)>lcgwf<e(xeqn4H9miDyxA<sf!kJc5+=
zD0p4zLimBNxECFR%r|ZjFYCmL)56rfh==TtupZCFA~bm-FXbv$;#25ttT~X4FKHQY
zjlG2<i~>pjJV5b^C+O8fA+kK53aJ&Z@s0H}nSN9we<dyYHS!$Fdz%p#u!Cm$O0ryk
zE{t5Gaf2%aCaS40d!a{{e;p$oR`V_VWlDQP9q8Bj62^kJ`k(Qswd4~cHSZg<V}$h&
zLyXq_>Of&`EY4`8L19)2Zu77;2L6RuU3(ihCSFHV@I}~lp2bh+W2}$X3n$eAAh<Ue
zJS7M5nf1z~or%X+6<17{w-)x-SK<EGD|jAX2?v2blw9N@fw*KevV8vFN+r7eVJG5>
zPmt@@qcnN97WK_gplRD$Va4&nP^$;r*ylW{voO#0#KMAFsK|_?khE^-o2H=QWj0(F
zmcdk{3h(d~qpSukv*shNE-l28lgE&Ca4|I8*5Nc8cb^{8q!f=Cl)Q^WBW0qbu*n8I
zF+P}FUyi*UeDv{`3JumW<XfXb8%^$TPVSgT%{30B^<gQ=-4dp;X4Y>vPM+G6_L5oO
z8p<8dO~1Lj@QvNWMn?0P-6Hq?!E?>F&&+(S@67mz+*FY-N>dNj0=$KIl3a#M8><oZ
zr2rkQzH!gD6DB)vA?AT1re-Q5(pdx!%RX_g$je}N#2FA9+huOeL2Pd(TdUyCdOJ7c
znc7Y?$os)$ZXu$t3et%HmT$}2h9}PkI730NS&k_ZJo8qNnSnF)2Ah-kfys10I0K8M
zFXPB&ABg2%g!2q<v@eLnP>=w*KV-F<9r=)0%ii0<L~MBEgEim1a6=^$k4w`qb6+Iv
z7szAE(LT<hkakXQ>jk`8Gz!TkL5j2!rMXkp$nn%OsICjaqMtR$zr#l<(`J%*svhMW
zSWjM?w$Vh-i4<IK!BGfgedXWx(CQ7F=`fcWtL<;3JD)a?Z_6Y~itNR>j4sGM{D7FP
z1<ZB*(tq%LKc$~pRy4%?{3TEQE5%4(w-A~K(;)jZAEP37Q9thncn=ienk^5RPDw?C
z)e<b;;ERhNJ@B+!8@JBqa#o}{bHw+kLT%jz6tX$9#ET;^v~xylARBLPzKXK`7Rb$(
zp(@RG+-aPHF##Wr{^(8^_8-N`5M{DRJ4E?u7Bs;6&~~o5iH%<ZQ8D>4u5e$%{?GtC
z?W;%h#|d<!>nE1@#$fY>T$XEzMQ_0+*ebA`{_4v(yUr8yl63I2NQcwpV1U=hSg(KM
zdq|dk!8zM;v}c(W`8celdd&&6U26;x>^fEV{s;U_gs5kq9I2f$ra57jq?=_(;jE{~
zb6A`+#eWvvPhUo|lUGyxWnU`OIYDw3BJ}q7ZRj3;#Ciw+VBg(hrv6v?KX@M8)Xy9-
z8fFq&ROoU(H+@}}4St<;+^egD@%UQss24zO@-QYXm82KSRmfi!fXIQ9IMT@A&zK;1
zGI<dDWIyse(%5(;7b-Vn!F=|`mzUl+)aQx-y<ohX!cVVniIc+J1hgya!Rf;jj%m3)
zs`JV)CZ$FT`yA+V^+Hl*V_mbkw{YNWFd`hUpyRqf8=qdq$3PJhXL&g@@p{DeltJ-J
zDOT<X!nNyO$WsY`c9<6)MeE_V%|}kJ;dxYXvpVt?De{xhr4uVyjcLvfdT@OcJ^jJv
zk43z6*7PkpblAMqyBS;hUtvm|C_T)XM&^#P)K)4_Tullrhn$8K`6{G2n9<tTr)d6s
zKPr;3rH3-66#Zg6y$s{04L2po;N@NB%eTsZ@YGrVl^NMP!d%(P)(N<B(~*tu;oQbW
z)9#DYc_|@EWao)M`xF{mDMIVcr@*ptE>fhFkulK&2l6kW{uVn<Rz62|8O#6gj6rPi
z9yEX0idmiZP|LH117A9RT$H8K!e5vo;fIqFJ5gD*8N1Dra8hy*a!M-n;jKQsJ0?cS
zbF&fS!t$?0A(%4zG8!vcy=ihXCg&^B-VF*gE3OzyelKC9lY^(WM-i^%f@;=#UUYj6
z78lQi_2@Miw07XEju1Q6WU1Lnl@2CZ(4@BQB(qeVT38PC<&sYPJWz$g1EqK?UWMBW
zKBF^fI*mRMr+xB2pw`X$_0O^OI0fyfZN0>asF+O~M*Yaf{0bFIoFuI^=gBMc2(905
zON_@VYI;z>ESdYi#^>HOgUs#TKTKh&JPmNmQA4>H%~BGh2BlxPu(1X2xJ2pLupDKx
zKI3yCt06Tvgwt`~6TV{=kog^nz{Xp!87zcQ;1#sV=^_7(79L))g^->-3>@6>hWUxR
zQ{{;>t{P|g1F?RGFKSl!K=*YBN(QrV%}j{Wi)Cpkt2xT;$;9WZB)G8r3R$>gMzt^0
z`)U!hS%SPi_d&mV0LpC>NU1@X4nK{<iIM#XSj)ltG+k6$-NZ@f5_oz(!f(T`Xq(E`
zq`jU%VpXQ}rp<u*CyUcHbzU0Z)`O@m$!wf|j`iq<qh!+*QopZDUTZp`#p>)&qJ`;t
zj|z>Ms?(tzyfl~9D=IXNsjB`EEq%Y2bi5r%W$`}RGQ)v>uG~nUG%}bIMn(UO&y?T7
z#uLvzX1pxM(bS*fBoM<zo>t$Wx2GIW9Sh*p!$)b!Phn{`5uC&W9LdT`&hJUL;p52q
z4j)9}yiGdh9*INDBzy2&+J^LxL72Gz4ur-PW4!kVn0m-lKp!7<_T7f5ax~((Z?bu0
z0@gO>!KSJmqF<%Snd=ii9b#*3SP$edtD!6r3B<0P-I(+43N(GHkemDn^19y;_k1?J
z*tnkbcW$9;PAms>BMPG19N-y}3az=1@z}T)ofjTs`i*YPH5f+BvdN_LnnMD4-*EVb
z8@A5%hmyyA_|2@slB#PcnVthN9#xuHzmRT@XV<f%<EhSElS+=8P{fXPB$i@K(R*$~
z%^(STHHfNT+R_{WGg35ILYw4_=+xXUQeRol7^wLFgXgXqVdJb`17`Us7p2t*QsNdt
z^4}>*<5-=kLQR4S5@qOp^aFg{x*8egdpJjSJ>gvB$;Z2%c1&k|@iW)g;YmRvsvktd
z`CuHrUQfilrUH~SHR9qPe(H#lBr5F0zVZ18zMqEKF?rx;8WH!j41#X$FiPU5i?`At
z*}&$V-%B9mT7+iV7`O<Xz`=lEY-DQ;ZWztMs}5^yd*g#S`@|_R+lu74^=Rp}GVqDr
zho#S31owBL)1VuL-^y`uTO>?=)ImO3n>J-FBqxg^#JusqQf)s}v);9o<%7^V`3qg$
zZ0*QdRhr^CpL$*|ro9ombc0Ke!dZXv(c7BjxK*EW$0m}5STKiox)@0x-cBx?ElJ}y
zo2&MykyM^DQ**zBX>R`WPd!5Ch9PF-m3sVA`38XzA^Mj83#K!>z(40VQgl>ED@&GY
zRbJx#{CK=>R)nNT9LK&c9d*-xA@bfZ=5M|az65VbDqq52_F;TKd<OGbZ{kQq1%~^7
z!C{b%nYB{TEyVgYR4?FtUJRN?%28aLhM1Quf9k>J&gXOCxtPsCz276Nvk|SFNQixM
zz+F`rY&p6VnZ>U;I;M*-XGR1Z{C*<LQk@<P%2K=QdmNwt32l<mIOpU8ofKCrHr#}>
zzB`~}+lR^R%jr|UAq6-+#=@Ql5VCJZq}fMoXKVlN+On~Ny$~&$Axkx%#uKlvI0Y&2
zl35rR&ATB?1FUc2mxCTvE?1#m(Ub6Qu7gYUE?U^Ljg(I<BAE^IXjCzbnKc~8{HJ8|
z?|PZ(O5d5cVt4TCcPTbHaFI>-7;0Fr_>bShBq7MIRS6?lV!%a1?p>?~bRDXxD)@6d
z4E~CuB-Hy8YD;dgxpg46np{NpGHaCeZo-+NW2_J7E)?Ik;EY-jKCyL6($n`MXU=hq
zYstj2jtb0J^#Y?-Tr@a$7!O&EE@6Es-e13t`RQ9RL1PMj|1rcyPYXP@P{YMSP4Mxt
zdj5)JgbX+0`<Wl`XSLV>L2q=<yM|d-`yt)B3z35CzYxf0YjU3zA<lRX?UQvR%glK+
z6rG0%I~VjDhGRlfA<qAL1vzN}5^fv8pN<z8pYsE^*&3Gj6`x>cUJaQ=jrj0+4lP|`
zOj9-r(#C-#jB%79GsdCI&Ij095mWkhYYTDZU1Mr36aR&$hw(Q?dSVi08fT)9=N-nX
zxar-4epnxsqv_&eba%>QB=$doOG^oSZ$-h)$sKOe$v7awOU`Xg(CqVu>l$B7UF40S
z56AJs?I>iU{1Bn>0!vs8uTE|sUQXG`#sz1<!+MPGx;|h%DlcJIGmO0o-8k5tiQ;Wn
zKwc)0ym$b{tAb$Y5(r)WAlQuagDSh{oZhwv9=qf5L9HC)#ypWT<1#9$!(g`YG8#A+
zk;(oZ?!kN<dX<Z#oLWTpDbfWmC(2i_qmPGK4+P7voIU1&KZE{|c1l1G_dV8U$4iBk
z+@v@@ldU0%MNHKr2rvJN^}pEpB-jJK)${4xE_Iss;U}6OeS%-nWhiZW2~C5;^sCf^
z=4W4|BYxq`TGNz&;aQKb%#XMjyx>j5^`*tIA9#vWT-+?*CC%o$U*Nf~3XxWsIBFY&
z-P%^D@>53aT9$7Neu(X*H*mVp9Y0Q<K-pPWh#Yf;{i)M<+Mb5x5zW|;n1LW&cVs<1
ziLgK)_=(2BO0N)FZfzJC{~R-T0-<p+5bg2FkX9^UYi{xpG%Fr*jqGplKZ5LQ>v3XW
zKR#3+N7s}nmY*xff{Sbo;C?B3wnLdz9}7`oNCO(yCxdHq4EA^yAxgF%%ls|LqRExk
z%dvU(y&n*{nvGj$0#W+-ICkwhjR1uV7*FiPxZE}@+ExntJ5i9St;B}levH(7K(oh0
zntFCR1=+qt)ZKCnJ?O&7t0A0^n}&ny*HeIaIvpGEB$a138E-q*YxCFR@A)3v|Bcz5
z`w+d4K4W8wEL~x}yE5PA(X$Q(x@PwV2ZF0X+6iE97Q9#TqpDhlQ(C(Lt51eQe$xd^
z3^<QG-~FtQX&<Uq9*6HT)~kKJ3Z-Z2pqG(?pOL{>{J{qyqW+k3G#)NrOEJ*F`a%*|
z&F5=A8}mE@>87Hn{1o_KZ2?cLA5vGvqIvWvyic!#!gO1xe+`CW=5sXXvOa?)Ms(dz
zpJsY0k)jnp72C3P$1CpQSqEDu_PrTTcB|14uLoVvJW6Go#3?T?0*9{X<MbR^w8xDH
zoc1Hs`!3{uzrenW<#<_Hi*qAS5UbORJGtBxmnu#n>{_wn6c^=oUBf$xn@Ds>huzc9
zSbjYj%hM*)9sWEj&N)U+TsN7*E36mnf8g1EaENKs7{S?3tWMfyOOHP6p*M@R(&uH9
zNNZLe*h*+5oeqVO{wfU4mP1NH5NF>xThs=*;*a!sq<eX=`H>3(A8beSO*@=1NX9Vg
zu)`n+6TfC6z&{q5m#!i%k=>8hhQr+P0X~FxK)37`^mLBF`TkW5EDy)qP<EYY*n`m?
zw&vq(Dx`-l<N3K;7~fZfwR8Ki^f&88d_R-ymE_6Akee(I4`b2NK6p%Q#Q)LWl?P+l
z{(F0titLF*N{fhABG2dh+)A55y`&U+mGnlXJt<MLl%++R7A=I3r6@_3vSrOq8_Jdz
zNrgJsncw@)d7F36oH=vO%=zcsGu+RZ`*~)3pX<Kv@AcWP!(XxxntcZmP6?RII7o?c
z9%Q*@B7Kp$fKQ@F5LBauY#C9cca4LC?rB({0<B}(FnPrr*uQSY9*F^{giDe8(~;zK
z6!dh{csg2;hERh?_^xypUY57eY>|$Q0@?628b&e{O%C&oN#T$4{4~1sOFi}FzVoAJ
zE78cXd33|vnyP*D$>-Ehc;-Yxknt_0jQim~&BM`Q@km*I63t0oa1yuz`x}usBXbTH
zkDkMWHm2)ZvJt!OHelXp#>WTe;wvu^Qf+rIj%n;a$s}NH&_(Rk-HnXiFxc5N;#MEa
z3;5xKp&Ah=a>#|zgM7qpOT%ZSJRCa3uxj-b?08&-LBm0OV_w)LGsMWHpc7+vKEn2x
zB#52OK)FX27RChP4+$eYU>OAEJ43KVz?5oQ)>89<ZX`%p!uF4=5VyMsnVCTdYB`GO
zo(Xu#yrL%6VzmB-B-t_j*-SSDx=9)o#Q62c(pF?2qeeH5MZslV0HlTlLbW9ZVw<y3
zHnkR^J;fL=9Y_Z9_O$qR7~jMv{FizLJQ(2jtyU&I1!Ib}T1NIo(p2{2ByI?v!8G$h
zJl#K*cBE>Nhtos~;3iRS!(?h-tV=aJHRyr9C}obR#eSOv@JfAgWlIneE=A)ctKSOm
zR2)>gh9$diqE7NIQs=S1mqy}{5aWRkH(|G|DD@o@BCGgL+zhG30l5meZ+Z=z13ie_
z*Mm=W-w|f}lj+@sss5cH9nXD-$sUh!clLGcZw-Xgm=jQmKZjn+K<JF}K(~Pnek9k!
z?A;3TtY1o!V?N-<-Z1QOj)OXvhgD3k)FH-t)sfj4ENsO)b#c->J)Fv3Nz($BL$Er3
zDowL=B424^QkH7QZR2w6ot%VzuVn0F?^XQSg(-cabmgWxuY~#Hw@S{X750Aos~cG_
z|F`>n?!F)V1LI_=Xv8wQ@!XlN$_%3qjH9g)e-7I<LX=!BLCc;f(5gBfIRv^;ulGI*
zDBD6^ldZ`6!6cgIDNaUPo?-E-2<9ITfyM!*$;h|>t%+I4bk0Q3<to?(_rt<LhI;js
zm_9~{_C6U&otqWti=!AFy8afzeeLM@CQ8k}tI^7{x}+&Eg$zB%(+1&T6n2Se|6lpS
zvD*~|Y4I>xn}+97Ctz{+96sKUMB$1@jOQxFPl-nus`Co-h8WYW(kV0|wF+D3hGO{J
zyI9fL4A-nSEYB`R;<I=>*cO8?F$KWLN)$XU#c%eKw7YT<UD&gWB0H66?fs{?ogq$V
z$`wiTrwr-#h>~BB2-SR&CXHS8*s>=Rk`<aHW_*ahu`}qGdU`nva_)?8i7(zlx2-mi
zd(ub9)aRpd;wLm+8^p9*pP?4QFqu(C^saR@&8-rmptDMpp`t@$_G{4j90^JbYDSG-
zCf4}{V~)uYNLyXS?!@PCzbr*grW&+kh#tL5n@Aa9BWTL7UYxtog)|AK88rQgJ8TU0
zsL0U4)tY3vV+@%p3?Y5TDlDFxgQ3m&7~vI*Jx5#-#-8!Lw-d49=mTuj%t3JSGi-IQ
zN1(`KxJ9L4%B^^qa0&SACQJD@CQ@M{%RPzA!a-g>dL7;(#<L27nlX^jS%}aZX1Jln
z{EBz8@K7TYy4FJ}pTC8IXRjvJk^#stZM{ZVGn|LLfpSA1df0lIN;PSZy&5eOKaI2V
zJs{gIPU2+;`KjUlztpqmz9=^wOq(OOjXI4N(8x<2xZeIAeFi^pzvVl^T|c1TMS)T~
zjL7pA({0;SVlYpRJe#yAcDV{ww*G+aVWv@k_Xt)-xp=En3M-9jmcQ@{CY4M(S0PNA
zHUiYRTac_1`=H+b1YvVyvH9{nXkUE+;R`vabG(6vp3L{@(u8owkI3xE#rDt?NEsHx
z%>E<93hJ?ZLIH##$}zIP4lVas|1#$l9-39+Ame?T5QE7@*I~bwX=!|3AfjH2zHevV
zf+r1#>%0tiQE%wW-9+NgD){C50z)<;kLi7M*5snCxEemjN;Eq00OcAkrOoHdF+-H;
zoBZOTvW<C$_Eh8CC?VQcGMi?v)1@*E#@EjJz$;A6LH~&Z{NRDZztl5floTg_lSlR8
zi|Aga7A5!yk<QMSXhb1&Z&qOVt9%4}lcDcvJ<z^BAM@&D$Wv(wxywwT*h8P;yfOts
z1sO1Sk&i`VtMUANJt(jo``U{!QN0|Wi)$d9_Z-G;9eC*T5(bqec&&O8({q<0>tX~_
zPM${`(>c7_cm$mv3-Qxki7M4}==wb^n)O<SuF4If!?pt}BN}nctpoG5M5()4lqO~O
zB8jb6_cs;bc1afYF&yaCo_RRC;u^9N)Tm*!8vU{5D;||UM&*)|7+UZUH6jw!(ftGR
z6YCLuF(2=2%P{?O14e96CB2WE$fA4-WmhL+Y42GKk6}7I{hNrIorTn&-7wG{OP+US
zsiChPrn|Cu_!5qUpaXoRmA=2!b8(wAx5{J?nvSn=BKR$&n_E%5Fb>*7ZXjTD5n^<n
z!f@~{y4o+Js*~y3E{-Qh+X*y3g89TOcffSZQJBw2#>a<`km>t`>BNeWTbP5WzzoJO
zK7?_9DoPtN(d&|p)j`SVK2ip~+m~>^X$oeFAIBw?8fZT4#XJiEDk~D9*NM#kD<(nf
zy~N2dv>(@ow6Xl0Z<sQ5IPt?r)4858<o=s7nKS+*<yisB;|^i-;3nn`WBB|vFZ^cn
z0DHEm(Pbqa+UBlBMWje&D&Mex`k|^bg5oYH&<<ZIQaw9>2i;%r`ovGvRd6I9=t!#u
zBuJ^3>De89F>urcPuv4w;gtfLvBN3gfhx6s|A=E|Z78f6gGSXj$R&95MIwJ4pI^&l
zx%Mtkcxs=<g~Zpe+x`lwH&bybC<~8OU!Yd>1sYTKBPm%HYPt_G&#C~|-p1pS$wO3D
zwd1qxd-$BL!Q0)eKc96UvBGze`0xGu5)N4lm%?O3B%1O{A)u28k<fTt-<6BY=d*Ai
zCl9_6WhnhtiDhk-Xc<zCQmVl=!wQ67DM71W8J5czU`0+9zLvbl%&k2bZv6}xn+L<Q
zCD7~tjuF0ZP%P??#O!sLKPwz^KN?{9Q-K<qr;)U^DfQ$Slg8r_bZqc8cKF4jMesS&
znLp);iZsO(O3)}Daq{1#Os5ST=zG;taz6JSBfByX7#{}FF)k2Y?2C~>d6;!eiE_H6
zD8RZ7k6tqV)y0slqZtOf;|RYZhIzUE_Ixj_kmJN|@5IjDSnN98hyC(x*sPt8wE7aH
zCAHy2cM-aWio#xDF1q2550`Jk)VvE@n0L-4P@7C*<mu|}YBap`#JkQ-2y8ox9Zyfd
zp>72xTXGnXU|1&$G!pfE41*^HC=-5$`e`k&WEw-8Ee!}gP>e~6kD*nR$+8(zaNH}2
z<y)jfL$d%W3(n!I?G~ix`ru9KK{#g}fTgQH&ON;iX?`2hJ_^#ow00~!P=TvQnsMNk
zFg=bMMU{n<=<3~B)DShDQuYYYJl+WuWd!2;;C1xcvrLSYJ?M$;!;#frkYjI3#_28;
zUpAE*uYE=no4cYbE~3)dAITxCXZ>A_mNQPK+pG)ALW^Ledmk<jCGb9MG-;gk=ZhyY
zFW2Add3BgPXRgkCXO?+T-2N7OWEuA1m<gFh&#-+Z(@W`;VDk$0xitxhu6qHUrjOX^
zo{zh@57D!aV8)&|4h%miNWO@xl_A)k#JmqP@-bDXfUU_H#&m&ck}qW7O6vn`HcZ3g
z*(?h~;3EtY`ta523+}x7g0BYOFj&-t5yIDSH}yC$|2VkRBv`IvJhgHFj&E{9hO-v}
z#)lzY;tG7%XJc4GJuao!K<q&Z-mbfiu`zA%G8;xJJ5*@YR5gk%(x+w0`gAq92af69
zsM$XUy1y?$(Q{wqUn)hc)gXqc3@7z&8*+ECq~Gl{DKJNzo+~pwg?}WZd(NS!zZNg%
z&!J`mEjm4g&6obc7~__TvPou$A1_6RNBQ$j4xRd?p6<sw`4-b3Bha-OOZ?km`Robv
z#uVYA(jZcn$dl~B``}++hB2alpxxsb6hvZhJ0ui4&aQ*|t&=Ew^aY146lsXB0tvIs
zKQD<-nDe3rVLdre3ax=cViSzsCSas}8lJe7z|Fe@kGK3l@SaX|UQ2@Fz7TBwbP7&u
z9%8}FN(eLETM?TZJY1`B=d~mi9}yt!FZmFk90w()ISiT+$TYjB!Bf2opEt>{VE_L~
z(=RY>9!3ewOFb+}j&!VrY0LqJI~&STtiu_Ub`Hmy0A+B)W}w4qA0pQzK|rey@tHiD
zF>yUD8>dHGXKGRI5k*ohYlKY69Z2$$5b{`!tUs~deI~<V5^E5ttVCkP%x|-CCw4qK
zfx&h^{_Nn;U+Nir#g4zVxD^lM`*G-94`xI^MdO?@oKP7-Ukpaj!_j*&Yw1{&=a?d<
za}5GYA~E0g6^4C&i&H<}K#Ad_N_!;8bNx^%X%nMiZN1odw;HV?X-MZPV4~W9##`B}
zcdkKNaV@T~@BUcj0UWM9k9%v`xQfoilH_V6#*5Rfag%82r*ZUR)=0|rk)!|@VTyYx
zOp)sw*w|z1Q>H<mSj^V?C(dK~!3cQ1&qZy72&IY(Q-oqRQl`a1W<8rHUtdMA$X66`
zq9oaC4*Lz?d8UQauzI^Awt561Pw^5|ODeF)WFe(~+(3qc2K3QYpO)X7K-Z)NXh4$b
zn~jRFn0bJ`1!ZV(vjct`-iu+SvQ)m+km?~oI{TC9tKvDn$|0v;>iLCV00-AfyonK_
z`l$mfFRv70jE{J=N`-E0kfMG)TLd$2xLrd6ma%?Mp{@<a;qMVY?Hx3Cv|>gh^K0*9
zT00bDb4DE&C%3|R+H0&^{RaK-KcT@h1#y!OV?(kxE{KOPe#8&rvsb}g{S@xGet}r}
zWa_gsr<tRyX`Y`YxfsqOk7Yy2>01p}?5e}mZC~(tY9F+eD`3j9x*k+s!S#+XNC$65
z)ENi#I2=OsD^J|=-i7lLXR-WH6gD}Q;LS{B(sdW67W+fUN?r&X{sLrAaDdw?mJO60
z455TKSZy+!3c{C?xvUwb3YyaT2onm9kf(sFweTux#c>lEsz31x-8OkxvRaU|7?#kj
zWJ75V4y1F!m7gGSk-zdx;lJ`N|HSJwEw5mAW)vp1e#iLY0kmzbhTYTmaC23mbnii=
zEj^Cij=SKc?*lvQQ1n%X<A-q=Je<zMO(F$GrV)_46bzRUq4-VxG?YfNxqRtO=-e*E
z0Miw&aZE?WcVA?(9%3Nl9DJW=VsKm<te=PD(DO)0GL2&4KpLh#f5<$DKjA_@V3tq=
zIiDo>b5Y=jKZM(m&q!GD5^aUIVB#7I{ew<e>FSKta{d^Q4TjA4>kK=31!+++`fZXp
z73V3^wZm$(YJDNX%(tTBl`}FfxIsAL094laVUyq~Y;S8rV&rJjZ8oRWfaPReVoZwD
z%xI0nI69r!1@Dhdu$bA3nR`{~)J-D_Z8xR^nQFACMTeSpTGPVcBKfPWPXB9c{_S~m
z?^!XDZCeoi@*N_XC%fQr4cZGu$UIe=hB}vHAm|Ft4h=@0-(p1U%HcVmnac}3c7-?G
z-4$n(yz#ltmhq5w@YrXA>?a!-ue2ZCPvRlH?lBy!87CC=5plBJ2<b0}dk*s;9ScE@
z_XgziGI;7-C-2_Lb?CTwm0>u|=#*@SPDKJz-bTZ$miZHN9PzE*3(-ofZwvT}>L0hU
zli_*K8z14llMqEz$x&#5EOp)=N(ZknE=Y3-9e6K7YnI57mH!Zm|8Nxs)|+tlxD#fY
zx<SXz2g||(5tE+=Ip=5i_OKh1HRe(8DMwnm&V)J|=g>%pN#vR*MQ@_p*!?I*qON1e
zra+DCq(o`fl6;Km{fV%++4Q0CCcoe1;y>!Syj_|?M*W2Kk4o?_7NeHwuAgc4V7_@T
zq&62IGdm7tMNFfiy%v#rLKq@-j<@2(B2>H##a(50?6+8oo<2LIr5Pf<L<zEI=i|)e
zhdBLNgjVQ{An`OQYIJSFQPx8^3AtiLt}RSOkE8wQF19W{js;H|k-vEWvu2h-aqe|&
zda@U$OS}-WeIJazTSBVX1@~9q#SX6?ywDXRxf)e^n#nM)7b>LFCqs)oq)4Ytlmt1p
z9_Sdr@+xUsRG~)s@~zk#8-!gyY!PBK2^q7;W8Mc_NU%&p`!iLD_>v3Z<K-x@nnoi%
z?MXr3jGTM)X*JVOs-zF67#}fGIIKh>%JQV`DnTtPo6t2S3Y`oa`!rmE`YjXqS&4sg
zM}HfiRSE-05eh|9!ZKXk69#LQ0?g0Khv1qA*x^tLYngb|8v4P))EHwI4Ml66IwW@n
z<9uf*E=~x<8Q1gB2x0v03}58Tm%`_{VZ8V{KNxj4;A3eA`cL)1XMZDfm{+HDiZ52u
zEohqMVa4{t*i^g%#)kXxq}vT^tCvGQY74Z(HX!WO6xbdg4W5}5{9V_=g6|6nv2?_-
z`S{x&Q7Uj0pg?{%E^K;<l6MVw&itJECe4UD_7(9#k|ciR1^jzFAo^e|VsGmqD`qN=
z>>JCnFqgyfXettIYS6?b!XfD?<f=x}9F{Xvt`5qI)S&9g<0xXgAf@IErHljPY4lTL
zYPdF!R;iDr*avlNZJml~OkY2(I*I>!+U<Xg&%^Bd)TsS}+x~Z<|MCu0G%De^AO=g-
zB609$3DX+KLO6wK5w$gtoUeybJ~z=jr;@F^o`JjC$g;iIbMjOcM3%&0zN#N)+nhtg
zoK<Y?>4S=ymtZE>i;3?CQS&JV`hpSA(Jn^6$2~};T4MwA7Q5MmK%C8`*-!4l%&ifs
z@?AJ@R1am=$7*Ld;Fz{77MDEZxf%$eD_stPu@iAZ+5$%kTrlwVATC(2Oq54!k#N`*
z#xa@DFsjA;cwZ=QoB{J`6L6{02;44fbZFb5!S)dAk@Ml(Rs!28ap;mR!0y|kRQX{#
zRfx@^QR-?m=DjdIxBkd-Ya~eK9>XlnS-#R}0rKe|O;`F=Y4XYYSYL4h-m!Q1-yE0^
z_OIj5=aJL*0z3|VjVR-4=70Onu>W`Xa_kI>RIWl<I~SJe5omu(nB^=7Zt-43$Cu;K
z_*U?1yYMLP4Gc<((c^LlRaJ>tZCVL+5jh(4WZrxGWZeFK7G}Ff(ENKM#QS<4E&1`t
zVwguO>v_}}_AOX?6U(M&z|$}r=Jy*|{@NEPy%M0PHQ%8t-V5iC84$@^i>BTRUg21E
z-X?<r(^fra(=EkUO!Eci@$zcZc+%s_dAx-Zu$ALboAeezc1*KdnTOp)w;&>a9zjih
znBnRH_ZM4O{>3dUpWB2dwT~dhV_Z&L0n(H{!?Z$$CS938O7auwhwm-Ci1x>rm<L#1
z^9H5?LUgi;_0j5@WS35Kd=bMOp+H?1GWjozGyhS~KFu2-{R);{#<DdlM$qtSotV^e
z5!=??f?O`ctL7cSg;xVS?PK?OYWxV66WPLY<(eSa{S1F-GS5V73@|zx$JU>LLTeP#
z&JAGw76FoQV>tViVr06iQm3C7O>VsdFN-{2UjoE*{ctub3OZHE$l;xaL{b1&viyf~
z`4$`%t;ea5NAOx$jlpu3na~)6d&>HluWZ44YgJ?Vqb}6cVecK&Ig>Sc{T{b?VMnWZ
zseVRyT9gaUV?14EUXE9yXR%K43f^YNFilc4<UANBGc*i7&k7+ql<D0z=b>goI<8ii
zq0No;VtJCZc*;bwVl};c(E}PgGf=dm8zOUNk;<7BBt2+K<MqtRNn#GY_Zdb3$;=}b
zp3jeKX1e6Rjz7oewL|jUWe+9J=&cI(s%A8|U{I5Dn>C*M<L(5`)p9a-Az=#VI75$n
z_DGMr{7sijjnd(siB95XCW>(3mpAekTsB8etu@n{o3We_rctyngS*2=bhaqc=M9EL
zq9Do1toOn@B-!RbW=8;&#`s~MMiQodO2T}Rb6B>I>6-5+VU|!fCbpl!g~lC_Wc{UQ
z*GYH`k42!~Vdmo!!!f@Lylr;wyhYzac?!tj2?q;eVTU%RaV9W#?Lz97p)8BV3jK43
zqD;mSHYR7_&Ae&0+j23trWSR-Gwi0a6&D_{d&DyqlTMakf%G8GhDg#oe-qleT#J^r
z2~ug5E^QMuCN*)6l%!|S)8A&(oBGZ4>75A))QVE+6*=mZt>C|0!+Z#T9e>tSDOjG1
z^i<*&6^!B<V#jdXO|>}jj`3V$_(bmhI2~@bwk{VpU61SU*5$H#bhzfhiJYvfHg{#;
zZ0?QNEG|c6ELZ(pgIn#O!FjFG<jU2xxc-gfxN+ylaRG<4xK=q$POfP*x4V55SA20K
z_x_U-_q&B6_gPq;TfcY&XMgx5UrN!GpR~Q0ai4eLynG}+EE9(8`T^dy&(0{Qi3Q0$
zg2b_NaC&45O{<wWGJHDAvN6GoBWn>*$g-M0hr!$A2g?92gT|T)-U4ea-m1_}Uiepz
z<!E@rZ(%g{CcK8u7<SKm=)id9Ei#vXj1OK-_}C{!#bdRp)WCvT9cGe0rjtRpF-0vU
zviRMQie@O$2PaMXptFY*<*aDfFl`E6%(OF34g4L!wg0H+U*-Mpo<o&g_`<&na;-vR
zxx<YbT$kk-&h_DFu2ya|*WIVa-AqyA>h#sPnyIRsLV^mHEX!(nS&7TqsmN^@p}=)N
zm*u1u$Z*r|sc|R2t8j0YNpMgS<G3XvT!66%*S}VT)72N{T+76`k+$Mo;86+gnWZH6
zJw%dwvs03*36$h^ER*62UJm8_zYgb;v}L)cGvD(K=G*g2H1G17)-C0QFMY`~wX}g1
z(^5rE>A;7KLU`qMAudFO=FF8O-OQoXbV!cAu48&IJ03}WVs(|7Mh{=;(ai>RD(>z?
z|M`<JuxvpK-<%AzhERTVD)ZbhuJPV$e)Fcr|FwGlckBP3UH|@$86R%-e1&o1T&0&N
zx8woKCA=rZnYsvZ{mDX{R;V!7peo8e8xZ3<rigQLA>y3N4RNlZMVyPiBf;4iN^whT
zhH>`xBe=+c@mx#a1g@`Ci7VQv%w<kh<{C>BxmXi<E-g=nGi2A5o{{1z^d&f#QZa5$
ztO$3?TZkjYAN*N|Kl8V2?&dq`b@NrnXY!9OoXX!Qponl|Bi`2y>2P6LHz!mjm@l>-
z#@89<6Vk?i^YZmS>-pz>|Ce6m|Ihz}{ZU2;{BP`bRJ-xN-t8!#UGrb<cGSkM`D?c$
z0oVW7?Pw}{{nu_sZS3{!KX*Htx=BF5`<jq|7#sT`?DfAN2i5=M=j9Lqfm3WpBn3@&
kW_$T4I0^{-wd2uW=YPHHKmYyzJ0Jbu?0huEkR1a529GH1cmMzZ

literal 0
HcmV?d00001

diff --git a/tests/test_zonal.py b/tests/test_zonal.py
new file mode 100644
index 0000000..b861edb
--- /dev/null
+++ b/tests/test_zonal.py
@@ -0,0 +1,30 @@
+# test zonal stats
+import os
+from rasterstats import raster_stats
+DATA = os.path.join(os.path.dirname(os.path.abspath(__file__)), "data")
+
+def test_zonal():
+    polygons = os.path.join(DATA, 'polygons.shp')
+    raster = os.path.join(DATA, 'slope.tif')
+    stats = raster_stats(polygons, raster)
+    assert sorted(stats[0].keys()) == sorted(
+        ['std', 'count', 'min', 'max', 'sum', 'fid', 'mean'])
+    assert len(stats) == 2
+    assert stats[0]['count'] == 75
+    assert stats[1]['count'] == 50
+
+def test_zonal_global_extent():
+    polygons = os.path.join(DATA, 'polygons.shp')
+    raster = os.path.join(DATA, 'slope.tif')
+    stats = raster_stats(polygons, raster, global_src_extent=True)
+    assert len(stats) == 2
+    assert stats[0]['count'] == 75
+    assert stats[1]['count'] == 50
+
+def test_zonal_nodata():
+    polygons = os.path.join(DATA, 'polygons.shp')
+    raster = os.path.join(DATA, 'slope.tif')
+    stats = raster_stats(polygons, raster, global_src_extent=True, nodata_value=0)
+    assert len(stats) == 2
+    assert stats[0]['count'] == 75
+    assert stats[1]['count'] == 50