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_mini_raster():
    polygons = os.path.join(DATA, 'polygons.shp')
    stats = zonal_stats(polygons, raster, raster_out=True)
    stats2 = zonal_stats(polygons, stats[0]['mini_raster'],
                         raster_out=True, transform=stats[0]['mini_raster_GT'])
    assert (stats[0]['mini_raster'] == stats2[0]['mini_raster']).sum() == \
        stats[0]['count']

def test_mini_raster():
    polygons = os.path.join(DATA, 'polygons.shp')
    stats = zonal_stats(polygons, raster, raster_out=True)
    stats2 = zonal_stats(polygons, stats[0]['mini_raster_array'],
                         raster_out=True, affine=stats[0]['mini_raster_affine'])
    assert (stats[0]['mini_raster_array'] == stats2[0]['mini_raster_array']).sum() == \
        stats[0]['count']

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 segQuality(inVector, inImage):
    # open the vector
    lyr = fiona.open(inVector)
    features = [x for x in lyr]
    values = np.zeros([len(features), 5], dtype=float)
    # loop over features
    for  i in range(len(features)):
        geometry1 = shape(features[i]['geometry'])
        
        restFeatures = features[:i] + features[(i+ 1):]
        
        value = zonal_stats(geometry1, inImage, stats=['count'], add_stats={'mymean':mymean, "myvarianza":varianza } )
            
        df =  pd.DataFrame.from_dict(value, orient='columns', dtype=None)
        
            
        for j in range(len(restFeatures)):        
            geometry2 = shape(features[j]['geometry'])   
            if geometry2.intersects(geometry1) == True:
                #print("They touch")
                value = zonal_stats(geometry2, inImage, stats=['count'], add_stats={'mymean':mymean, "myvarianza":varianza } )
                df = df.append(pd.DataFrame.from_dict(value, orient='columns', dtype=None))
            
        values[i,0] = df.iloc[0,0] # count
        values[i,1] = df.iloc[0,1] # mean
        values[i,2] = df.iloc[0,2] # myvarianza
        values[i,3] = np.var(df.iloc[:,1]) # varianza between
        values[i,4] = len(df.iloc[1:]) # neighbours
    # get overal values   
    intraVarWeighted = np.nansum( values[:,0]*values[:,2] ) / np.nansum(values[:,0])
    interVarWeighted = np.nansum( values[:,4]*values[:,3] ) / np.nansum(values[:,4])
    normVariance = (intraVarWeighted - interVarWeighted) / (intraVarWeighted + interVarWeighted)
    numberSegments =  len(values[:,4])
    
    return( intraVarWeighted, interVarWeighted, normVariance, numberSegments )

def test_nopoints():
    with pytest.raises(TypeError):
        shapely_to_ogr_type('Point')
    with pytest.raises(TypeError):
        shapely_to_ogr_type('MultiPoint')

        zonal_stats(geoms, raster, global_src_extent=True)

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_mini_raster():
    from geopandas import GeoDataFrame
    polygons = os.path.join(DATA, 'polygons.shp')
    df = GeoDataFrame.from_file(polygons)
    stats = zonal_stats(df.geometry, raster, raster_out=True)
    stats2 = zonal_stats(df.geometry, stats[0]['mini_raster'],
                         raster_out=True, transform=stats[0]['mini_raster_GT'])
    assert (stats[0]['mini_raster'] == stats2[0]['mini_raster']).sum() == \
        stats[0]['count']

def test_transform():
    with rasterio.open(raster) as src:
        arr = src.read(1)
        affine = src.transform
    polygons = os.path.join(DATA, 'polygons.shp')

    stats = zonal_stats(polygons, arr, affine=affine)
    stats2 = zonal_stats(polygons, arr, transform=affine.to_gdal())
    assert stats == stats2
    pytest.deprecated_call(zonal_stats, polygons, raster, transform=affine.to_gdal())

def test_range():
    polygons = os.path.join(DATA, 'polygons.shp')
    stats = zonal_stats(polygons, raster, stats="range min max")
    for stat in stats:
        assert stat['range'] == stat['max'] - stat['min']
    ranges = [x['range'] for x in stats]
    # without min/max specified
    stats = zonal_stats(polygons, raster, stats="range")
    assert 'min' not in stats[0]
    assert ranges == [x['range'] for x in stats]

def test_direct_features_collections():
    polygons = os.path.join(DATA, 'polygons.shp')
    features = read_features(polygons)
    collection = read_featurecollection(polygons)

    stats_direct = zonal_stats(polygons, raster)
    stats_features = zonal_stats(features, raster)
    stats_collection = zonal_stats(collection, raster)

    assert stats_direct == stats_features == stats_collection

def test_featurecollection():
    from geopandas import GeoDataFrame
    polygons = os.path.join(DATA, 'polygons.shp')
    df = GeoDataFrame.from_file(polygons)
    assert df.__geo_interface__['type'] == 'FeatureCollection'

    stats = zonal_stats(polygons, raster)

    # geointerface featurecollection
    stats2 = zonal_stats(df, raster)
    assert stats == stats2

def test_ndarray_affine():
    polygons = os.path.join(DATA, 'polygons.shp')
    arr, gt = _get_raster_array_gt(raster)
    stats1 = zonal_stats(polygons, arr, transform=gt)

    from affine import Affine
    atrans = Affine.from_gdal(*gt)
    stats2 = zonal_stats(polygons, arr, transform=atrans)
    assert stats1[0]['count'] == stats2[0]['count']

    stats3 = zonal_stats(polygons, arr, affine=gt)
    assert stats1[0]['count'] == stats3[0]['count']

def test_geodataframe_zonal():
    polygons = os.path.join(DATA, 'polygons.shp')

    try:
        import geopandas as gpd
        df = gpd.read_file(polygons)
        if not hasattr(df, '__geo_interface__'):
            pytest.skip("This version of geopandas doesn't support df.__geo_interface__")
    except ImportError:
        pytest.skip("Can't import geopands")

    expected = zonal_stats(polygons, raster)
    assert zonal_stats(df, raster) == expected

def test_ndarray():
    arr, gt = _get_raster_array_gt(raster)
    polygons = os.path.join(DATA, 'polygons.shp')
    stats = zonal_stats(polygons, arr, transform=gt)
    assert stats == zonal_stats(polygons, raster)
    assert stats[0]['count'] == 75
    assert stats[1]['count'] == 50

    points = os.path.join(DATA, 'points.shp')
    stats = zonal_stats(points, arr, transform=gt)
    assert stats == zonal_stats(points, raster)
    assert sum([x['count'] for x in stats]) == 3
    assert round(stats[0]['mean'], 3) == 11.386
    assert round(stats[1]['mean'], 3) == 35.547

    def rasterstats_statistics(self,admin_vect_file):

        from rasterstats import zonal_stats
        stats = zonal_stats(admin_vect_file, self.population_raster)
        print stats

        return "Calcoli statistici effettuati per categorie cicloni e popolazione....\n"

def test_categorical():
    polygons = os.path.join(DATA, 'polygons.shp')
    categorical_raster = os.path.join(DATA, 'slope_classes.tif')
    stats = zonal_stats(polygons, categorical_raster, categorical=True)
    assert len(stats) == 2
    assert stats[0][1.0] == 75
    assert 5.0 in stats[1]

def coverage(raster, regions, key_property=None):
    if key_property is None:
        key_property = 'id'

    with rasterio.open(raster) as r:
        res = r.res
        cell_area = res[0] * res[1]

    def _not_nodata(x):
        return (~x.mask).sum()

    with fiona.open(regions) as features:
        key_mapping = {}

        areas = {}
        keys = []

        stats = zonal_stats(
            features,
            raster,
            stats='',
            add_stats={'count':_not_nodata},
            geojson_out=True)

    def _coverage(item):
        count = item['properties']['count']
        total_area = shapely.geometry.shape(item['geometry']).area
        return count * cell_area / total_area

    result = {
        item['properties'][key_property]: _coverage(item)
        for item in stats}

    return result

def main():
    for t in tifs:
        stats = rasterstats.zonal_stats(inshp,t,stats=['sum'],copy_properties=True)
        df = pd.DataFrame(stats)
        df.set_index("BasinID",inplace=True)
        df.sort_index(inplace=True)
        df.to_csv(t[:-4]+'.csv')

def calculate_pop_value(node, raster_array, affine):
    stats = zonal_stats(node.polygon, raster_array, affine=affine, stats="sum", nodata=-1)
    total = stats[0]['sum']
    if total:
        return total
    else:
        return 0