现有一套全国的建筑物数据集其数据格式为tif格式全国范围的数据划分为310个瓦片每个瓦片的大小为2°×2°存储在一个文件夹内每个瓦片的大小为2°×2°对于每个瓦片文件命名方式为:CNBH10m_X中心经度_Y中心纬度tif每个中心经纬度都是奇数以CNBH10m_X101_Y21tif这个文件名为例表示该瓦片覆盖的范围是东经100°—102°北纬20°—22°现有一个全国分省市的shp遍历整个shp
要实现上述需求,可以按照以下步骤进行:
- 读取全国分省市的shp文件,获取每个市的边界范围(经纬度)和市名称。
- 遍历每个市的边界范围,确定其所涵盖的瓦片文件。
- 读取对应瓦片文件的建筑物数据。
- 将建筑物数据进行拼接。
- 将拼接后的数据保存为对应市名称的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进行安装。
原文地址: https://www.cveoy.top/t/topic/i5PB 著作权归作者所有。请勿转载和采集!