基于gdal的矢量裁剪栅格

矢量与栅格

# -*- coding: utf-8 -*-
import os
from osgeo import gdal, gdalnumeric, ogr, osr, gdal_array
gdal.UseExceptions()def world2Pixel(geoMatrix, x, y):"""Uses a gdal geomatrix (gdal.GetGeoTransform()) to calculatethe pixel location of a geospatial coordinate"""ulX &＃61; geoMatrix[0]ulY &＃61; geoMatrix[3]xDist &＃61; geoMatrix[1]yDist &＃61; geoMatrix[5]rtnX &＃61; geoMatrix[2]rtnY &＃61; geoMatrix[4]pixel &＃61; int((x - ulX) / xDist)line &＃61; int((ulY - y) / xDist)return (pixel, line)#
# EDIT: this is basically an overloaded
# version of the gdal_array.OpenArray passing in xoff, yoff explicitly
# so we can pass these params off to CopyDatasetInfo
#
def OpenArray( array, prototype_ds &＃61; None, xoff&＃61;0, yoff&＃61;0 ):# ds &＃61; gdal.Open( gdalnumeric.GetArrayFilename(array))ds &＃61; gdal_array.OpenArray(array)if ds is not None and prototype_ds is not None:if type(prototype_ds).__name__ &＃61;&＃61; &＃39;str&＃39;:prototype_ds &＃61; gdal.Open( prototype_ds )if prototype_ds is not None:gdalnumeric.CopyDatasetInfo( prototype_ds, ds, xoff&＃61;xoff, yoff&＃61;yoff )return dsdef write_img(filename,im_proj,im_geotrans,im_data):if &＃39;int8&＃39; in im_data.dtype.name:datatype &＃61; gdal.GDT_Byteelif &＃39;int16&＃39; in im_data.dtype.name:datatype &＃61; gdal.GDT_UInt16else:datatype &＃61; gdal.GDT_Float32if len(im_data.shape) &＃61;&＃61; 3:im_bands, im_height, im_width &＃61; im_data.shapeelse:im_bands, (im_height, im_width) &＃61; 1,im_data.shapedriver &＃61; gdal.GetDriverByName("GTiff")dataset &＃61; driver.Create(filename, im_width, im_height, im_bands, datatype)dataset.SetGeoTransform(im_geotrans)dataset.SetProjection(im_proj)if im_bands &＃61;&＃61; 1:dataset.GetRasterBand(1).WriteArray(im_data)else:for i in range(im_bands):dataset.GetRasterBand(i&＃43;1).WriteArray(im_data[i])del datasetdef shpClipRaster(shapefile_path, raster_path, save_path):# Load the source data as a gdalnumeric array# srcArray &＃61; gdalnumeric.LoadFile(raster_path)# Also load as a gdal image to get geotransform# (world file) infosrcImage &＃61; gdal.Open(raster_path)geoTrans &＃61; srcImage.GetGeoTransform()geoProj &＃61; srcImage.GetProjection()# Create an OGR layer from a boundary shapefileshapef &＃61; ogr.Open(shapefile_path)lyr &＃61; shapef.GetLayer( os.path.split( os.path.splitext( shapefile_path )[0] )[1] )poly &＃61; lyr.GetNextFeature()# Convert the layer extent to image pixel coordinatesminX, maxX, minY, maxY &＃61; lyr.GetExtent()ulX, ulY &＃61; world2Pixel(geoTrans, minX, maxY)lrX, lrY &＃61; world2Pixel(geoTrans, maxX, minY)# Calculate the pixel size of the new imagepxWidth &＃61; int(lrX - ulX)pxHeight &＃61; int(lrY - ulY)# clip &＃61; srcArray[:, ulY:lrY, ulX:lrX]clip &＃61; srcImage.ReadAsArray(ulX,ulY,pxWidth,pxHeight) #***只读要的那块***## EDIT: create pixel offset to pass to new image Projection info#xoffset &＃61; ulXyoffset &＃61; ulYprint("Xoffset, Yoffset &＃61; ( %f, %f )" % ( xoffset, yoffset ))# Create a new geomatrix for the imagegeoTrans &＃61; list(geoTrans)geoTrans[0] &＃61; minXgeoTrans[3] &＃61; maxYwrite_img(save_path, geoProj, geoTrans, clip)gdal.ErrorReset()if __name__ &＃61;&＃61; "__main__":shp &＃61; r"D:\2021\6\ShanDong\00000001_V1\00000001_V1_POLY.shp"img &＃61; r"D:\2021\6\ShanDong\data\jiyang.tif"out &＃61; r"D:\2021\6\ShanDong\temp.tif"shpClipRaster(shp, img, out)

结果展示