GDAL whl文件安装详细教程

一、准备工作

1. 确认Python和系统信息

# 查看Python版本（非常重要！）
python --version

# 查看详细版本信息
python -c "import sys; print(f'Python {sys.version_info.major}.{sys.version_info.minor}.{sys.version_info.micro}')"
print(f'系统架构: {platform.architecture()[0]}')

2. GDAL版本对应关系

Python版本	推荐GDAL版本	whl文件名示例
Python 3.7	GDAL 3.4.3	GDAL-3.4.3-cp37-cp37m-win_amd64.whl
Python 3.8	GDAL 3.4.3	GDAL-3.4.3-cp38-cp38-win_amd64.whl
Python 3.9	GDAL 3.4.3	GDAL-3.4.3-cp39-cp39-win_amd64.whl
Python 3.10	GDAL 3.6.4	GDAL-3.6.4-cp310-cp310-win_amd64.whl
Python 3.11	GDAL 3.6.4	GDAL-3.6.4-cp311-cp311-win_amd64.whl
Python 3.12	GDAL 3.8.0	GDAL-3.8.0-cp312-cp312-win_amd64.whl

二、下载whl文件

目前直接用pip install gdal安装不上，只能通过whl安装。首先为了顺利安装这个模块，我们需要下载whl文件，可以尝试去github仓库pythonlibs_whl_mirror，然后找到对应版本whl文件，之后进行一个步骤安装

文件命名规则理解：

GDAL-{gdal版本}-cp{python主版本}{python次版本}-cp{python主版本}{python次版本}-win_amd64.whl
示例：GDAL-3.4.3-cp39-cp39-win_amd64.whl
表示：GDAL 3.4.3 for Python 3.9 64位

三、详细安装步骤

步骤1：安装Microsoft Visual C++ Redistributable（必需）

GDAL需要VC++运行库，下载并安装：

https://aka.ms/vs/17/release/vc_redist.x64.exe

步骤2：按顺序安装whl文件

# 打开CMD（管理员权限推荐）
# 切换到下载目录
cd C:\Users\你的用户名\Downloads

# 方法A：直接安装（如果只有一个whl文件）
pip install GDAL-3.6.4-cp311-cp311-win_amd64.whl

# 方法B：完整路径安装
pip install "C:\Users\你的用户名\Downloads\GDAL-3.6.4-cp311-cp311-win_amd64.whl"

# 方法C：如果有依赖问题，强制安装
pip install GDAL-3.6.4-cp311-cp311-win_amd64.whl --no-deps --force-reinstall

# 方法D：用户级安装（无管理员权限）
pip install --user GDAL-3.6.4-cp311-cp311-win_amd64.whl

步骤3：验证基本安装

# 验证安装
python -c "import osgeo; print('GDAL导入成功')"

# 查看版本
python -c "from osgeo import gdal; print(f'GDAL版本: {gdal.__version__}')"

四、验证安装的测试代码

测试1：基础功能测试

# test_gdal_basic.py
"""
GDAL 基础功能测试
"""

import sys
import os

print("=" * 60)
print("GDAL 安装测试")
print("=" * 60)

# 测试1：检查Python版本
print(f"Python版本: {sys.version}")
print(f"Python执行文件: {sys.executable}")
print(f"当前工作目录: {os.getcwd()}")
print()

# 测试2：尝试导入GDAL
try:
    from osgeo import gdal, ogr, osr, gdal_array
    print("✅ 成功导入GDAL核心模块")
    print(f"  - gdal: {gdal.__version__}")
    
    # 测试版本信息
    print(f"  - GDAL版本: {gdal.VersionInfo('RELEASE_NAME')}")
    print(f"  - 构建信息: {gdal.VersionInfo('BUILD_INFO')}")
    
except ImportError as e:
    print(f"❌ 导入GDAL失败: {e}")
    print("\n可能的原因：")
    print("1. GDAL未正确安装")
    print("2. Python版本不匹配")
    print("3. 缺少VC++运行库")
    print("4. 系统路径问题")
    sys.exit(1)

# 测试3：检查驱动支持
print("\n✅ 检查GDAL驱动支持：")
try:
    drivers = []
    for i in range(gdal.GetDriverCount()):
        driver = gdal.GetDriver(i)
        drivers.append(driver.ShortName)
    
    print(f"  共发现 {len(drivers)} 个数据驱动")
    print(f"  前10个驱动: {', '.join(drivers[:10])}")
    
    # 检查常用驱动
    essential_drivers = ['GTiff', 'HFA', 'JPEG', 'PNG', 'VRT']
    available = []
    missing = []
    
    for driver in essential_drivers:
        if driver in drivers:
            available.append(driver)
        else:
            missing.append(driver)
    
    print(f"  ✓ 可用驱动: {', '.join(available)}")
    if missing:
        print(f"  ⚠ 缺少驱动: {', '.join(missing)}")
    
except Exception as e:
    print(f"  ⚠ 检查驱动时出错: {e}")

# 测试4：检查OGR（矢量）驱动
print("\n✅ 检查OGR驱动支持：")
try:
    from osgeo import ogr
    
    ogr_drivers = []
    for i in range(ogr.GetDriverCount()):
        driver = ogr.GetDriver(i)
        ogr_drivers.append(driver.GetName())
    
    print(f"  共发现 {len(ogr_drivers)} 个矢量驱动")
    print(f"  前10个驱动: {', '.join(ogr_drivers[:10])}")
    
except Exception as e:
    print(f"  ⚠ 检查OGR驱动时出错: {e}")

# 测试5：检查PROJ（坐标系统）
print("\n✅ 检查PROJ支持：")
try:
    from osgeo import osr
    
    # 创建WGS84坐标系
    srs = osr.SpatialReference()
    srs.ImportFromEPSG(4326)  # WGS84
    print(f"  PROJ版本: {osr.SRS_VERSION}")
    print(f"  WGS84坐标系: {srs.GetName()}")
    
    # 测试坐标转换
    srs_utm = osr.SpatialReference()
    srs_utm.ImportFromEPSG(32650)  # UTM Zone 50N
    
    transform = osr.CoordinateTransformation(srs, srs_utm)
    print(f"  坐标转换功能正常")
    
except Exception as e:
    print(f"  ⚠ 检查PROJ时出错: {e}")

# 测试6：内存读写测试
print("\n✅ 内存读写测试：")
try:
    # 创建内存数据集
    driver = gdal.GetDriverByName('MEM')
    dataset = driver.Create('test_mem', 100, 100, 1, gdal.GDT_Byte)
    
    # 写入数据
    band = dataset.GetRasterBand(1)
    data = [[i % 256 for i in range(100)] for _ in range(100)]
    band.WriteArray(data)
    
    # 读取数据
    read_data = band.ReadAsArray()
    print(f"  内存数据集创建成功")
    print(f"  数据形状: {read_data.shape}")
    print(f"  数据类型: {read_data.dtype}")
    
    # 清理
    dataset = None
    
except Exception as e:
    print(f"  ⚠ 内存测试时出错: {e}")

# 测试7：文件格式支持测试
print("\n✅ 文件格式支持测试：")
try:
    test_formats = {
        'GeoTIFF': 'GTiff',
        'JPEG': 'JPEG',
        'PNG': 'PNG',
        'HFA/Erdas': 'HFA',
        'VRT': 'VRT'
    }
    
    for format_name, driver_name in test_formats.items():
        driver = gdal.GetDriverByName(driver_name)
        if driver:
            metadata = driver.GetMetadata()
            if 'DCAP_CREATE' in metadata:
                capability = "读写"
            elif 'DCAP_CREATECOPY' in metadata:
                capability = "复制"
            else:
                capability = "只读"
            print(f"  {format_name:15} ({driver_name:8}): ✓ 支持 [{capability}]")
        else:
            print(f"  {format_name:15} ({driver_name:8}): ✗ 不支持")

except Exception as e:
    print(f"  ⚠ 格式测试时出错: {e}")

print("\n" + "=" * 60)
print("测试完成！")
print("=" * 60)

if 'gdal' in locals():
    print(f"\n🎉 GDAL {gdal.__version__} 安装成功！")
    print("所有基本功能测试通过。")
else:
    print("\n❌ GDAL安装存在问题，请检查以上错误信息。")

测试2：栅格数据处理示例

# gdal_raster_demo.py
"""
GDAL 栅格数据处理示例
创建、读取、处理栅格数据
"""

import numpy as np
import matplotlib.pyplot as plt
from osgeo import gdal, osr
import tempfile
import os

def create_test_raster():
    """创建测试栅格数据"""
    print("1. 创建测试栅格数据...")
    
    # 创建临时文件
    temp_dir = tempfile.mkdtemp()
    output_file = os.path.join(temp_dir, "test_raster.tif")
    
    # 栅格参数
    cols = 500
    rows = 500
    bands = 3  # RGB三个波段
    data_type = gdal.GDT_Byte
    
    # 创建数据集
    driver = gdal.GetDriverByName('GTiff')
    dataset = driver.Create(output_file, cols, rows, bands, data_type)
    
    if dataset is None:
        print("创建数据集失败！")
        return None
    
    # 设置地理参考信息
    # 设置仿射变换参数
    # 左上角X坐标, X方向分辨率, 旋转参数, 左上角Y坐标, 旋转参数, Y方向分辨率
    geotransform = (100.0, 0.1, 0.0, 50.0, 0.0, -0.1)
    dataset.SetGeoTransform(geotransform)
    
    # 设置坐标系（WGS84）
    srs = osr.SpatialReference()
    srs.ImportFromEPSG(4326)  # WGS84
    dataset.SetProjection(srs.ExportToWkt())
    
    # 创建测试数据（渐变效果）
    print("2. 生成渐变数据...")
    for band_num in range(1, bands + 1):
        band = dataset.GetRasterBand(band_num)
        
        # 创建不同的渐变模式
        if band_num == 1:  # 红色波段 - 水平渐变
            data = np.fromfunction(lambda i, j: (j / cols * 255).astype(np.uint8), (rows, cols))
        elif band_num == 2:  # 绿色波段 - 垂直渐变
            data = np.fromfunction(lambda i, j: (i / rows * 255).astype(np.uint8), (rows, cols))
        else:  # 蓝色波段 - 对角渐变
            data = np.fromfunction(lambda i, j: ((i + j) / (rows + cols) * 255).astype(np.uint8), (rows, cols))
        
        band.WriteArray(data)
        band.SetNoDataValue(0)
        
        # 设置统计信息
        band.ComputeStatistics(False)
        
        # 设置颜色解释
        if band_num == 1:
            band.SetColorInterpretation(gdal.GCI_RedBand)
        elif band_num == 2:
            band.SetColorInterpretation(gdal.GCI_GreenBand)
        elif band_num == 3:
            band.SetColorInterpretation(gdal.GCI_BlueBand)
    
    # 构建金字塔（缩略图）
    print("3. 构建金字塔...")
    gdal.SetConfigOption('COMPRESS_OVERVIEW', 'DEFLATE')
    dataset.BuildOverviews("NEAREST", [2, 4, 8, 16])
    
    # 关闭数据集
    dataset = None
    
    print(f"4. 测试栅格已创建: {output_file}")
    return output_file

def read_and_analyze_raster(filepath):
    """读取并分析栅格数据"""
    print("\n读取和分析栅格数据...")
    
    # 打开数据集
    dataset = gdal.Open(filepath, gdal.GA_ReadOnly)
    if dataset is None:
        print(f"无法打开文件: {filepath}")
        return
    
    # 基本信息
    print(f"文件: {filepath}")
    print(f"尺寸: {dataset.RasterXSize} x {dataset.RasterYSize}")
    print(f"波段数: {dataset.RasterCount}")
    
    # 地理信息
    geotransform = dataset.GetGeoTransform()
    print(f"地理变换参数: {geotransform}")
    print(f"投影: {dataset.GetProjection()}")
    
    # 读取所有波段数据
    all_bands = []
    for i in range(1, dataset.RasterCount + 1):
        band = dataset.GetRasterBand(i)
        
        # 波段信息
        print(f"\n波段 {i}:")
        print(f"  数据类型: {gdal.GetDataTypeName(band.DataType)}")
        print(f"  无数据值: {band.GetNoDataValue()}")
        
        # 统计信息
        min_val, max_val, mean_val, std_val = band.GetStatistics(True, True)
        print(f"  最小值: {min_val:.2f}, 最大值: {max_val:.2f}")
        print(f"  平均值: {mean_val:.2f}, 标准差: {std_val:.2f}")
        
        # 读取数据
        data = band.ReadAsArray()
        all_bands.append(data)
        
        print(f"  数据形状: {data.shape}")
    
    # 创建可视化
    print("\n创建可视化图表...")
    fig, axes = plt.subplots(2, 3, figsize=(15, 10))
    
    # 显示原始波段
    titles = ['Red Band (水平渐变)', 'Green Band (垂直渐变)', 'Blue Band (对角渐变)']
    for i, (data, title) in enumerate(zip(all_bands, titles)):
        ax = axes[0, i]
        im = ax.imshow(data, cmap='viridis', vmin=0, vmax=255)
        ax.set_title(title)
        ax.axis('off')
        plt.colorbar(im, ax=ax, shrink=0.8)
    
    # 计算NDVI（模拟，使用波段差异）
    if len(all_bands) >= 2:
        # 模拟NDVI计算（实际上需要近红外和红光波段）
        simulated_ndvi = (all_bands[1].astype(float) - all_bands[0].astype(float)) / \
                         (all_bands[1].astype(float) + all_bands[0].astype(float) + 0.001)
        
        ax = axes[1, 0]
        im = ax.imshow(simulated_ndvi, cmap='RdYlGn', vmin=-1, vmax=1)
        ax.set_title('模拟NDVI')
        ax.axis('off')
        plt.colorbar(im, ax=ax, shrink=0.8)
    
    # 显示RGB合成
    if len(all_bands) >= 3:
        rgb = np.stack([all_bands[0], all_bands[1], all_bands[2]], axis=-1)
        ax = axes[1, 1]
        ax.imshow(rgb)
        ax.set_title('RGB合成')
        ax.axis('off')
    
    # 显示直方图
    ax = axes[1, 2]
    colors = ['red', 'green', 'blue']
    for i, (data, color) in enumerate(zip(all_bands, colors)):
        ax.hist(data.flatten(), bins=50, alpha=0.5, color=color, 
                label=f'Band {i+1}', density=True)
    ax.set_title('波段直方图')
    ax.set_xlabel('像素值')
    ax.set_ylabel('频率')
    ax.legend()
    ax.grid(True, alpha=0.3)
    
    plt.suptitle('GDAL 栅格数据处理演示', fontsize=16, fontweight='bold')
    plt.tight_layout()
    
    # 保存图片
    output_img = "gdal_raster_analysis.png"
    plt.savefig(output_img, dpi=300, bbox_inches='tight')
    print(f"分析图表已保存: {output_img}")
    
    plt.show()
    
    # 关闭数据集
    dataset = None
    
    return all_bands

def reproject_raster(input_file):
    """重投影栅格数据"""
    print("\n执行栅格重投影...")
    
    # 创建输出文件
    temp_dir = tempfile.mkdtemp()
    output_file = os.path.join(temp_dir, "reprojected.tif")
    
    # 定义源和目标坐标系
    src_srs = osr.SpatialReference()
    src_srs.ImportFromEPSG(4326)  # WGS84
    
    dst_srs = osr.SpatialReference()
    dst_srs.ImportFromEPSG(3857)  # Web Mercator
    
    # 重投影选项
    warp_options = gdal.WarpOptions(
        format='GTiff',
        srcSRS=src_srs,
        dstSRS=dst_srs,
        resampleAlg=gdal.GRA_Bilinear,
        dstNodata=0
    )
    
    # 执行重投影
    gdal.Warp(output_file, input_file, options=warp_options)
    
    print(f"重投影完成: {output_file}")
    return output_file

if __name__ == "__main__":
    try:
        print("=" * 60)
        print("GDAL 栅格数据处理演示")
        print("=" * 60)
        
        # 步骤1：创建测试栅格
        raster_file = create_test_raster()
        
        if raster_file and os.path.exists(raster_file):
            # 步骤2：读取和分析
            data = read_and_analyze_raster(raster_file)
            
            # 步骤3：重投影演示
            reprojected_file = reproject_raster(raster_file)
            
            print("\n" + "=" * 60)
            print("✅ 栅格数据处理演示完成！")
            print("=" * 60)
            
            # 清理临时文件
            import shutil
            temp_dir1 = os.path.dirname(raster_file)
            temp_dir2 = os.path.dirname(reprojected_file)
            shutil.rmtree(temp_dir1, ignore_errors=True)
            shutil.rmtree(temp_dir2, ignore_errors=True)
            print("临时文件已清理")
        else:
            print("❌ 创建测试栅格失败")
            
    except Exception as e:
        print(f"❌ 发生错误: {e}")
        import traceback
        traceback.print_exc()

测试3：矢量数据处理示例

# gdal_vector_demo.py
"""
GDAL 矢量数据处理示例
创建、读取、处理矢量数据
"""

from osgeo import ogr, osr
import tempfile
import os
import random

def create_test_vector():
    """创建测试矢量数据"""
    print("1. 创建测试矢量数据...")
    
    # 创建临时文件
    temp_dir = tempfile.mkdtemp()
    output_file = os.path.join(temp_dir, "test_vector.gpkg")
    
    # 创建数据源（使用GeoPackage格式）
    driver = ogr.GetDriverByName('GPKG')
    if driver is None:
        print("GPKG驱动不可用！")
        return None
    
    # 创建数据源
    datasource = driver.CreateDataSource(output_file)
    
    # 创建坐标系（WGS84）
    srs = osr.SpatialReference()
    srs.ImportFromEPSG(4326)
    
    # 创建点图层
    print("2. 创建点图层...")
    point_layer = datasource.CreateLayer(
        "cities",
        srs,
        ogr.wkbPoint,
        ['OVERWRITE=YES']
    )
    
    # 添加属性字段
    point_layer.CreateField(ogr.FieldDefn("name", ogr.OFTString))
    point_layer.CreateField(ogr.FieldDefn("population", ogr.OFTInteger))
    point_layer.CreateField(ogr.FieldDefn("capital", ogr.OFTInteger))  # 0/1表示是否首都
    
    # 创建一些城市点
    cities = [
        ("北京", 116.4074, 39.9042, 21540000, 1),
        ("上海", 121.4737, 31.2304, 24870000, 0),
        ("广州", 113.2644, 23.1291, 18680000, 0),
        ("深圳", 114.0579, 22.5431, 17560000, 0),
        ("成都", 104.0668, 30.5728, 16580000, 0),
        ("武汉", 114.2986, 30.5844, 12330000, 0),
        ("西安", 108.9480, 34.2632, 12950000, 0),
        ("杭州", 120.1551, 30.2741, 11940000, 0),
        ("南京", 118.7969, 32.0603, 8500000, 0),
        ("重庆", 106.5516, 29.5630, 31240000, 0)
    ]
    
    for city_name, lon, lat, population, is_capital in cities:
        feature = ogr.Feature(point_layer.GetLayerDefn())
        feature.SetField("name", city_name)
        feature.SetField("population", population)
        feature.SetField("capital", is_capital)
        
        # 创建几何点
        point = ogr.Geometry(ogr.wkbPoint)
        point.AddPoint(lon, lat)
        feature.SetGeometry(point)
        
        point_layer.CreateFeature(feature)
        feature = None
    
    print(f"  添加了 {len(cities)} 个城市点")
    
    # 创建线图层（河流）
    print("3. 创建线图层...")
    line_layer = datasource.CreateLayer(
        "rivers",
        srs,
        ogr.wkbLineString,
        ['OVERWRITE=YES']
    )
    
    line_layer.CreateField(ogr.FieldDefn("name", ogr.OFTString))
    line_layer.CreateField(ogr.FieldDefn("length_km", ogr.OFTReal))
    
    # 创建一些河流线
    rivers = [
        ("长江", [
            (91.11, 33.45), (94.5, 32.2), (101.0, 30.0),
            (104.0, 30.5), (108.0, 31.0), (112.0, 30.5),
            (115.0, 30.0), (118.0, 32.0), (121.0, 31.5)
        ], 6300),
        ("黄河", [
            (95.8, 35.0), (100.2, 36.5), (104.5, 36.0),
            (108.0, 34.5), (110.0, 35.0), (112.5, 34.8),
            (114.5, 34.9), (116.5, 37.5), (118.5, 37.8)
        ], 5464)
    ]
    
    for river_name, coordinates, length in rivers:
        feature = ogr.Feature(line_layer.GetLayerDefn())
        feature.SetField("name", river_name)
        feature.SetField("length_km", length)
        
        # 创建几何线
        line = ogr.Geometry(ogr.wkbLineString)
        for lon, lat in coordinates:
            line.AddPoint(lon, lat)
        
        feature.SetGeometry(line)
        line_layer.CreateFeature(feature)
        feature = None
    
    print(f"  添加了 {len(rivers)} 条河流")
    
    # 创建面图层（省份）
    print("4. 创建面图层...")
    polygon_layer = datasource.CreateLayer(
        "provinces",
        srs,
        ogr.wkbPolygon,
        ['OVERWRITE=YES']
    )
    
    polygon_layer.CreateField(ogr.FieldDefn("name", ogr.OFTString))
    polygon_layer.CreateField(ogr.FieldDefn("area_km2", ogr.OFTReal))
    
    # 创建一些省份面（简化版）
    provinces = [
        ("四川省", [
            (101.0, 32.0), (101.0, 26.0), (108.0, 26.0),
            (108.0, 32.0), (101.0, 32.0)
        ], 485000),
        ("广东省", [
            (109.0, 25.0), (109.0, 20.0), (117.0, 20.0),
            (117.0, 25.0), (109.0, 25.0)
        ], 179800),
        ("新疆", [
            (73.0, 49.0), (73.0, 34.0), (96.0, 34.0),
            (96.0, 49.0), (73.0, 49.0)
        ], 1660000)
    ]
    
    for province_name, coordinates, area in provinces:
        feature = ogr.Feature(polygon_layer.GetLayerDefn())
        feature.SetField("name", province_name)
        feature.SetField("area_km2", area)
        
        # 创建几何面
        ring = ogr.Geometry(ogr.wkbLinearRing)
        for lon, lat in coordinates:
            ring.AddPoint(lon, lat)
        
        polygon = ogr.Geometry(ogr.wkbPolygon)
        polygon.AddGeometry(ring)
        feature.SetGeometry(polygon)
        
        polygon_layer.CreateFeature(feature)
        feature = None
    
    print(f"  添加了 {len(provinces)} 个省份")
    
    # 关闭数据源
    datasource = None
    
    print(f"5. 测试矢量数据已创建: {output_file}")
    return output_file

def analyze_vector_data(filepath):
    """分析矢量数据"""
    print("\n分析矢量数据...")
    
    # 打开数据源
    datasource = ogr.Open(filepath, 0)  # 0表示只读
    if datasource is None:
        print(f"无法打开文件: {filepath}")
        return
    
    print(f"文件: {filepath}")
    print(f"图层数量: {datasource.GetLayerCount()}")
    
    # 遍历所有图层
    for layer_idx in range(datasource.GetLayerCount()):
        layer = datasource.GetLayerByIndex(layer_idx)
        layer_name = layer.GetName()
        
        print(f"\n图层 {layer_idx + 1}: {layer_name}")
        print(f"  要素数量: {layer.GetFeatureCount()}")
        print(f"  几何类型: {layer.GetGeomType()} ({ogr.GeometryTypeToName(layer.GetGeomType())})")
        
        # 获取空间范围
        extent = layer.GetExtent()
        print(f"  空间范围: X({extent[0]:.2f}, {extent[1]:.2f}), Y({extent[2]:.2f}, {extent[3]:.2f})")
        
        # 获取属性字段
        layer_defn = layer.GetLayerDefn()
        print(f"  字段数量: {layer_defn.GetFieldCount()}")
        
        print("  字段信息:")
        for i in range(layer_defn.GetFieldCount()):
            field_defn = layer_defn.GetFieldDefn(i)
            print(f"    {i+1}. {field_defn.GetName():20} ({field_defn.GetTypeName():10})")
        
        # 读取前3个要素
        print(f"\n  前3个要素示例:")
        layer.ResetReading()
        for i in range(min(3, layer.GetFeatureCount())):
            feature = layer.GetNextFeature()
            if feature is not None:
                geom = feature.GetGeometryRef()
                print(f"    要素 {i+1}: ", end="")
                
                # 显示属性
                for j in range(feature.GetFieldCount()):
                    field_name = layer_defn.GetFieldDefn(j).GetName()
                    field_value = feature.GetField(j)
                    print(f"{field_name}: {field_value}, ", end="")
                
                # 显示几何信息
                if geom:
                    geom_type = geom.GetGeometryName()
                    if geom_type == 'POINT':
                        print(f"位置: ({geom.GetX():.4f}, {geom.GetY():.4f})")
                    elif geom_type == 'LINESTRING':
                        print(f"长度: {geom.Length():.1f}度")
                    elif geom_type == 'POLYGON':
                        print(f"面积: {geom.Area():.1f}平方度")
                print()
        
        # 空间分析示例
        if layer.GetGeomType() == ogr.wkbPoint:
            print(f"\n  空间分析:")
            
            # 计算中心点
            layer.ResetReading()
            points = []
            feature = layer.GetNextFeature()
            while feature:
                geom = feature.GetGeometryRef()
                if geom:
                    points.append((geom.GetX(), geom.GetY()))
                feature = layer.GetNextFeature()
            
            if points:
                center_x = sum(p[0] for p in points) / len(points)
                center_y = sum(p[1] for p in points) / len(points)
                print(f"    中心点: ({center_x:.4f}, {center_y:.4f})")
                print(f"    点数: {len(points)}")
    
    # 关闭数据源
    datasource = None

def query_vector_data(filepath):
    """查询矢量数据"""
    print("\n执行矢量数据查询...")
    
    datasource = ogr.Open(filepath, 0)
    if datasource is None:
        return
    
    # 查询点图层
    point_layer = datasource.GetLayerByName("cities")
    if point_layer:
        print("1. 查询人口超过2000万的城市:")
        point_layer.SetAttributeFilter("population > 20000000")
        for feature in point_layer:
            name = feature.GetField("name")
            population = feature.GetField("population")
            print(f"  - {name}: {population:,} 人")
        
        print("\n2. 查询省会城市:")
        point_layer.SetAttributeFilter("capital = 1")
        for feature in point_layer:
            name = feature.GetField("name")
            print(f"  - {name}")
    
    # 查询面图层
    polygon_layer = datasource.GetLayerByName("provinces")
    if polygon_layer:
        print("\n3. 查询面积超过50万km²的省份:")
        polygon_layer.SetAttributeFilter("area_km2 > 500000")
        for feature in polygon_layer:
            name = feature.GetField("name")
            area = feature.GetField("area_km2")
            print(f"  - {name}: {area:,.0f} km²")
    
    datasource = None

def export_to_shapefile(filepath):
    """导出为Shapefile格式"""
    print("\n导出为Shapefile格式...")
    
    temp_dir = tempfile.mkdtemp()
    shapefile_path = os.path.join(temp_dir, "cities_export.shp")
    
    # 打开源数据
    src_datasource = ogr.Open(filepath, 0)
    src_layer = src_datasource.GetLayerByName("cities")
    
    # 创建Shapefile驱动
    driver = ogr.GetDriverByName("ESRI Shapefile")
    
    # 创建目标数据源
    dst_datasource = driver.CreateDataSource(shapefile_path)
    
    # 复制图层
    dst_layer = dst_datasource.CopyLayer(src_layer, "cities")
    
    # 关闭数据源
    dst_datasource = None
    src_datasource = None
    
    # 检查创建的文件
    files = [f for f in os.listdir(temp_dir) if f.startswith("cities_export")]
    print(f"  导出的文件: {', '.join(files)}")
    print(f"  导出路径: {shapefile_path}")
    
    return shapefile_path

if __name__ == "__main__":
    try:
        print("=" * 60)
        print("GDAL 矢量数据处理演示")
        print("=" * 60)
        
        # 步骤1：创建测试矢量数据
        vector_file = create_test_vector()
        
        if vector_file and os.path.exists(vector_file):
            # 步骤2：分析矢量数据
            analyze_vector_data(vector_file)
            
            # 步骤3：查询数据
            query_vector_data(vector_file)
            
            # 步骤4：导出为其他格式
            exported_file = export_to_shapefile(vector_file)
            
            print("\n" + "=" * 60)
            print("✅ 矢量数据处理演示完成！")
            print("=" * 60)
            
            # 清理
            import shutil
            temp_dir1 = os.path.dirname(vector_file)
            temp_dir2 = os.path.dirname(exported_file)
            shutil.rmtree(temp_dir1, ignore_errors=True)
            shutil.rmtree(temp_dir2, ignore_errors=True)
            print("临时文件已清理")
        else:
            print("❌ 创建测试矢量数据失败")
            
    except Exception as e:
        print(f"❌ 发生错误: {e}")
        import traceback
        traceback.print_exc()

五、常见错误解决

错误1：ImportError: DLL load failed

ImportError: DLL load failed while importing _gdal: 找不到指定的模块。

解决：

1. 安装VC++ Redistributable
2. 检查Python版本匹配
3. 重启计算机

错误2：版本不匹配

ERROR: GDAL-3.6.4-cp311-cp311-win_amd64.whl is not a supported wheel on this platform

解决：

# 检查Python版本
python --version

# 下载对应版本的whl文件
# Python 3.11 → cp311
# Python 3.10 → cp310
# Python 3.9 → cp39

错误3：缺少依赖

ERROR: Could not install packages due to an OSError

解决：

# 使用--user选项
pip install --user GDAL-*.whl

# 或使用--no-deps跳过