Python 使用 GDAL 和 geopandas 从城市建筑物高度数据集提取分市栅格数据

使用 Python 从城市建筑物高度数据集提取分市栅格数据

现有一套城市建筑物高度数据集，其分瓦片存储，数据格式为 tif 格式，全国范围的数据划分为 310 个瓦片，每个瓦片的大小为 2°×2°，其命名规则为 CNBH10m_X<中心经度>_Y<中心纬度>.tif，以 CNBH10m_X101_Y21.tif 这个文件名为例，表示该瓦片覆盖的范围是东经 100°——102°，北纬 20°——22°。

现在还有一个全国的分市 shp，要求遍历所有市，获取覆盖市的建筑物数据瓦片，并将其进行拼接，然后用 shp 裁剪拼接完的 tif，得到分市的建筑物高度栅格内容。

实现这个任务需要使用一些工具和库来进行地理空间数据处理。可以使用 Python 编程语言，并使用 GDAL（Geospatial Data Abstraction Library）和 geopandas 库来处理矢量和栅格数据。以下是一个大致的步骤：

1. 导入所需的库：

import os
import glob
import geopandas as gpd
from osgeo import gdal
from shapely.geometry import box

2. 读取分市的 shp 文件：

shp_path = 'path_to_shapefile.shp'
cities = gpd.read_file(shp_path)

3. 遍历每个市的边界，获取对应的建筑物数据瓦片：

tile_folder = 'path_to_tile_folder'
output_folder = 'path_to_output_folder'

for index, city in cities.iterrows():
    city_name = city['city_name']  # 假设 shp 中有一个名为'city_name'的列来存储市的名称
    
    # 获取市的边界框
    bbox = city.geometry.bounds
    xmin, ymin, xmax, ymax = bbox
    
    # 根据边界框计算对应的瓦片
    start_x = int(xmin // 2) * 2
    start_y = int(ymin // 2) * 2
    end_x = int(xmax // 2) * 2 + 2
    end_y = int(ymax // 2) * 2 + 2
    
    # 遍历瓦片并拼接
    tiles = []
    for x in range(start_x, end_x, 2):
        for y in range(start_y, end_y, 2):
            tile_name = f'CNBH10m_X{x}_Y{y}.tif'
            tile_path = os.path.join(tile_folder, tile_name)
            tiles.append(tile_path)
    
    # 使用 GDAL 库来进行瓦片拼接
    mosaic_path = os.path.join(output_folder, f'{city_name}_mosaic.tif')
    gdal.BuildVRT(mosaic_path, tiles)
    gdal.Translate(mosaic_path, mosaic_path, format='GTiff')

4. 使用 shp 裁剪拼接完的 tif：

for index, city in cities.iterrows():
    city_name = city['city_name']
    mosaic_path = os.path.join(output_folder, f'{city_name}_mosaic.tif')
    output_path = os.path.join(output_folder, f'{city_name}_height.tif')
    
    # 使用 GDAL 库进行裁剪
    gdal.Warp(output_path, mosaic_path, cutlineDSName=shp_path, cropToCutline=True)

这样，就可以遍历所有市，获取覆盖市的建筑物数据瓦片，并将其进行拼接，然后用 shp 裁剪拼接完的 tif，得到分市的建筑物高度栅格数据。请注意，代码中的路径需要根据实际情况进行修改。