用 Python 批量提取 1990–2019 年高层建筑并按城市导出 Shapefile
好的，下面是删除“创建虚拟环境”步骤后的完整可复制博客内容，结构已严格按你要求整理。
用 Python 批量提取 1990–2019 年高层建筑并按城市导出 Shapefile
背景介绍
如果你手里有一批多波段建筑高度栅格数据（TIF），每个波段代表一年（比如 B1=2019，B30=1990），同时还有各城市的行政区 Shapefile，你可能会遇到这样的问题：
👉 如何快速提取“高度大于 21 米”的建筑，并按“年份 + 城市”分别导出成矢量文件？
手动操作不仅耗时，而且非常容易出错。
本文这段代码可以帮你：
- 自动遍历 1990–2019 年数据
- 提取 高度 > 21m 的建筑
- 将栅格转为矢量多边形
- 按城市范围裁剪
- 输出为 每个城市每年的 Shapefile
- 自动生成缺失数据日志
运行完成后，你会看到类似这样的目录结构：
矢量化/├── 2019/│ ├── 郑州_2019.shp│ ├── 洛阳_2019.shp├── 2018/├── ...└── Missing_Data_Log.txt
代码功能说明
一句话总结这段代码：
它是一个“自动筛选 + 自动裁剪 + 自动导出”的批量处理工具。
整个流程可以理解为三步：
- 把所有栅格拼成完整地图
- 只保留高度大于 21 米的区域
- 按城市边界切开并分别保存
适合场景：
- 城市高层建筑时序分析
- 城市扩张研究
- 建筑高度分布研究
- 批量栅格矢量化处理
运行环境准备
1️⃣ Python 版本要求
建议使用 Python 3.8 及以上
查看版本：
python --version
如果低于 3.8，请升级 Python。
2️⃣ 安装依赖库
在命令行（cmd / PowerShell）输入：
pip install numpy pandas rasterio geopandas
如果安装 geopandas 报错，可以使用国内源：
pip install geopandas -i https://pypi.tuna.tsinghua.edu.cn/simple
为什么要安装这些？
- numpy：处理数组
- rasterio：读取栅格
- geopandas：处理矢量
- pandas：表格处理
详细运行步骤
下面按「从零到成功运行」一步步讲清楚。
第一步：准备数据文件夹
你的文件结构应该是：
CMAB_origin/├── 90_19/ # 存放所有 TIF 文件├── 城区/ # 存放城市 shp 文件
修改代码中的路径为你自己的路径：
raster_folder = r"D:\【your_path】\90_19"
shp_folder = r"D:\【your_path】\城区"
output_root = r"D:\【your_path】\矢量化"
注意：
- 路径前面必须加
r - 文件夹必须真实存在
- 不要写错盘符
第二步：保存代码文件
创建一个文件：
extract_building.py
把下面代码完整复制进去（已模糊私人路径）：
import os
import shutil
import numpy as np
import rasterio
from rasterio.merge import merge
from rasterio.features import shapes
import geopandas as gpd
import pandas as pd# ================= 1. 路径配置 =================
raster_folder = r"D:\【your_path】\90_19"
shp_folder = r"D:\【your_path】\城区"
output_root = r"D:\【your_path】\矢量化"log_file_path = os.path.join(output_root, "Missing_Data_Log.txt")# ================= 2. 环境清理 =================
if os.path.exists(output_root):shutil.rmtree(output_root)
os.makedirs(output_root)def get_file_size_mb(file_path):base = os.path.splitext(file_path)[0]total = 0for ext in ['.shp', '.shx', '.dbf', '.prj']:f = base + extif os.path.exists(f):total += os.path.getsize(f)return round(total / (1024 * 1024), 2)print("--- 开始处理 (波段倒序: B1=2019, B30=1990) ---")city_shps = [f for f in os.listdir(shp_folder) if f.endswith('.shp')]
tif_paths = [os.path.join(raster_folder, f) for f in os.listdir(raster_folder) if f.endswith('.tif')]city_gdf_list = []
base_crs = None
for c_shp in city_shps:path = os.path.join(shp_folder, c_shp)gdf = gpd.read_file(path)if base_crs is None: base_crs = gdf.crsif 'name' not in gdf.columns: gdf['name'] = "未知区"city_gdf_list.append((c_shp.split('_')[0], gdf.to_crs(base_crs)))missing_data_records = []for band_idx in range(1, 31):current_year = 2019 - (band_idx - 1)print(f"\n[波段 {band_idx} -> 年份 {current_year}] 正在合并并处理...")src_files = [rasterio.open(p) for p in tif_paths]mosaic_data, mosaic_transform = merge(src_files, indexes=[band_idx])out_crs = src_files[0].crsfor src in src_files: src.close()band_array = mosaic_data[0].astype('float32')mask = band_array > 21if not np.any(mask):for city_name, _ in city_gdf_list:missing_data_records.append(f"{current_year}年 | {city_name} | 全区域无符合高度数据")continueresults = ({'properties': {'height': v}, 'geometry': s}for s, v in shapes(band_array, mask=mask, transform=mosaic_transform))gdf_all_builds = gpd.GeoDataFrame.from_features(list(results), crs=out_crs)if gdf_all_builds.crs != base_crs:gdf_all_builds = gdf_all_builds.to_crs(base_crs)year_dir = os.path.join(output_root, str(current_year))if not os.path.exists(year_dir): os.makedirs(year_dir)for city_name, city_gdf in city_gdf_list:res_intersection = gpd.overlay(gdf_all_builds, city_gdf, how='intersection')if res_intersection.empty:missing_data_records.append(f"{current_year}年 | {city_name} | 城市范围内无符合高度数据")continueres_intersection['Year'] = current_yearres_intersection['City'] = city_nameres_intersection['District'] = res_intersection['name']out_path = os.path.join(year_dir, f"{city_name}_{current_year}.shp")final_cols = ['height', 'Year', 'City', 'District', 'geometry']res_intersection[final_cols].to_file(out_path, encoding='utf-8')print(f" -> 已导出: {city_name} | 大小: {get_file_size_mb(out_path)} MB")with open(log_file_path, 'w', encoding='utf-8') as f:f.write("=== 建筑高度矢量化缺失数据记录 ===\n")for line in missing_data_records:f.write(line + "\n")print("\n--- 任务全部完成 ---")
第三步：运行代码
进入代码所在目录：
cd D:\【your_script_folder】
执行：
python extract_building.py
第四步：检查输出结果
查看：
D:\【your_path】\矢量化
确认是否生成：
- 年份文件夹
- 每个城市的 shp 文件
- Missing_Data_Log.txt
核心代码解析（新手版）
① 栅格合并
mosaic_data, mosaic_transform = merge(...)
就像把很多地图小块拼成一整张完整大地图。
② 高度筛选
mask = band_array > 21
意思是：
“只保留高度大于 21 米的地方。”
③ 栅格转矢量
shapes(...)
就像给高层建筑区域描边，变成多边形。
④ 城市裁剪
gpd.overlay(...)
用城市边界把建筑图层切开。
常见问题解决
❌ No module named rasterio
解决：
pip install rasterio
❌ CRS mismatch 报错
原因：坐标系不一致
解决：确保 shp 和 tif 投影一致。
❌ 输出文件为空
可能原因：
- 没有高度 > 21m
- 波段顺序不对
测试方法：
mask = band_array > 10
❌ Windows 路径报错
错误写法：
"E:\data\test"
正确写法：
r"E:\data\test"
总结
这段代码实现了：
✔ 多年数据自动遍历
✔ 高层建筑自动提取
✔ 栅格转矢量
✔ 按城市裁剪
✔ 自动导出 + 日志记录