python – How to Calculate Mean Value of a Raster for Each Polygon in a Shapefile

geopandaspolygonpythonrastershapefile

I wrote a Python script who calculate, for rasters, mean value for each polygon of a shapefile. But I think that it's a "war machine" (it takes many treatment time) and I would like to simplify it.

Is somebody know a simple way to do it?

I'm using geopandas for shapefile manipulations.

Best Answer

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)

Related Solutions

[GIS] Get shapefile polygon attribute value at a specific point using Python (e.g. via PyShp)

At the simplest level, you must loop through each polygon in your shapefile and check to see if they point is inside it. You can stop the loop when you find the right polygon. You can speed this process up a little bit by checking if the point is inside the extent of the polygon. By checking the extent first you are quickly narrowing the possibilies by quickly checking the 4 points of a rectangle.

The following function checks if a point is inside a polygon. This function alone will answer your question. If you use the pyshp shapeRecords() method, it will keep the geometry and the records associated so when you find the right polygon, you'll have easy access to the attributes:

def point_in_poly(x,y,poly):

   # check if point is a vertex
   if (x,y) in poly: return "IN"

   # check if point is on a boundary
   for i in range(len(poly)):
      p1 = None
      p2 = None
      if i==0:
         p1 = poly[0]
         p2 = poly[1]
      else:
         p1 = poly[i-1]
         p2 = poly[i]
      if p1[1] == p2[1] and p1[1] == y and x > min(p1[0], p2[0]) and x < max(p1[0], p2[0]):
         return "IN"

   n = len(poly)
   inside = False

   p1x,p1y = poly[0]
   for i in range(n+1):
      p2x,p2y = poly[i % n]
      if y > min(p1y,p2y):
         if y <= max(p1y,p2y):
            if x <= max(p1x,p2x):
               if p1y != p2y:
                  xints = (y-p1y)*(p2x-p1x)/(p2y-p1y)+p1x
               if p1x == p2x or x <= xints:
                  inside = not inside
      p1x,p1y = p2x,p2y

   if inside: return "IN"
   else: return "OUT"

# Test a vertex for inclusion
poligono = [(-33.416032,-70.593016), (-33.415370,-70.589604),
(-33.417340,-70.589046), (-33.417949,-70.592351),
(-33.416032,-70.593016)]
lat= -33.416032
lon= -70.593016

print point_in_poly(lat, lon, poligono)

# test a boundary point for inclusion
poly2 = [(1,1), (5,1), (5,5), (1,5), (1,1)]
x = 3
y = 1
print point_in_poly(x, y, poly2)

Python GDAL – How to Generate Euclidean Distance for Each Polygon in Shapefile

Here is an example that

opens vector dataset and loops through each feature,
gdal.Rasterize each feature
gdal.ComputeProximity each rasterized feature

You should reproject your vector if they are in geographic lat/lon CRS so the proximity (distance) values are in metres, I used QGIS to reproject this data to UTM Zone 18N.

from osgeo import gdal
gdal.UseExceptions()

out_raster_template = "out_{}.tif"
shape_file = "boroughs_utm.shp"

pixel_size = 10
nodata = -9999

id_field = 'boro_code'

drv = gdal.GetDriverByName("ESRI Shapefile")

shp_ds = gdal.OpenEx(shape_file, gdal.OF_VECTOR)
lyr = shp_ds.GetLayer()

xmin, xmax, ymin, ymax = lyr.GetExtent()
srs = lyr.GetSpatialRef()

feat_def = lyr.GetLayerDefn()

lyr.ResetReading()
for feat in lyr:
    id = int(feat.GetField(id_field))

    tmp_feat = feat.Clone()

    out_raster = out_raster_template.format(id)
    tmp_fn = '/vsimem/tmp.shp'
    tmp_ds = drv.Create(tmp_fn, 0, 0, 0, gdal.GDT_Unknown )
    tmp_lyr = tmp_ds.CreateLayer(tmp_fn, None, feat_def.GetGeomType())
    tmp_lyr.CreateFeature(tmp_feat)
    tmp_feat, tmp_lyr, tmp_ds = None, None, None

    out_ds = gdal.Rasterize(out_raster, tmp_fn,
                   outputType=gdal.GDT_Float32, format='GTIFF', creationOptions=["COMPRESS=DEFLATE"],
                   noData=nodata, initValues=nodata,
                   xRes=pixel_size, yRes=-pixel_size, outputBounds=(xmin, ymin, xmax, ymax), outputSRS=srs,
                   allTouched=True, burnValues=0)

    out_ds = None

    src_ds = gdal.OpenEx(out_raster, gdal.OF_RASTER)
    dst_ds = gdal.OpenEx(out_raster, gdal.OF_UPDATE)

    src_band = src_ds.GetRasterBand(1)
    dst_band = dst_ds.GetRasterBand(1)

    gdal.ComputeProximity(src_band, dst_band, options=[f'VALUES=0'])

    dst_band, src_band, dst_ds, src_ds = None, None, None, None

    drv.Delete(tmp_fn)

If you want to rasterize some value from your shapefile to use later, but still do the euc distance to another:

from osgeo import gdal
gdal.UseExceptions()


#Create euclidean distance for each polygon and store "Values"
out_raster_template = "out_{}.tif"
out_proximity_template = "prox_{}.tif"
shape_file = "boroughs_utm.shp"

pixel_size = 10
nodata = -9999

id_field = 'boro_code'
value_field = 'value'

drv = gdal.GetDriverByName("ESRI Shapefile")

shp_ds = gdal.OpenEx(shape_file, gdal.OF_VECTOR)
lyr = shp_ds.GetLayer()

xmin, xmax, ymin, ymax = lyr.GetExtent()
srs = lyr.GetSpatialRef()

feat_def = lyr.GetLayerDefn()

lyr.ResetReading()
for feat in lyr:
    id = int(feat.GetField(id_field))
    val = feat.GetField(value_field)

    tmp_feat = feat.Clone()

    out_raster = out_raster_template.format(id)
    prox_raster = out_proximity_template.format(id)
    tmp_fn = '/vsimem/tmp.shp'
    tmp_raster = '/vsimem/tmp.tif'
    tmp_ds = drv.Create(tmp_fn, 0, 0, 0, gdal.GDT_Unknown )
    tmp_lyr = tmp_ds.CreateLayer(tmp_fn, None, feat_def.GetGeomType())
    tmp_lyr.CreateFeature(tmp_feat)
    tmp_feat, tmp_lyr, tmp_ds = None, None, None

    out_ds = gdal.Rasterize(out_raster, tmp_fn,
                            outputType=gdal.GDT_Float32, format='GTIFF', creationOptions=["COMPRESS=DEFLATE"],
                            noData=nodata, initValues=nodata,
                            xRes=pixel_size, yRes=-pixel_size, outputBounds=(xmin, ymin, xmax, ymax), outputSRS=srs,
                            allTouched=True, burnValues=val)

    out_ds = None

    out_ds = gdal.Rasterize(tmp_raster, tmp_fn,
                            outputType=gdal.GDT_Int32, format='GTIFF', creationOptions=["COMPRESS=DEFLATE"],
                            noData=nodata, initValues=nodata,
                            xRes=pixel_size, yRes=-pixel_size, outputBounds=(xmin, ymin, xmax, ymax), outputSRS=srs,
                            allTouched=True, burnValues=id)

    out_ds = None

    gdal.Translate(prox_raster, out_raster, creationOptions=["COMPRESS=DEFLATE"])
    src_ds = gdal.OpenEx(tmp_raster, gdal.OF_RASTER)
    dst_ds = gdal.OpenEx(prox_raster, gdal.OF_UPDATE)

    src_band = src_ds.GetRasterBand(1)
    dst_band = dst_ds.GetRasterBand(1)

    gdal.ComputeProximity(src_band, dst_band, options=[f'VALUES={id}'])

    dst_band, src_band, dst_ds, src_ds = None, None, None, None

    drv.Delete(tmp_fn)

Best Answer

Related Solutions

[GIS] Get shapefile polygon attribute value at a specific point using Python (e.g. via PyShp)

Python GDAL – How to Generate Euclidean Distance for Each Polygon in Shapefile

Related Question