I found the problem....
The problem is in defining the geo_transform. I had the following:
x_rotation = 0
y_rotation = 0
cell_width_meters = 50
cell_height_meters = 50
geo_transform = [ top_left_x, cell_width_meters, x_rotation, top_left_y, y_rotation, cell_height_meters ]
dataset.SetGeoTransform(geo_transform)
The Gdal documentation isn't real clear about what these values are. (See SetGeoTransform) Searching around the internets I derived that the passed values should be (in order):
- top_left_x
- cell_width_meters
- x_rotation
- top_left_y
- y_rotation
- cell_height_meters
Which seems right, BUT re-reviewing the GDAL API Tutorial I noticed that the last value, cell_height_meters was shown being given in a negative value. It seems that this was all that was needed to properly output the data in the expected orientation.
So now I've changed the geo_transform definition line to:
(Notice the added "-")
geo_transform = [ top_left_x, cell_width_meters, x_rotation, top_left_y, y_rotation, -cell_height_meters ]
You can use the rasterstats package for this.
For example (assuming you have a geopandas.GeoDataFrame called gdf):
from rasterstats import zonal_stats
with rasterio.open("/path/to/raster.tif") as src:
affine = src.transform
array = src.read(1)
df_zonal_stats = pd.DataFrame(zonal_stats(gdf, array, affine=affine))
# adding statistics back to original GeoDataFrame
gdf2 = pd.concat([gdf, df_zonal_stats], axis=1)
Best Answer
I believe that a call to "data.frame" is all you need. If you have several rasters, just stack them. If you want the rasters as columns and the statistic as rows you can just transpose on the fly (see second example)