利用Python中的GDAL和OGR模块实现shapefile对栅格DEM数据的裁剪

之前已经把31景DEM数据拼接完成，这篇文章，利用shapefile文件实现对山东DEM数据的裁剪。

参考：https://blog.csdn.net/GISShiXiSheng/article/details/72843070?ops_request_misc=%257B%2522request%255Fid%2522%253A%2522161423692216780271562057%2522%252C%2522scm%2522%253A%252220140713.130102334.pc%255Fblog.%2522%257D&request_id=161423692216780271562057&biz_id=0&utm_medium=distribute.pc_search_result.none-task-blog-2~blog~first_rank_v1~rank_blog_v1-1-72843070.pc_v1_rank_blog_v1&utm_term=%E6%A0%85%E6%A0%BC%E6%95%B0%E6%8D%AE%E8%A3%81%E5%89%AA&spm=1018.2226.3001.4450

Python、GIS基础教程：https://www.osgeo.cn/pygis/ogr-ogrintro.html

1.代码

from osgeo import gdal, gdalnumeric, ogr,gdal_array
from PIL import Image, ImageDraw
import os
import operator

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.tobytes(), '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.frombytes('L', (a.shape[1], a.shape[0]),
                        (a.astype('b')).tobytes())
    return i


def world2Pixel(geoMatrix, x, y):
    """
    Uses a gdal geomatrix (gdal.GetGeoTransform()) to calculate
    the pixel location of a geospatial coordinate
    """
    # GetGeoTransform返回值：左上角x坐标， 水平分辨率，旋转参数， 左上角y坐标，旋转参数，竖直分辨率。

    ulX = geoMatrix[0]
    ulY = geoMatrix[3]
    xDist = geoMatrix[1]  # 水平分辨率
    pixel = int((x - ulX) / xDist)
    line = int((ulY - y) / xDist)
    return (pixel, line)


#
# EDIT: this is basically an overloaded
# version of the gdal_array.