要实现上述需求,可以按照以下步骤进行:

  1. 读取全国分省市的shp文件,获取每个市的边界范围(经纬度)和市名称。
  2. 遍历每个市的边界范围,确定其所涵盖的瓦片文件。
  3. 读取对应瓦片文件的建筑物数据。
  4. 将建筑物数据进行拼接。
  5. 将拼接后的数据保存为对应市名称的tif文件。

下面是一个示例代码,可以通过gdal库来实现tif文件的读取和保存:

import os
import glob
from osgeo import gdal, gdal_array, ogr, osr

# 读取shp文件获取市边界范围和名称
def read_shp(shp_path):
    shp_ds = ogr.Open(shp_path)
    shp_lyr = shp_ds.GetLayer()
    
    cities = []
    for feature in shp_lyr:
        city_name = feature.GetField("市名称字段名")
        city_bbox = feature.geometry().GetEnvelope()
        cities.append((city_name, city_bbox))
    
    shp_ds = None
    return cities

# 获取指定瓦片文件的建筑物数据
def read_tile_data(tile_path):
    ds = gdal.Open(tile_path)
    data = ds.GetRasterBand(1).ReadAsArray()
    ds = None
    return data

# 将拼接后的数据保存为tif文件
def save_as_tif(data, extent, output_path):
    rows, cols = data.shape
    geotransform = (extent[0], (extent[2] - extent[0]) / cols, 0, extent[3], 0, (extent[1] - extent[3]) / rows)
    
    driver = gdal.GetDriverByName("GTiff")
    ds = driver.Create(output_path, cols, rows, 1, gdal.GDT_Byte)
    ds.SetGeoTransform(geotransform)
    
    srs = osr.SpatialReference()
    srs.ImportFromEPSG(4326)
    ds.SetProjection(srs.ExportToWkt())
    
    ds.GetRasterBand(1).WriteArray(data)
    ds = None

# 主函数
def process_tiles(shp_path, tiles_folder):
    cities = read_shp(shp_path)
    
    for city_name, city_bbox in cities:
        xmin, xmax, ymin, ymax = city_bbox
        
        # 确定所涵盖的瓦片文件
        tile_files = glob.glob(os.path.join(tiles_folder, f"CNBH10m_X*_Y*.tif"))
        tiles_to_process = []
        for tile_file in tile_files:
            tile_name = os.path.basename(tile_file)
            lon = int(tile_name.split("_")[1][1:])
            lat = int(tile_name.split("_")[2][1:-4])
            
            if lon >= xmin and lon <= xmax and lat >= ymin and lat <= ymax:
                tiles_to_process.append(tile_file)
        
        if len(tiles_to_process) == 0:
            continue
        
        # 读取并拼接瓦片数据
        merged_data = None
        for tile_file in tiles_to_process:
            tile_data = read_tile_data(tile_file)
            if merged_data is None:
                merged_data = tile_data
            else:
                merged_data = np.concatenate((merged_data, tile_data), axis=0)
        
        # 保存为tif文件
        output_path = f"{city_name}.tif"
        save_as_tif(merged_data, city_bbox, output_path)

# 示例使用
shp_path = "分省市shp文件路径"
tiles_folder = "瓦片文件夹路径"
process_tiles(shp_path, tiles_folder)

请注意,上述代码中的部分内容需要根据实际情况进行修改,比如替换相应的字段名、文件路径等。另外,代码中使用的是gdal库来处理tif文件,如果未安装该库,可以通过pip install gdal进行安装。

现有一套全国的建筑物数据集其数据格式为tif格式全国范围的数据划分为310个瓦片每个瓦片的大小为2°×2°存储在一个文件夹内每个瓦片的大小为2°×2°对于每个瓦片文件命名方式为:CNBH10m_X中心经度_Y中心纬度tif每个中心经纬度都是奇数以CNBH10m_X101_Y21tif这个文件名为例表示该瓦片覆盖的范围是东经100°—102°北纬20°—22°现有一个全国分省市的shp遍历整个shp

原文地址: https://www.cveoy.top/t/topic/i5PB 著作权归作者所有。请勿转载和采集!

免费AI点我,无需注册和登录