获取TIF栅格的数据范围边框（不是外接矩形）保存为shp

        在论文写作中，我们经常需要绘制研究区和研究数据分布图（如图1所示），那这样的数据分布框（图1 中蓝色框或图2中黄色框）怎么获取呐？作者在制作研究数据分布图时，通过搜索发现大部分的教程都是生成TIF的外接矩形框（图2中红色框）而不是数据范围框（图1 中蓝色框或图2中黄色框）。接下来看一看如何获取数据的范围边框。

 图1  研究区和研究数据分布图

图2 外接矩形框和数据范围框

方法

 图3 方法示意图

方法思路其实很简答，获取顶点1/2/3/4的坐标，生成面或者多边形写入shp就行。

代码

import pandas as pd
from osgeo import gdal, ogr
import os
import numpy as np
from tqdm import tqdm
def get_extend(tifname,shpname,no_data):
    # 打开TIFF文件
    dataset = gdal.Open(tifname)

    # 获取地理参考信息
    geotransform = dataset.GetGeoTransform()
    xmin = geotransform[0]
    ymax = geotransform[3]
    xmax = xmin + geotransform[1] * dataset.RasterXSize
    ymin = ymax + geotransform[5] * dataset.RasterYSize

    # 获取四点坐标 创建Polygon对象
    band = dataset.GetRasterBand(1)
    array = band.ReadAsArray()
    for i in range(0,array.shape[0]):
        a1=array[i,:]
        if np.sum(a1!=no_data)>0:
            index1=np.argmax(a1!=no_data)
            break
    for i in range(0, array.shape[0]):
        a2=array[-i,:]
        if np.sum(a2 != no_data) > 0:
            index2=np.argmax(a2!=no_data)
            break
    for i in range(0, array.shape[1]):
        a3=array[:,i]
        if np.sum(a3 != no_data) > 0:
            index3=np.argmax(a3!=no_data)
            break
    for i in range(0, array.shape[1]):
        a4=array[:,-i]
        if np.sum(a4 != no_data) > 0:
            index4=np.argmax(a4!=no_data)
            break

    ring = ogr.Geometry(ogr.wkbLinearRing)
    ring.AddPoint(xmin+index1*geotransform[1], ymax)
    ring.AddPoint(xmax, ymax+geotransform[5]*index4)
    ring.AddPoint(xmin+index2*geotransform[1], ymin)
    ring.AddPoint(xmin, ymax+geotransform[5]*index3)
    ring.AddPoint(xmin+index1*geotransform[1], ymax)
    polygon = ogr.Geometry(ogr.wkbPolygon)
    polygon.AddGeometry(ring)

    # 创建Shapefile文件
    shp_path = shpname
    driver = ogr.GetDriverByName('ESRI Shapefile')
    if os.path.exists(shp_path):
        driver.DeleteDataSource(shp_path)
    shp_ds = driver.CreateDataSource(shp_path)
    layer = shp_ds.CreateLayer('extent', geom_type=ogr.wkbPolygon)
    layer_defn = layer.GetLayerDefn()
    field_defn = ogr.FieldDefn('id', ogr.OFTInteger)
    layer.CreateField(field_defn)
    feature = ogr.Feature(layer_defn)
    feature.SetGeometry(polygon)
    feature.SetField('id', 1)
    layer.CreateFeature(feature)
    feature = None
    # 关闭文件
    shp_ds = None
    dataset = None

代码中主要变量与方法示意图中一致，大家可以相互对应理解，有帮助的话别忘了一键三连呀！