python随笔二：批量栅格图像的裁剪（单独的应该也没问题）

from osgeo import ogr 
#这个就是前面说的gdal相关的包，前面直接用的import，就一直加载不上，但是换成这个后就可以加载了
import os
import numpy
from tqdm import tqdm
from osgeo import gdal
import time
from osgeo import gdalnumeric
from PIL import Image,ImageDraw

# 应该是使gdal库支持中文编码、中文字段、中文路径
gdal.SetConfigOption("GDAL_ARRAY_OPEN_BY_FILENAME","TRUE")
# 栅格文件夹的路径
rasterpath=r"E:\NDVI\青藏高原NDVI1982_2020逐月数据5km\result"
# 栅格后缀名的设置(.dat/.tif)
lastname=".tif"
#矢量文件的路径
shp=r"E:\QTPSHP\DBATP_Polygon.shp"
#命名字段(矢量文件的字段)，矢量文件的名字
filename="DBATP_Polygon.shp"
#裁剪后文件存放位置（文件夹路径）
outputpath = r"E:\NDVI\青藏高原NDVI1982_2020逐月数据5km\cut"

#将一个Python图像库的数组转换为一个gdal_array图片
def image2Array(i):
    a=gdalnumeric.fromstring(i.tobytes(),'b')
    a.shape=i.im.size[1],i.im.size[0]
    return a

# 数组写入dataset
def OpenArray(array,prototype_ds = None,xoff=0,yoff=0):
    ds = gdal.Open(gdalnumeric.GetArrayFilename(array))
    print(ds)
    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 world2Pixel(geoMatrix, x, y):
    ulx = geoMatrix[0]
    uly = geoMatrix[3]
    xDist = geoMatrix[1]
    yDist = geoMatrix[5]
    pixel = int((x - ulx) // xDist)
    line = int((uly - y) // abs(yDist))
    return (pixel, line)

strattime=time.time()

# 读取栅格
rasters = os.listdir(rasterpath)
rasterlist=list(filter(lambda x: x[-4:] == lastname, rasters))
print(rasterlist)
# 裁剪
with tqdm(total=len(rasterlist), iterable='iterable') as pbar:
    for ra in rasterlist:
        raster = str(rasterpath + '/' + ra)
        print(raster)

        # 将数据源作为gdal_array载入
        srcArray=gdalnumeric.LoadFile(raster)
        # 同时载入gdal库的图片从而获取geotransform
        srcImage=gdal.Open(raster)
        geoTrans=srcImage.GetGeoTransform()
        proj = srcImage.GetProjection()

        # 打开前面设置的shp文件，下面的也是一样的，就是获取shape文件
        r = ogr.Open(shp)
        lyr = r.GetLayer()
        # 获取要素
        feature=lyr.GetNextFeature()
        while feature:
            geometry=feature.geometry()
            name = feature.GetField("id")
            # 这个name的设置可能只是得到一个field对象，之后继续往后输出
            print(name)
            # 将图层扩展转换为图片像素坐标，需要每一个shp点的所在像素的左上角坐标
            minX,maxX,minY,maxY = geometry.GetEnvelope()

            # 计算要素四至对应的图片四至（左上、右下经纬度）
            ulX, ulY = world2Pixel(geoTrans, minX, maxY)
            lrX, lrY = world2Pixel(geoTrans, maxX, minY)
            # 计算新图片的尺寸
            pxWidth = lrX - ulX
            pxHeight = lrY - ulY
            # 获取新图片影像数组
            clip = srcArray[ulY: lrY, ulX: lrX]
            print("clip.shape",clip.shape)
            # 计算新图片四至
            newOriginX = geoTrans[0]+ulX*geoTrans[1]
            newOriginY = geoTrans[3]+ulY*geoTrans[5]
            newEndX = geoTrans[0]+lrX*geoTrans[1]
            newEndY = geoTrans[3]+lrY*geoTrans[5]
            # 为新图片创建一个新的geomatrix对象以便附加地理参照数据
            newgeoTrans=list(geoTrans)
            newgeoTrans[0]=newOriginX
            newgeoTrans[3]=newOriginY

            # 使用PIL创建一个空白图片用于绘制多边形
            rasterPoly = Image.new('L', (pxWidth, pxHeight), 1)
            draw=ImageDraw.Draw(rasterPoly)
            # 在一个空白的8字节黑白掩膜图片上把点映射为像元绘制要素
            points = []
            pixels = []
            geom = feature.GetGeometryRef()
            cnt = geom.GetGeometryCount()

            if cnt > 1:#一个要素有多个元素(例如中国shp上的大陆与海南岛)
                for n in range(0,cnt):
                    pos = geom.GetGeometryRef(n)#依次获取要素的部分元素polygon
                    pts = pos.GetGeometryRef(0)# 获取linestring
                    print("Geometry:", n, "points:", pts.GetPointCount())
                    for i in range(pts.GetPointCount()):#获取point
                        points.append((pts.GetX(i), pts.GetY(i)))
    #                 print(len(points))
                    for p in points:
                        p1,p2 = world2Pixel(newgeoTrans, p[0], p[1])
                        pixels.append((p1,p2))
    #                 print(len(pixels))
                    # 用像元位置绘制多边形
                    draw.polygon(pixels,fill=0)
                    print("Geometry", n, "PILsize:",draw)
                    points = []
                    pixels = []

            else:  # 要素只有一个元素
                pts = geom.GetGeometryRef(0)
                # 获取所有点
                for i in range(pts.GetPointCount()):
                    points.append((pts.GetX(i), pts.GetY(i)))
                # 点坐标转图片坐标
                for p in points:
                    p1, p2 = world2Pixel(newgeoTrans, p[0], p[1])
                    pixels.append((p1, p2))
                # print(len(pixels))
                # 用像元位置绘制多边形
                draw.polygon(pixels, fill=0)
                print("PILsize:", draw)

            # 使用PIL图片转换为Numpy掩膜数组
            mask = numpy.array(rasterPoly)
            # 根据掩膜图层对图像进行裁剪
            clip = gdalnumeric.numpy.choose(mask,(clip,0)).astype(gdalnumeric.numpy.float32)

            # 保存裁剪后的成果
            gtiffDriver = gdal.GetDriverByName('GTiff')
            if gtiffDriver is None:
                raise ValueError("Can't found Geotiff Driver")
# output这里卡了很久，提示的TypeError: can only concatenate str (not "int") to str，这个的意思是字符的类别不对，应该是使用字符串对应字符串才对，所以后面name本来没有str这个的，是在后面加上的，这样name的字符类型就和前面的字符类型对上了，就都是字符串类型，就可以运行了。
            output = outputpath+(raster.split("/")[-1]).split(lastname)[0]+'_'+str(name)+"_clip.tiff"
            # 数组写入dataset
            out = OpenArray(clip, prototype_ds=raster, xoff=ulX, yoff=ulY)
            print("out:", out)
            # 保存dataset
            outds = gtiffDriver.CreateCopy(output, out)
            # 写入坐标
            outds.SetGeoTransform(newgeoTrans)
            # 设置投影
            outds.SetProjection(proj)
            del out, outds, geometry, clip, mask
            feature = lyr.GetNextFeature()

        del feature, srcArray, srcImage, geoTrans, proj, r, lyr
        pbar.update(1)

endtime = time.time()
print("spend:", endtime-strattime)

 还有很多代码的意思都不懂，只是会运行，但是很多都还需要继续学才行，要逐渐积累，积累的够多，后面才会进步的够快。

原代码链接：(23条消息) 批量裁剪栅格 python gdal_一只大笨猪的博客-CSDN博客