【Python】栅格处理：读取tif文件，将边界转为shapely.geometry的Polygon

原文作者：我辈李想
版权声明：文章原创，转载时请务必加上原文超链接、作者信息和本声明。

---

一、常见用法

from osgeo import gdal, ogr, osr
import shapely.wkb as wkblib
import shapely
from shapely import wkt
from shapely.geometry import Polygon

def get_image_extent(image_path):
    """
    获取影像文件的范围（边界）
    :param image_path: str, 影像文件路径
    :return: list, 范围（xmin, ymin, xmax, ymax）
    """
    ds = gdal.Open(image_path)
    if ds is None:
        print(f"Failed to open {image_path}")
        return None

    # 获取文件的边界
    print('获取文件的边界',ds.GetGeoTransform())
    xmin, xpixel, _, ymax, _, ypixel = ds.GetGeoTransform()
    xmax = xmin + (ds.RasterXSize * xpixel)
    ymin = ymax + (ds.RasterYSize * ypixel)

    srs = osr.SpatialReference(wkt=ds.GetProjection())
    target_srs = osr.SpatialReference()
    target_srs.ImportFromEPSG(4326)
    
    return [xmin, ymin, xmax, ymax],srs,target_srs


def create_poly_by_extent(extent, srs=None, target_srs=None):
    """
    根据范围（边界）创建多边形几何对象
    :param extent: list, 范围（xmin, ymin, xmax, ymax）
    :param srs: osr.SpatialReference, 文件的投影和地理坐标系信息，可选
    :param target_srs: osr.SpatialReference, 目标投影和地理坐标系信息，可选
    :return: ogr.Geometry, 多边形几何对象
    """
    if len(extent) != 4:
        print("Invalid extent")
        return None

    xmin, ymin, xmax, ymax = extent
    poly = ogr.Geometry(ogr.wkbPolygon)
    ring = ogr.Geometry(ogr.wkbLinearRing)
    ring.AddPoint(xmin, ymax)
    ring.AddPoint(xmin, ymin)
    ring.AddPoint(xmax, ymin)
    ring.AddPoint(xmax, ymax)
    ring.AddPoint(xmin, ymax)
    poly.AddGeometry(ring)

    print('poly',poly)
    if srs is not None and target_srs is not None:
        transform = osr.CoordinateTransformation(srs, target_srs)
        poly.Transform(transform)

    return poly

调用和打印

tif = os.path.join(os.path.abspath('.'),'Landcover','50R_20220101-20230101.tif')

extent, srs, target_srs = get_image_extent(tif)
poly = create_poly_by_extent(extent, srs, target_srs)
print('extent',extent)
print(type(poly),poly)
print(poly.ExportToWkt())
print(wkt.loads(poly.ExportToIsoWkt()))
shapely_geom = shapely.wkt.loads(poly.ExportToIsoWkt())
def _drop_m(x, y, z):
    return y,x
print(111,shapely.ops.transform(_drop_m, shapely_geom))

输出如下：

获取文件的边界 (194776.85543130862, 10.0, 0.0, 3544361.5057707, 0.0, -10.0)
poly POLYGON ((194776.855431309 3544361.5057707 0,194776.855431309 2654231.5057707 0,805226.855431309 2654231.5057707 0,805226.855431309 3544361.5057707 0,194776.855431309 3544361.5057707 0))
extent [194776.85543130862, 2654231.5057707, 805226.8554313086, 3544361.5057707]
<class 'osgeo.ogr.Geometry'> POLYGON ((31.9942608734697 113.769556027249 0,23.9707114161094 114.000718766114 0,23.9707107036005 119.999317659786 0,31.9942598738755 120.230483194614 0,31.9942608734697 113.769556027249 0))
POLYGON ((31.9942608734697 113.769556027249 0,23.9707114161094 114.000718766114 0,23.9707107036005 119.999317659786 0,31.9942598738755 120.230483194614 0,31.9942608734697 113.769556027249 0))
POLYGON Z ((31.9942608734697 113.769556027249 0, 23.9707114161094 114.000718766114 0, 23.9707107036005 119.999317659786 0, 31.9942598738755 120.230483194614 0, 31.9942608734697 113.769556027249 0))
111 POLYGON ((113.769556027249 31.9942608734697, 114.000718766114 23.9707114161094, 119.999317659786 23.9707107036005, 120.230483194614 31.9942598738755, 113.769556027249 31.9942608734697))

二、

	def create_poly_by_extent(extent, srs=None, target_srs=None):
	    """
	    根据范围（边界）创建多边形几何对象
	    :param extent: list, 范围（xmin, ymin, xmax, ymax）
	    :param srs: osr.SpatialReference, 文件的投影和地理坐标系信息，可选
	    :param target_srs: osr.SpatialReference, 目标投影和地理坐标系信息，可选
	    :return: ogr.Geometry, 多边形几何对象
	    """
	    if len(extent) != 4:
	        print("Invalid extent")
	        return None
	
	    xmin, ymin, xmax, ymax = extent
	    image_extent  = 'POLYGON ((%f %f,%f %f,%f %f,%f %f,%f %f))' %(
	                    min_x, min_y, min_x, max_y, max_x, max_y, max_x, min_y, min_x, min_y)
	    poly = ogr.CreateGeometryFromWkt(image_extent)  # 定义tif范围

	    return poly