rasterFunc.py 文件源码

def filterByCoverage(vectorFile, rasterFile, covPerc):

    srcVector = ogr.Open(vectorFile)
    srcLayer = srcVector.GetLayer()
    # merge all features in one geometry (multi polygone)
    multi  = ogr.Geometry(ogr.wkbMultiPolygon)
    for feature in srcLayer:
        geom = feature.GetGeometryRef()
        multi.AddGeometry(geom)

    #attributes of raster file
    rasList = raster2array(rasterFile)

    xsize = rasList[4][0]
    ysize = abs(rasList[4][1])

    pixel_area = xsize*ysize

    rows = rasList[0].shape[0]
    cols = rasList[0].shape[1]

    x1 = rasList[2][0]
    y1 = rasList[2][3]

    #iterate over raster cells
    for i in range(rows):
        for j in range(cols):
            ring = ogr.Geometry(ogr.wkbLinearRing)

            ring.AddPoint(x1, y1)
            ring.AddPoint(x1 + xsize, y1)
            ring.AddPoint(x1 + xsize, y1 - ysize)
            ring.AddPoint(x1, y1 - ysize)
            ring.AddPoint(x1, y1)

            poly = ogr.Geometry(ogr.wkbPolygon)
            poly.AddGeometry(ring)

            intersect = multi.Intersection(poly)

            if intersect.ExportToWkt() != 'GEOMETRYCOLLECTION EMPTY':
                perc = (intersect.GetArea()/pixel_area)*100
                if perc > covPerc:
                    rasList[0][i][j] = numpy.nan     
            x1 += xsize
        x1 = rasList[2][0]
        y1 -= ysize

    return(rasList[0]) #return the filtered array


# numpy array to geo raster