python - 使用 GDAL 裁剪数字高程模型

Question

因此，我正在 Arcgis 中开发一个 Arcpy 工具箱，该工具箱将 DEM 栅格文件纳入特定处理。但是我需要剪辑这些图像，因为原始图像太大且太长而无法处理。但问题是，Arcgis 裁剪工具改变了我当时无法使用的数据类型。

我开始寻找代码来做到这一点，看起来 GDAL 库可能有助于用 shapefile 剪辑 geotiff。以下是我为适应我的 1 通道 DEM 进行了一些小改动的代码：< https://pcjericks.github.io/py-gdalogr-cookbook/raster_layers.html >

from osgeo import gdal, gdalnumeric, ogr, osr
from PIL import Image, ImageDraw
gdal.UseExceptions()


# This function will convert the rasterized clipper shapefile
# to a mask for use within GDAL.
def imageToArray(i):
    """
    Converts a Python Imaging Library array to a
    gdalnumeric image.
    """
    a=gdalnumeric.fromstring(i.tostring(),'b')
    a.shape=i.im.size[1], i.im.size[0]
    return a

def arrayToImage(a):
    """
    Converts a gdalnumeric array to a
    Python Imaging Library Image.
    """
    i=Image.fromstring('L',(a.shape[1],a.shape[0]),
            (a.astype('b')).tostring())
    return i

def world2Pixel(geoMatrix, x, y):
  """
  Uses a gdal geomatrix (gdal.GetGeoTransform()) to calculate
  the pixel location of a geospatial coordinate
  """
  ulX = geoMatrix[0]
  ulY = geoMatrix[3]
  xDist = geoMatrix[1]
  yDist = geoMatrix[5]
  rtnX = geoMatrix[2]
  rtnY = geoMatrix[4]
  pixel = int((x - ulX) / xDist)
  line = 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 = None, xoff=0, yoff=0 ):
    ds = gdal.Open( gdalnumeric.GetArrayFilename(array) )

    if ds is not None and prototype_ds is not None:
        if type(prototype_ds).__name__ == 'str':
            prototype_ds = gdal.Open( prototype_ds )
        if prototype_ds is not None:
            gdalnumeric.CopyDatasetInfo( prototype_ds, ds, xoff=xoff, yoff=yoff )
    return ds

def histogram(a, bins=range(0,256)):
  """
  Histogram function for multi-dimensional array.
  a = array
  bins = range of numbers to match
  """
  fa = a.flat
  n = gdalnumeric.searchsorted(gdalnumeric.sort(fa), bins)
  n = gdalnumeric.concatenate([n, [len(fa)]])
  hist = n[1:]-n[:-1]
  return hist

def stretch(a):
  """
  Performs a histogram stretch on a gdalnumeric array image.
  """
  hist = histogram(a)
  im = arrayToImage(a)
  lut = []
  for b in range(0, len(hist), 256):
    # step size
    step = reduce(operator.add, hist[b:b+256]) / 255
    # create equalization lookup table
    n = 0
    for i in range(256):
      lut.append(n / step)
      n = n + hist[i+b]
  im = im.point(lut)
  return imageToArray(im)

def main( shapefile_path, raster_path ):
    # Load the source data as a gdalnumeric array
    srcArray = gdalnumeric.LoadFile(raster_path)

    # Also load as a gdal image to get geotransform
    # (world file) info
    srcImage = gdal.Open(raster_path)
    geoTrans = srcImage.GetGeoTransform()

    # Create an OGR layer from a boundary shapefile
    shapef = ogr.Open("%s.shp" % shapefile_path)
    lyr = shapef.GetLayer( os.path.split( os.path.splitext( shapefile_path )[0] )[1] )
    poly = lyr.GetNextFeature()

    # Convert the layer extent to image pixel coordinates
    minX, maxX, minY, maxY = lyr.GetExtent()
    ulX, ulY = world2Pixel(geoTrans, minX, maxY)
    lrX, lrY = world2Pixel(geoTrans, maxX, minY)

    # Calculate the pixel size of the new image
    pxWidth = int(lrX - ulX)
    pxHeight = int(lrY - ulY)

    clip = srcArray[ulY:lrY, ulX:lrX]

    #
    # EDIT: create pixel offset to pass to new image Projection info
    #
    xoffset =  ulX
    yoffset =  ulY
    print "Xoffset, Yoffset = ( %f, %f )" % ( xoffset, yoffset )

    # Create a new geomatrix for the image
    geoTrans = list(geoTrans)
    geoTrans[0] = minX
    geoTrans[3] = maxY

    # Map points to pixels for drawing the
    # boundary on a blank 8-bit,
    # black and white, mask image.
    points = []
    pixels = []
    geom = poly.GetGeometryRef()
    pts = geom.GetGeometryRef(0)
    for p in range(pts.GetPointCount()):
      points.append((pts.GetX(p), pts.GetY(p)))
    for p in points:
      pixels.append(world2Pixel(geoTrans, p[0], p[1]))
    rasterPoly = Image.new("L", (pxWidth, pxHeight), 1)
    rasterize = ImageDraw.Draw(rasterPoly)
    rasterize.polygon(pixels, 0)
    mask = imageToArray(rasterPoly)

    # Clip the image using the mask
    clip = gdalnumeric.choose(mask, \
        (clip, 0)).astype(gdalnumeric.uint8)


      clip[:,:] = stretch(clip[:,:])

    # Save new tiff
    #
    #  EDIT: instead of SaveArray, let's break all the
    #  SaveArray steps out more explicity so
    #  we can overwrite the offset of the destination
    #  raster
    #
    ### the old way using SaveArray
    #
    # gdalnumeric.SaveArray(clip, "OUTPUT.tif", format="GTiff", prototype=raster_path)
    #
    ###
    #
    gtiffDriver = gdal.GetDriverByName( 'GTiff' )
    if gtiffDriver is None:
        raise ValueError("Can't find GeoTiff Driver")
    gtiffDriver.CreateCopy( "OUTPUT.tif",
        OpenArray( clip, prototype_ds=raster_path, xoff=xoffset, yoff=yoffset )
    )

    # Save as an 8-bit jpeg for an easy, quick preview
    clip = clip.astype(gdalnumeric.uint8)
    gdalnumeric.SaveArray(clip, "OUTPUT.jpg", format="JPEG")

    gdal.ErrorReset()


if __name__ == '__main__':


    main( "shapefile", "DEM.tiff" )

但是我得到了一个“形状不匹配 ValueError”：

<ipython-input-22-32e4e8197a02> in main(shapefile_path, raster_path, region_shapefile_path)
    166 
    167     # Clip the image using the mask
--> 168     clip = gdalnumeric.choose(mask,         (clip, 0)).astype(gdalnumeric.uint8)
    169 


/home/edgar/anaconda3/envs/gis2/lib/python2.7/site-packages/numpy/core/fromnumeric.pyc in choose(a, choices, out, mode)
    399 
    400     """
--> 401     return _wrapfunc(a, 'choose', choices, out=out, mode=mode)
    402 
    403 

/home/edgar/anaconda3/envs/gis2/lib/python2.7/site-packages/numpy/core/fromnumeric.pyc in _wrapfunc(obj, method, *args, **kwds)
     49 def _wrapfunc(obj, method, *args, **kwds):
     50     try:
---> 51         return getattr(obj, method)(*args, **kwds)
     52 
     53     # An AttributeError occurs if the object does not have

ValueError: shape mismatch: objects cannot be broadcast to a single shape

我试图在代码中查看它可能来自哪里，我意识到这部分可能无法正常工作：

minX, maxX, minY, maxY = lyr.GetExtent()
ulX, ulY = world2Pixel(geoTrans, minX, maxY)
lrX, lrY = world2Pixel(geoTrans, maxX, minY)
print("ulX, lrX, ulY, lrY : " , ulX, lrX, ulY, lrY) #image pixel coordinates of the shapefile
print(srcArray.shape)  #shape of the raster image 
clip = srcArray[ ulY:lrY, ulX:lrX] #extracting the shapefile zone from the raster image?
print(clip)

并返回：

('ulX, lrX, ulY, lrY : ', 35487, 37121, 3844, 5399)
(5041, 5041)
[]

似乎这些索引超出了界限（但奇怪的是 python 并没有那么麻烦）并且没有任何内容被复制。因此，我尝试通过使用与我的总光栅图像相对应的 shapefile 来更改代码以获取与我希望提取的区域相对应的“真实”像素值：

#shapefile corresponding to the whole raster image
region_shapef = ogr.Open("%s.shp" % region_shapefile_path)
region_lyr = region_shapef.GetLayer( os.path.split( os.path.splitext( region_shapefile_path )[0] )[1] )

RminX, RmaxX, RminY, RmaxY = region_lyr.GetExtent()
RulX, RulY = world2Pixel(geoTrans, RminX, RmaxY)
RlrX, RlrY = world2Pixel(geoTrans, RmaxX, RminY)

#linear regression to find the equivalent pixel values of the clipping zone
pX = float(srcArray.shape[1])/(RlrX - RulX)
X0 = -(RulX*pX)
pY = float(srcArray.shape[0])/(RlrY - RulY)
Y0 = -(RulY*pY)

idXi = int(ulX*pX+X0)
idXf = int(lrX*pX+X0)
idYi = int(ulY*pY+Y0)
idYf = int(lrY*pY+Y0)

clip = srcArray[idYi:idYf, idXi:idXf]
print(clip)

返回一个真正提取值的数组：

[[169.4 171.3 173.7 ... 735.6 732.8 729.7]
 [173.3 176.4 179.9 ... 734.3 731.5 728.7]
 [177.8 182.  186.5 ... 733.1 730.3 727.5]
 ...
 [ 73.3  77.5  83.  ... 577.4 584.9 598.1]
 [ 72.8  76.5  81.5 ... 583.1 593.  606.2]
 [ 71.3  74.7  79.  ... 588.9 599.1 612.3]]

虽然我还有那个该死的：

<ipython-input-1-d7714555354e> in main(shapefile_path, raster_path, region_shapefile_path)
    170 
    171     # Clip the image using the mask
--> 172     clip = gdalnumeric.choose(mask,         (clip, 0)).astype(gdalnumeric.uint8)
    173 
    174     # This image has 3 bands so we stretch each one to make them

/home/edgar/anaconda3/envs/gis2/lib/python2.7/site-packages/numpy/core/fromnumeric.pyc in choose(a, choices, out, mode)
    399 
    400     """
--> 401     return _wrapfunc(a, 'choose', choices, out=out, mode=mode)
    402 
    403 

/home/edgar/anaconda3/envs/gis2/lib/python2.7/site-packages/numpy/core/fromnumeric.pyc in _wrapfunc(obj, method, *args, **kwds)
     49 def _wrapfunc(obj, method, *args, **kwds):
     50     try:
---> 51         return getattr(obj, method)(*args, **kwds)
     52 
     53     # An AttributeError occurs if the object does not have

ValueError: shape mismatch: objects cannot be broadcast to a single shape

我是否遗漏或误解了它？我真的开始缺乏想法，所以如果有人有想法，那将不胜感激。

否则，如果您知道另一种方法来剪辑我的 DEM 而不改变它，我也很好。

标签： pythonarcgisgdalshapefile

Answer 1

您可以使用 gdal.Warp() 和 shapefile 作为切割线以更简单的方式实现这一点

from osgeo import gdal

input_raster = "path/to/yourDEM.tif" 
# or as an alternative if the input is already a gdal raster object you can use that gdal object
input_raster=gdal.Open("path/to/yourDEM.tif")
input_shape = "path/to/yourShapefile.shp" # or any other format
output_raster="path/to/outputDEM.tif" #your output raster file

ds = gdal.Warp(output_raster,
              input_raster,
              format = 'GTiff',
              cutlineDSName = input_shape, # or any other file format
              cutlineWhere="FIELD = 'whatever'" # optionally you can filter your cutline (shapefile) based on attribute values
              dstNodata = -9999) # select the no data value you like
ds=None #do other stuff with ds object, it is your cropped dataset. in this case we only close the dataset.