A basic introduction to Raster Zonal Stats
Basic Usage
Generate zonal stats for a GeoDataframe containing areas of interest using raster data
Terms:
- aoi - (area of interest) a geodataframe which we are interested in generating zonal statistics for
- raster data - the source raster containing the features which we are interested in collecting zonal stats for our aoi.
import geopandas as gpd
import matplotlib.pyplot as plt
import geowrangler.raster_zonal_stats as rzs
%%time
# area multipolygons for regions 3,4,ncr of the philippines
aoi = gpd.read_file("../data/region34ncr_admin.geojson")
ax = aoi.plot(
ax=plt.axes(),
facecolor="none",
edgecolor=[
"#C62828",
"#283593",
"#FF9800",
],
)
aoi
Download Philippine Population Data
see the Humanitarian Data Exchange World Population Counts - Philippines
We download our raster Data as GeoTiff files from the Humanitarian Data Exchange site.
Note: This maybe slow as the file is about 180 Mb and depending on your internet download speed may take more than 5 minutes
phil_pop_link = "https://data.worldpop.org/GIS/Population/Global_2000_2020/2020/PHL/phl_ppp_2020.tif"
phil_pop_dset = "phl_pop_2020.tif"
%%time
![ ! -e ../data/{phil_pop_dset} ] && curl -o ../data/{phil_pop_dset} {phil_pop_link}
To create our raster zonal stats, we just need to
set the aggregations, as well as some extra arguments,
such as the nodata value in the raster.
%%time
results = rzs.create_raster_zonal_stats(
aoi,
f"../data/{phil_pop_dset}",
aggregation=dict(
func=["sum", "count"],
column="population",
output=["population_count", "samples"],
),
extra_args=dict(nodata=-99999), # nodata value is -99999
)
results
grid_aoi_file = "../data/region3_admin_grids.geojson"
%%time
grid_aoi_results = rzs.create_raster_zonal_stats(
grid_aoi_file,
f"../data/{phil_pop_dset}",
aggregation=dict(
func=["sum", "count"],
column="population",
output=["population_count", "samples"],
fillna=[True, True],
),
extra_args=dict(nodata=-99999), # nodata value is -99999
)
grid_aoi_results.head()
ax = aoi[aoi.Reg_Name == "Region III"].plot(
ax=plt.axes(), facecolor="none", edgecolor="black"
)
ax = grid_aoi_results.plot(
ax=ax, column="population_count", edgecolor="blue", alpha=0.5
)
bingtile_grid_aoi_file = "../data/region3_bingtile_grid13.geojson"
%%time
bingtile_grid_aoi_results = rzs.create_raster_zonal_stats(
bingtile_grid_aoi_file,
f"../data/{phil_pop_dset}",
aggregation=dict(
func=["sum", "count"],
column="population",
output=["population_count", "samples"],
fillna=[True, True],
),
extra_args=dict(nodata=-99999), # nodata value is -99999
)
bingtile_grid_aoi_results.head()
ax = aoi[aoi.Reg_Name == "Region III"].plot(
ax=plt.axes(), facecolor="none", edgecolor="black"
)
ax = bingtile_grid_aoi_results.plot(
ax=ax, column="population_count", edgecolor="blue", alpha=0.5
)
