import gdal, ogr, osr, numpy import sys def zonal_stats(feat, input_zone_polygon, input_value_raster): # Open data raster = gdal.Open(input_value_raster) shp = ogr.Open(input_zone_polygon) lyr = shp.GetLayer() # Get raster georeference info transform = raster.GetGeoTransform() xOrigin = transform[0] yOrigin = transform[3] pixelWidth = transform[1] pixelHeight = transform[5] # Reproject vector geometry to same projection as raster sourceSR = lyr.GetSpatialRef() targetSR = osr.SpatialReference() targetSR.ImportFromWkt(raster.GetProjectionRef()) coordTrans = osr.CoordinateTransformation(sourceSR,targetSR) feat = lyr.GetNextFeature() geom = feat.GetGeometryRef() geom.Transform(coordTrans) # Get extent of feat geom = feat.GetGeometryRef() if (geom.GetGeometryName() == 'MULTIPOLYGON'): count = 0 pointsX = []; pointsY = [] for polygon in geom: geomInner = geom.GetGeometryRef(count) ring = geomInner.GetGeometryRef(0) numpoints = ring.GetPointCount() for p in range(numpoints): lon, lat, z = ring.GetPoint(p) pointsX.append(lon) pointsY.append(lat) count += 1 elif (geom.GetGeometryName() == 'POLYGON'): ring = geom.GetGeometryRef(0) numpoints = ring.GetPointCount() pointsX = []; pointsY = [] for p in range(numpoints): lon, lat, z = ring.GetPoint(p) pointsX.append(lon) pointsY.append(lat) else: sys.exit("ERROR: Geometry needs to be either Polygon or Multipolygon") xmin = min(pointsX) xmax = max(pointsX) ymin = min(pointsY) ymax = max(pointsY) # Specify offset and rows and columns to read xoff = int((xmin - xOrigin)/pixelWidth) yoff = int((yOrigin - ymax)/pixelWidth) xcount = int((xmax - xmin)/pixelWidth)+1 ycount = int((ymax - ymin)/pixelWidth)+1 # Create memory target raster target_ds = gdal.GetDriverByName('MEM').Create('', xcount, ycount, 1, gdal.GDT_Byte) target_ds.SetGeoTransform(( xmin, pixelWidth, 0, ymax, 0, pixelHeight, )) # Create for target raster the same projection as for the value raster raster_srs = osr.SpatialReference() raster_srs.ImportFromWkt(raster.GetProjectionRef()) target_ds.SetProjection(raster_srs.ExportToWkt()) # Rasterize zone polygon to raster gdal.RasterizeLayer(target_ds, [1], lyr, burn_values=[1]) # Read raster as arrays banddataraster = raster.GetRasterBand(1) dataraster = banddataraster.ReadAsArray(xoff, yoff, xcount, ycount).astype(numpy.float) bandmask = target_ds.GetRasterBand(1) datamask = bandmask.ReadAsArray(0, 0, xcount, ycount).astype(numpy.float) # Mask zone of raster zoneraster = numpy.ma.masked_array(dataraster, numpy.logical_not(datamask)) # Calculate statistics of zonal raster return numpy.average(zoneraster),numpy.mean(zoneraster),numpy.median(zoneraster),numpy.std(zoneraster),numpy.var(zoneraster) def loop_zonal_stats(input_zone_polygon, input_value_raster): shp = ogr.Open(input_zone_polygon) lyr = shp.GetLayer() featList = range(lyr.GetFeatureCount()) statDict = {} for FID in featList: feat = lyr.GetFeature(FID) meanValue = zonal_stats(feat, input_zone_polygon, input_value_raster) statDict[FID] = meanValue return statDict def main(input_zone_polygon, input_value_raster): return loop_zonal_stats(input_zone_polygon, input_value_raster) if __name__ == "__main__": # # Returns for each feature a dictionary item (FID) with the statistical values in the following order: Average, Mean, Medain, Standard Deviation, Variance # # example run : $ python grid.py <full-path><output-shapefile-name>.shp xmin xmax ymin ymax gridHeight gridWidth # if len( sys.argv ) != 3: print "[ ERROR ] you must supply two arguments: input-zone-shapefile-name.shp input-value-raster-name.tif " sys.exit( 1 ) print 'Returns for each feature a dictionary item (FID) with the statistical values in the following order: Average, Mean, Medain, Standard Deviation, Variance' print main( sys.argv[1], sys.argv[2] )
此方法计算一个矢量数据集内的栅格的值的统计数据。它返回的每个特征词典项目(FID)依次为:平均统计值,均值,标准差,方差Medain,
而这个配方的作品,是一个很好的例子,通常建议使用[ rasterstats ](https://github.com/perrygeo/python-raster-stats)计算分区统计Python。
关注公众号
获取免费资源
Copyright © Since 2014.
开源地理空间基金会中文分会
吉ICP备05002032号
Powered by TorCMS