【GDAL应用】按矢量数据的空间四至范围裁剪影像

RSer_gis

已于 2024-05-06 00:26:08 修改

阅读量503

点赞数

分类专栏： GDAL 文章标签： python gdal

于 2022-07-05 00:00:37 首次发布

本文链接：https://blog.youkuaiyun.com/m0_37609748/article/details/124581073

版权

GDAL 专栏收录该内容

5 篇文章

订阅专栏

本文介绍了一种使用Python GDAL库按矢量数据（如shp文件）的空间范围裁剪栅格影像的方法。该方法适用于线要素和面要素数据，并确保裁剪后的影像与原数据坐标系一致。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

1 实现效果

数据准备： 栅格数据、矢量数据
在这里插入图片描述
裁剪效果：

2 实现功能

按图层的空间范围（shp数据的空间范围四至）进行裁剪，该方法无论是线要素还是面要素数据都能按照图层范围四至裁剪影像。

3 实现代码

# -*- coding: utf-8 -*-
"""
@time: 2022-8-21 12:30
@author: RSer_gis
"""

import os
from osgeo import gdal, ogr
import re
import sys

gdal.UseExceptions()


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]
    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 shpClipRaster(shapefile_path, raster_path, save_path):
    # open a gdal image to get geotransform (world file) info
    srcImage = gdal.Open(raster_path)
    geoTrans = srcImage.GetGeoTransform()
    geoProj = srcImage.GetProjection()

    # Create an OGR layer from a boundary shapefile
    shapef = ogr.Open(shapefile_path)
    lyr = shapef.GetLayer(os.path.split(os.path.splitext(shapefile_path)[0])[1])
    # poly = lyr.GetNextFeature()
    spatial = lyr.GetSpatialRef()
    # 判断数输入数据坐标系是否一致
    # proj = osr.SpatialReference(wkt=geoProj)
    # epsg = proj.GetAttrValue('AUTHORITY', 1)
    # print(epsg)
    # epsg1 = spatial.GetAttrValue('AUTHORITY', 1)
    # print(epsg1)
    geoProj1 = geoProj.replace(" ", "")
    # print(geoProj1)
    spt = str(spatial)
    spt = spt.replace("\n", "")
    spt = spt.replace(" ", "")
    # print(spt)
    # print(geoProj1[0:6], spt[0:6])
    print(geoProj1[7:44], spt[7:44])
    r1 = re.findall(r"\d+", geoProj1[7:44])
    r2 = re.findall(r"\d+", spt[7:44])
    # print(r1,r2)
    if (r1 == r2):
        print(True)
    else:
        print(False, "\nERROR！数据坐标系不一致，请重新输入")
        sys.exit()

    # 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 = srcImage.ReadAsArray(ulX, ulY, pxWidth, pxHeight)  # ***只读要的那块***

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

    write_img(save_path, geoProj, geoTrans, clip)
    gdal.ErrorReset()


if __name__ == "__main__":
    shp = r"D:\Desktop\test\POLY.shp"

    img = r"F:\1115.img"
    out = r"D:\Desktop\test\out.tif"

    shpClipRaster(shp, img, out)