From 70f0871fdc981a30960c9a03bf64a933c17a0a2b Mon Sep 17 00:00:00 2001 From: Tom Kenda <75068847+TomKenda@users.noreply.github.com> Date: Fri, 15 Sep 2023 13:37:57 +0200 Subject: [PATCH 1/3] add gdf_zonal_stat to main.py here is a function to get the result of the zonal stat directly as a Geodataframe to keep all the properties and the geometry together with the initial vectors CRS. I you think there's a more efficient way to do this I would be happy to ear it. --- src/rasterstats/main.py | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) diff --git a/src/rasterstats/main.py b/src/rasterstats/main.py index ad77377..455b747 100644 --- a/src/rasterstats/main.py +++ b/src/rasterstats/main.py @@ -35,6 +35,24 @@ def zonal_stats(*args, **kwargs): return a list rather than a generator.""" return list(gen_zonal_stats(*args, **kwargs)) +def gdf_zonal_stats(vectors, *args, **kwargs): + """The primary zonal statistics entry point. + All arguments are passed directly to ``gen_zonal_stats``. + See its docstring for details. + + The difference is that ``gdf_zonal_stats`` will + return a GeoDataFrame rather than a generator. + Besides this, we make sure the new gdf as the right + metadata by conserving the CRS + """ + + gdf = gpd.GeoDataFrame.from_features( + list(gen_zonal_stats(vectors, geojson_out=True, *args, **kwargs)) + ) + gdf.crs = vectors.crs + + return gdf + def gen_zonal_stats( vectors, From aeb0227d7a6101bbb0c09b0e814396a41f283981 Mon Sep 17 00:00:00 2001 From: Tom Kenda <75068847+TomKenda@users.noreply.github.com> Date: Thu, 5 Oct 2023 15:39:53 +0200 Subject: [PATCH 2/3] removing `list` from gdf_zonal_stat in main.py function `list` was not necessary --- src/rasterstats/main.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/rasterstats/main.py b/src/rasterstats/main.py index 455b747..72adf08 100644 --- a/src/rasterstats/main.py +++ b/src/rasterstats/main.py @@ -47,7 +47,7 @@ def gdf_zonal_stats(vectors, *args, **kwargs): """ gdf = gpd.GeoDataFrame.from_features( - list(gen_zonal_stats(vectors, geojson_out=True, *args, **kwargs)) + gen_zonal_stats(vectors, geojson_out=True, *args, **kwargs) ) gdf.crs = vectors.crs From 19e9f6027734981d7b641c96edd72b01de5f9d80 Mon Sep 17 00:00:00 2001 From: Tom Kenda Date: Mon, 20 Nov 2023 19:27:56 +0100 Subject: [PATCH 3/3] adapt to geolike object or file path + unit test --- src/rasterstats/main.py | 15 ++++++++++++- tests/test_zonal.py | 49 ++++++++++++++++++++++++++++++++++++++++- 2 files changed, 62 insertions(+), 2 deletions(-) diff --git a/src/rasterstats/main.py b/src/rasterstats/main.py index 72adf08..c524964 100644 --- a/src/rasterstats/main.py +++ b/src/rasterstats/main.py @@ -2,6 +2,7 @@ import warnings import numpy as np +import geopandas as gpd from affine import Affine from shapely.geometry import shape @@ -49,7 +50,19 @@ def gdf_zonal_stats(vectors, *args, **kwargs): gdf = gpd.GeoDataFrame.from_features( gen_zonal_stats(vectors, geojson_out=True, *args, **kwargs) ) - gdf.crs = vectors.crs + + # if the input is a file path : open has gdf and get crs + if isinstance(vectors, str): + gdf.crs = gpd.read_file(vectors).crs + + # if the vectors has a 'crs' attribute use it to keep the crs + elif hasattr(vectors, "crs"): + gdf.crs = vectors.crs + + # otherwise, we cannot store the crs ? (not sure if that is possible) + else: + gdf.crs = None + print("Warning : the crs is not stored in the output GeoDataFrame") return gdf diff --git a/tests/test_zonal.py b/tests/test_zonal.py index 02ef6c5..4f33b4b 100644 --- a/tests/test_zonal.py +++ b/tests/test_zonal.py @@ -4,13 +4,14 @@ import sys import numpy as np +import geopandas as gpd import pytest import rasterio import simplejson from affine import Affine from shapely.geometry import Polygon -from rasterstats import raster_stats, zonal_stats +from rasterstats import raster_stats, zonal_stats, gen_zonal_stats #, gdf_zonal_stats from rasterstats.io import read_featurecollection, read_features from rasterstats.utils import VALID_STATS @@ -564,6 +565,52 @@ def test_geodataframe_zonal(): assert zonal_stats(df, raster) == expected +#%% # add a unit test for the function gdf_zonal_stats + +def gdf_zonal_stats(vectors, *args, **kwargs): + """The primary zonal statistics entry point. + All arguments are passed directly to ``gen_zonal_stats``. + See its docstring for details. + + The difference is that ``gdf_zonal_stats`` will + return a GeoDataFrame rather than a generator. + Besides this, we make sure the new gdf as the right + metadata by conserving the CRS + """ + + gdf = gpd.GeoDataFrame.from_features( + gen_zonal_stats(vectors, geojson_out=True, *args, **kwargs) + ) + + # if the input is a file path : open has gdf and get crs + if isinstance(vectors, str): + gdf.crs = gpd.read_file(vectors).crs + + # if the vectors has a 'crs' attribute use it to keep the crs + elif hasattr(vectors, "crs"): + gdf.crs = vectors.crs + + # otherwise, we cannot store the crs ? (not sure if that is possible) + else: + gdf.crs = None + print("Warning : the crs is not stored in the output GeoDataFrame") + + return gdf + + +def test_gdf_zonal_stats(): + gpd = pytest.importorskip("geopandas") + polygons = os.path.join(DATA, "polygons.shp") + gdf = gpd.read_file(polygons) + if not hasattr(gdf, "__geo_interface__"): + pytest.skip("This version of geopandas doesn't support df.__geo_interface__") + + expected_with_path = gdf_zonal_stats(polygons, raster) + expected_with_gdf = gdf_zonal_stats(gdf, raster) + assert expected_with_gdf.equals(expected_with_path) # Use equals() to compare DataFrames + assert expected_with_gdf.crs == gdf.crs # Check that the crs is conserved + + # TODO # gen_zonal_stats() # TODO # gen_zonal_stats(stats=nodata) # TODO # gen_zonal_stats()