GDAL库的安装、矢量和栅格数据的加载、数据文件信息输出、文件坐标系转换

1.GDAL库

GDAL（Geospatial Data Abstraction Library）是一个开源的栅格与矢量地理空间数据转换库，在功能上类似于arcgis系列的arcpy包。GDAL是C++编写的库，但它提供了非常完善的Python绑定，可以使用python调用。

1.1GDAL库安装

Windows系统下GDAL库的安装：建议使用conda install -c conda-forge gdal进行安装。不使用pip是因为，GDAL本质上是一个C++库，而不是python库，直接用pip下载会让PyPi下载C++的源码包，而不是.py的文件，需要额外的多余操作。使用-c conda-forge在conda-forge社区频道下载的包版本更新，更稳定。

1.2GDAL的安装检验

from osgeo import gdal, ogr, osr

print("GDAL 安装验证")
print(f"GDAL 版本: {gdal.__version__}")

# 获取 GDAL 数据目录
data_dir = gdal.GetConfigOption('GDAL_DATA')
print(f"GDAL 数据目录: {data_dir}")

# 获取驱动数量
driver_count = gdal.GetDriverCount()
print(f"GDAL 驱动数量: {driver_count}")

# OGR 驱动数量
ogr_driver_count = ogr.GetDriverCount()
print(f"OGR 驱动数量: {ogr_driver_count}")

print("模块 (gdal, ogr, osr) 导入成功")

2栅格数据

2.1读取栅格数据文件

from osgeo import gdal  # 主要处理栅格数据

gdal.UseExceptions()  # 开启异常模式（gdal4之后默认）：出错时抛出异常
# gdal.DontUseExceptions() 不开启异常模式

if __name__=="__main__":
 
     # 读取影像文件
    raster_path = r""
    ds = gdal.Open(raster_path)

    # 关闭影像文件
    # GDAL 的 Python 绑定是通过 SWIG 生成的，其资源释放机制依赖于引用计数，而不是 Python 的 del
    ds = None  

2.2获取栅格数据文件基本信息

from osgeo import gdal  # 主要处理栅格数据

# ===================================
gdal.UseExceptions()  # 开启异常模式（gdal4之后默认）：出错时抛出异常
# gdal.DontUseExceptions() 不开启异常模式

# ---------- 检查栅格----------
def check_raster(ds):
    print("\n【栅格数据检查】")
    if ds is None:
        print("ds为None")
    else:
        print(f"文件路径:{ds.GetDescription()}")  # GetDescription() 返回文件路径
        print(f"栅格尺寸: {ds.RasterXSize} 列 × {ds.RasterYSize} 行")
        print(f"波段数量: {ds.RasterCount}")

        # 获取投影信息
        proj = ds.GetProjection()
        print(f"投影信息: {proj[:200]}..." if len(proj) > 200 else f"投影信息: {proj}")

        # 获取地理变换参数
        gt = ds.GetGeoTransform()
        if gt:
            print(f"左上角X坐标: {gt[0]:.6f}")  # 左上角像元的X坐标
            print(f"像元宽度: {gt[1]:.8f}")  # 像元宽度（X 方向分辨率）
            print(f"行旋转参数：{gt[2]:6f}")  # 行旋转参数（X 方向随行的变化）
            print(f"左上角Y坐标: {gt[3]:.6f}")  # 左上角像元的 Y 坐标
            print(f"列旋转参数：{gt[4]:6f}")  # 列旋转参数（Y 方向随列的变化）
            print(f"像元高度: {gt[5]:.8f}")  # 像元高度（Y 方向分辨率），通常为负值表示Y轴向下

        # 计算并打印坐标范围
        xmin = gt[0]
        xmax = gt[0] + ds.RasterXSize * gt[1]
        ymax = gt[3]  # 左上角Y是最大值
        ymin = gt[3] + ds.RasterYSize * gt[5]  # gt[5]是负数
        print(f"坐标范围: X[{xmin:.6f}, {xmax:.6f}], Y[{ymin:.6f}, {ymax:.6f}]")


if __name__ == "__main__":
    # 栅格数据文件读取
    raster_path = r""
    ds = gdal.Open(raster_path)

    # 基础信息读取
    check_raster(ds)

    # GDAL 的 Python 绑定是通过 SWIG 生成的，其资源释放机制依赖于引用计数，而不是 Python 的 del
    ds = None  # 关闭文件

2.3转换栅格数据文件坐标系（重投影）

if __name__ == "__main__":
    # 栅格数据文件读取
    raster_path = r""
    output_raster = r""

    target_srs = "EPSG:4326"
    gdal.Warp(
        output_raster,      # 输出文件路径
        raster_path,        # 输入文件路径
        dstSRS=target_srs,  # 目标坐标系
        format="GTiff",     # 输出格式
        resampleAlg=gdal.GRA_Bilinear  # 重采样算法
    )

3.矢量数据

3.1读取矢量数据文件

from osgeo import gdal  # 主要处理栅格数据
from osgeo import ogr   # 主要处理矢量数据
gdal.UseExceptions()    #开启异常模式


if __name__=="__main__":
    # ========== 路径 ==========
    vector_path = r""
    output_vector = r""

    ds = ogr.Open(vector_path)
    check_voctor(ds)
    ds = None  # 关闭文件

3.2输出主要信息

from osgeo import gdal  # 主要处理栅格数据
from osgeo import ogr   # 主要处理矢量数据
gdal.UseExceptions()    #开启异常模式

def check_vector(ds):
    print("\n【矢量数据检查】")
    if ds is None:
        print("ds为None")
    else:
        layer = ds.GetLayer(0)
        print(f"文件路径: {ds.GetDescription()}")
        print(f"图层名称: {layer.GetName()}")
        print(f"要素数量: {layer.GetFeatureCount()} 个站点")

        # 获取投影信息
        spatial_ref = layer.GetSpatialRef()
        if spatial_ref:
            print(f"投影信息: {spatial_ref.GetName()}")
        else:
            print("投影信息: 无（可能是地理坐标系）")

        # 获取范围
        extent = layer.GetExtent()
        print(f"坐标范围: X[{extent[0]:.6f}, {extent[1]:.6f}], Y[{extent[2]:.6f}, {extent[3]:.6f}]")

        # 获取字段信息
        layer_defn = layer.GetLayerDefn()
        print("字段列表:")
        for i in range(layer_defn.GetFieldCount()):
            field_defn = layer_defn.GetFieldDefn(i)
            print(f"  [{i}] {field_defn.GetName()}")

if __name__=="__main__":
    # ========== 路径 ==========
    vector_path = r""
    output_vector = r""

    ds = ogr.Open(vector_path)
    check_vector(ds)
    ds = None  # 关闭文件

3.3坐标系转换

    target_srs = "EPSG:32649"
    gdal.VectorTranslate(
        output_vector,
        vector_path,
        format="ESRI Shapefile",
        layerCreationOptions=["ENCODING=UTF-8"],  # 避免中文乱码
        dstSRS=target_srs  # 目标坐标系
    )