Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- 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] )
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement