rasterFunc.py 文件源码

def clipRaster(raster, newRaster, vector):
    raster = gdal.Open(raster)

    vect = ogr.Open(vector)
    lyr = vect.GetLayer()
    ext = lyr.GetExtent()

    gTrans = raster.GetGeoTransform()
    #get the x start of the left most pixel
    xlP = int((ext[0] - gTrans[0])/gTrans[1])*gTrans[1] - abs(gTrans[0])
    #get the x end of the right most pixel
    xrP = math.ceil((ext[1] - gTrans[0])/gTrans[1])*gTrans[1] - abs(gTrans[0])
    #get the y start of the upper most pixel
    yuP = abs(gTrans[3]) - int((gTrans[3] - ext[3])/gTrans[5])*gTrans[5]
    #get the y end of the lower most pixel
    ylP = abs(gTrans[3]) - math.floor((gTrans[3] - ext[2])/gTrans[5])*gTrans[5]

    gdal.Translate('tmp.tif', raster, projWin = [xlP, yuP, xrP, ylP])
    del raster
    tRas = gdal.Open('tmp.tif')
    band = tRas.GetRasterBand(1)
    noDat = band.GetNoDataValue()
    # store a copy before rasterize
    fullRas = band.ReadAsArray().astype(numpy.float)

    gdal.RasterizeLayer(tRas, [1], lyr, None, None, [-9999], ['ALL_TOUCHED=TRUE']) # now tRas is modified

    finRas = tRas.GetRasterBand(1).ReadAsArray().astype(numpy.float)

    for i, row in enumerate(finRas):
        for j, col in enumerate(row):
            if col == -9999.:
                finRas[i, j] = fullRas[i, j]
            else:
                finRas[i, j] = noDat

    array2raster(newRaster, 'tmp.tif', finRas, noDat, "float32")
    os.remove('tmp.tif')
    del fullRas, finRas, band, tRas

# create a reference raster with random values