标签：shp 裁剪 栅格 im GDAL geoTrans data gdal

　　本节将介绍如何在Python中用GDAL实现根据矢量边界裁剪栅格数据。

from osgeo import gdal, gdal_array
import shapefile
import numpy as np
import os

#批量shp裁剪tiff影像
try:
    import Image
    import ImageDraw
except:
    from PIL import Image, ImageDraw

def read_tiff(inpath):
  ds=gdal.Open(inpath)
  row=ds.RasterXSize
  col=ds.RasterYSize
  band=ds.RasterCount

  data=np.zeros([row,col,band])
  for i in range(band):
   dt=ds.GetRasterBand(1)
   data[:,:,i]=dt.ReadAsArray(0,0,col,row)
  return data

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

def world2Pixel(geoMatrix, x, y):
    """
    使用GDAL库的geomatrix对象((gdal.GetGeoTransform()))计算地理坐标的像素位置
    """
    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) / abs(yDist))
    return (pixel, line)


def write_img(filename,im_proj,im_geotrans,im_data):
    if 'int8' in im_data.dtype.name:
        datatype = gdal.GDT_Byte
    elif 'int16' in im_data.dtype.name:
        datatype = gdal.GDT_UInt16
    else:
        datatype = gdal.GDT_Float32

    if len(im_data.shape) == 3:
        im_bands, im_height, im_width = im_data.shape
    else:
        im_bands, (im_height, im_width) = 1,im_data.shape

    driver = gdal.GetDriverByName("GTiff")
    dataset = driver.Create(filename, im_width, im_height, im_bands, datatype)

    dataset.SetGeoTransform(im_geotrans)
    dataset.SetProjection(im_proj)
    if im_bands == 1:

        dataset.GetRasterBand(1).WriteArray(im_data)
    else:
        for i in range(im_bands):
            dataset.GetRasterBand(i+1).WriteArray(im_data[i])

    del dataset

def sha_raster(raster,shp,output):
    srcArray = gdal_array.LoadFile(raster)
    # 同时载入gdal库的图片从而获取geotransform
    srcImage = gdal.Open(raster)
    geoProj = srcImage.GetProjection()
    geoTrans = srcImage.GetGeoTransform()
    r = shapefile.Reader(shp)
    # 将图层扩展转换为图片像素坐标
    minX, minY, maxX, maxY = r.bbox
    ulX, ulY = world2Pixel(geoTrans, minX, maxY)
    lrX, lrY = world2Pixel(geoTrans, maxX, minY)
    pxWidth = int(lrX - ulX)
    pxHeight = int(lrY - ulY)
    clip = srcArray[:, ulY:lrY, ulX:lrX]
    # 为图片创建一个新的geomatrix对象以便附加地理参照数据
    geoTrans = list(geoTrans)
    geoTrans[0] = minX
    geoTrans[3] = maxY
    # 在一个空白的8字节黑白掩膜图片上把点映射为像元绘制市县
    # 边界线
    pixels = []
    for p in r.shape(0).points:
        pixels.append(world2Pixel(geoTrans, p[0], p[1]))
    rasterPoly = Image.new("L", (pxWidth, pxHeight), 1)
    # 使用PIL创建一个空白图片用于绘制多边形
    rasterize = ImageDraw.Draw(rasterPoly)
    rasterize.polygon(pixels, 0)
    # 使用PIL图片转换为Numpy掩膜数组
    mask = image2Array(rasterPoly)
    name = os.path.basename(raster).split(".tif")[0]
    outfile = output + "\\" + name+ "_cut.tif"  # 对输出文件命名
    # 根据掩膜图层对图像进行裁剪
    clip = gdal_array.numpy.choose(mask, (clip, 0)).astype(gdal_array.numpy.uint16)
    write_img(outfile, geoProj, geoTrans, clip)
    gdal.ErrorReset()

if __name__ == "__main__":
    raster = r'D:\test\裁剪实验\image\15.tif'
    # 用于裁剪的多边形shp文件
    shp = r'D:\test\裁剪实验\shp\2.shp'
    # 裁剪后的栅格数据
    output = r'D:\test\裁剪实验\out'

    #依据shp创建掩膜进行对tiff文件的裁剪
    sha_raster(raster,shp,output)

标签：shp,裁剪,栅格,im,GDAL,geoTrans,data,gdal
来源： https://www.cnblogs.com/suoyike1001/p/15164720.html

本站声明： 1. iCode9 技术分享网（下文简称本站）提供的所有内容，仅供技术学习、探讨和分享；
2. 关于本站的所有留言、评论、转载及引用，纯属内容发起人的个人观点，与本站观点和立场无关；
3. 关于本站的所有言论和文字，纯属内容发起人的个人观点，与本站观点和立场无关；
4. 本站文章均是网友提供，不完全保证技术分享内容的完整性、准确性、时效性、风险性和版权归属；如您发现该文章侵犯了您的权益，可联系我们第一时间进行删除；
5. 本站为非盈利性的个人网站，所有内容不会用来进行牟利，也不会利用任何形式的广告来间接获益，纯粹是为了广大技术爱好者提供技术内容和技术思想的分享性交流网站。