GIS 与遥感影像分析技能

环境配置

所有 Python 脚本必须使用隔离 venv 中的解释器：

Python: C:/Users/an/.workbuddy/binaries/python/envs/default/Scripts/python.exe
Pip:    C:/Users/an/.workbuddy/binaries/python/envs/default/Scripts/pip.exe

核心依赖库（已安装，2026-08-03 验证通过）：

• rasterio 1.5.0 — 栅格/遥感影像读写（GeoTIFF, ENVI, 等）
• geopandas 1.1.4 — 矢量数据处理（Shapefile, GeoJSON, 等）
• shapely 2.1.2 — 几何运算引擎
• fiona 1.10.1 — 矢量格式 I/O
• pyproj 3.7.2 — 坐标系转换
• folium 0.20.0 — 交互式 Web 地图
• matplotlib 3.11.1 — 静态制图
• earthpy 0.9.4 — 遥感可视化辅助
• xarray 2026.7.0 — 多维栅格数据
• rioxarray 0.23.0 — xarray 栅格扩展
• numpy 2.5.1 / scipy 1.18.0 — 数值计算
• pandas 2.3.3 — 数据分析（注意: 必须用 2.x，不能用 3.x，否则 xarray 不兼容）
• setuptools 80.10.2 — earthpy 依赖 pkg_resources（注意: 必须 <81，否则 pkg_resources 不可用）
• earthengine-api 1.7.37 — Google Earth Engine Python API（用于下载卫星影像）
  • 依赖: google-auth, google-api-python-client, google-cloud-storage, httplib2 等

GEE 影像下载（首次使用必读）

认证配置（只需一次）

通过 Google Earth Engine 下载影像前，需要先注册并认证：

1. 注册 GEE 账号：访问 https://code.earthengine.google.com ，用 Google 账号登录并完成注册（学术用途免费）
2. 命令行认证：运行以下命令，会打开浏览器进行 OAuth 授权

# 方式一：命令行（推荐）
"C:/Users/an/.workbuddy/binaries/python/envs/default/Scripts/earthengine.exe" authenticate

# 方式二：Python 代码
python -c "import ee; ee.Authenticate()"

3. 认证成功后，凭证会保存在 ~/.config/earthengine/credentials，之后无需重复认证

下载脚本

脚本路径：scripts/download_satellite_imagery.py

PYTHON="C:/Users/an/.workbuddy/binaries/python/envs/default/Scripts/python.exe"
SCRIPT="C:/Users/an/.workbuddy/skills/gis-remote-sensing/scripts/download_satellite_imagery.py"

# 下载太原市 Sentinel-2 影像（10m分辨率，2024年夏季，云量<20%）
$PYTHON $SCRIPT --place 太原市 --satellite sentinel2 \
    --start 2024-06-01 --end 2024-09-30 --cloud 20 --output taiyuan_s2.tif

# 下载 Landsat 8 影像（指定经纬度范围）
$PYTHON $SCRIPT --bbox 112.0 37.5 113.0 38.0 \
    --satellite landsat8 --start 2024-01-01 --end 2024-12-31 \
    --cloud 30 --output taiyuan_l8.tif

# 下载假彩色组合（NIR-Red-Green，适合植被分析）
$PYTHON $SCRIPT --place 太原市 --satellite landsat8 \
    --start 2024-06-01 --end 2024-09-30 --cloud 20 \
    --bands B5,B4,B3 --output taiyuan_nir.tif

# 用 shapefile 裁剪下载
$PYTHON $SCRIPT --shapefile boundary.shp \
    --satellite sentinel2 --start 2024-06-01 --end 2024-09-30 \
    --cloud 20 --output clipped.tif

# 查看可用影像列表（不下载）
$PYTHON $SCRIPT --place 太原市 --satellite sentinel2 \
    --start 2024-06-01 --end 2024-09-30 --cloud 20 --list-only

# 查看卫星波段信息
$PYTHON $SCRIPT --info

支持的卫星

| 卫星 | 数据集 | 分辨率 | 重访周期 | 常用场景 | |------|--------|--------|---------|---------| | Sentinel-2 | COPERNICUS/S2_SR_HARMONIZED | 10m | 5天 | 土地利用、植被监测、城市变化 | | Landsat 8 | LANDSAT/LC08/C02/T1_SR | 30m | 16天 | 长时间序列分析、热红外 | | Landsat 9 | LANDSAT/LC09/C02/T1_SR | 30m | 16天 | 与 Landsat 8 互补 | | MODIS NDVI | MODIS/061/MOD13A2 | 1km | 16天 | 大区域植被动态监测 |

常用波段组合

| 组合 | Landsat 8/9 | Sentinel-2 | 用途 | |------|-------------|------------|------| | 真彩色 RGB | B4,B3,B2 | B4,B3,B2 | 自然色彩，直观可视 | | 假彩色 NIR | B5,B4,B3 | B8,B4,B3 | 植被呈红色，适合植被分析 | | 短波红外 | B7,B5,B4 | B12,B8,B4 | 区分建筑/裸地/植被 | | 农业组合 | B6,B5,B2 | B11,B8,B2 | 农作物监测 |

下载策略

• 小区域（<8000×8000 像素）：单块下载，速度快
• 大区域（>8000×8000 像素）：自动分块下载并拼接，避免 GEE 单次请求限制
• 影像会自动做中值合成（median composite），减少云和阴影影响
• 反射率值会自动转换为物理意义值（0~1 范围）

核心工作流

1. 读取遥感影像（栅格数据）

import rasterio
import numpy as np

with rasterio.open("path/to/image.tif") as src:
    # 基本信息print(f"波段数: {src.count}")
    print(f"尺寸: {src.width} x {src.height}")
    print(f"坐标系: {src.crs}")
    print(f"分辨率: {src.res}")
    print(f"边界范围: {src.bounds}")

    # 读取所有波段
    data = src.read()  # shape: (bands, rows, cols)

    # 读取指定波段
    band1 = src.read(1)  # 第一波段

    # 获取 nodata 值
    nodata = src.nodata

2. 波段统计与直方图

import rasterio
import numpy as np
import matplotlib.pyplot as plt

with rasterio.open("image.tif") as src:
    for i in range(1, src.count + 1):
        band = src.read(i)
        valid = band[band != src.nodata] if src.nodata else band
        print(f"波段{i}: min={valid.min():.2f}, max={valid.max():.2f}, "
              f"mean={valid.mean():.2f}, std={valid.std():.2f}")

        # 直方图
        plt.figure(figsize=(8, 4))
        plt.hist(valid.flatten(), bins=100, color='steelblue', edgecolor='none')
        plt.title(f"Band {i} Histogram")
        plt.xlabel("Pixel Value")
        plt.ylabel("Frequency")
        plt.savefig(f"band{i}_hist.png", dpi=150, bbox_inches='tight')
        plt.close()

3. 遥感指数计算

NDVI（归一化植被指数）

# NIR = 近红外波段, RED = 红光波段
# NDVI = (NIR - RED) / (NIR + RED)
# 值域: -1 ~ 1, >0.5 密集植被, 0.2~0.5 稀疏植被, <0 水体/裸地

ndvi = (nir.astype(float) - red.astype(float)) / (nir + red + 1e-10)

NDWI（归一化水体指数）

# GREEN = 绿光波段, NIR = 近红外波段
# NDWI = (GREEN - NIR) / (GREEN + NIR)
# 值域: -1 ~ 1, >0 为水体

ndwi = (green.astype(float) - nir.astype(float)) / (green + nir + 1e-10)

NDBI（归一化建筑指数）

# SWIR = 短波红外, NIR = 近红外
# NDBI = (SWIR - NIR) / (SWIR + NIR)
# 值域: -1 ~ 1, >0 为建筑区

ndbi = (swir.astype(float) - nir.astype(float)) / (swir + nir + 1e-10)

保存计算结果为 GeoTIFF

import rasterio

with rasterio.open("input.tif") as src:
    profile = src.profile
    profile.update(dtype=rasterio.float32, count=1)

    with rasterio.open("ndvi.tif", 'w', **profile) as dst:
        dst.write(ndvi.astype(rasterio.float32), 1)

常用卫星波段对照 详见 references/satellite_bands.md

4. 矢量数据处理（GeoPandas）

import geopandas as gpd

# 读取 Shapefile / GeoJSON
gdf = gpd.read_file("boundary.shp")
# 或读取 GeoJSON
gdf = gpd.read_file("data.geojson")

# 查看信息
print(gdf.crs)          # 坐标系
print(gdf.columns)      # 属性字段
print(gdf.total_bounds) # 边界范围
print(gdf.head())       # 前几行

# 坐标系转换
gdf_4326 = gdf.to_crs(epsg=4326)   # WGS84
gdf_3857 = gdf.to_crs(epsg=3857)   # Web Mercator

# 空间筛选
point = gdf.geometry.iloc[0]
nearby = gdf[gdf.geometry.distance(point) < 0.01]

# 写出
gdf.to_file("output.shp", encoding='utf-8')
gdf.to_file("output.geojson", driver='GeoJSON')

5. 空间分析

缓冲区分析

# 注意: 缓冲区距离单位与 CRS 一致
# 度坐标系 (EPSG:4326) 中 0.01 ≈ 1km
# 投影坐标系 (如 EPSG:3857) 中单位为米
buffered = gdf.to_crs(epsg=3857).buffer(1000)  # 1km 缓冲区

空间叠加

# 交集
intersection = gpd.overlay(gdf1, gdf2, how='intersection')
# 并集
union = gpd.overlay(gdf1, gdf2, how='union')
# 差集
difference = gpd.overlay(gdf1, gdf2, how='difference')

空间连接（属性合并）

# 将点要素与面要素做空间连接，获取所在区域属性
joined = gpd.sjoin(points_gdf, polygons_gdf, how='left', predicate='within')

6. 栅格-矢量交互

用矢量裁剪栅格

import rasterio
from rasterio.mask import mask
import geopandas as gpd

shapefile = gpd.read_file("boundary.shp")
geoms = shapefile.geometry.values

with rasterio.open("image.tif") as src:
    # 确保 CRS 一致
    shapefile = shapefile.to_crs(src.crs)
    geoms = shapefile.geometry.values
    out_image, out_transform = mask(src, geoms, crop=True)
    out_meta = src.meta.copy()
    out_meta.update({
        "height": out_image.shape[1],
        "width": out_image.shape[2],
        "transform": out_transform
    })

with rasterio.open("clipped.tif", 'w', **out_meta) as dst:
    dst.write(out_image)

分区统计（Zonal Statistics）

# 需要 rasterstats 库: pip install rasterstats
from rasterstats import zonal_stats

stats = zonal_stats(
    "zones.shp",           # 矢量分区
    "raster.tif",           # 栅格数据
    stats=['mean', 'sum', 'min', 'max', 'count', 'std']
)
# 返回每个要素的统计字典列表

7. 坐标转换

from pyproj import Transformer, CRS

# WGS84 经纬度 -> 投影坐标
transformer = Transformer.from_crs("EPSG:4326", "EPSG:3857", always_xy=True)
x, y = transformer.transform(112.55, 37.87)  # 太原市大致坐标

# 批量转换
lons = [112.55, 112.87, 113.12]
lats = [37.87, 38.02, 37.75]
xs, ys = transformer.transform(lons, lats)

# 常用 EPSG 代码:
# 4326  - WGS84 经纬度 (GPS 默认)
# 3857  - Web Mercator (在线地图)
# 4490  - CGCS2000 (中国国家标准)
# 4547  - CGCS2000 / 3度带 Gauss-Kruger CM 114E (山西区域)
# 4480  - Beijing 1954
# 4610  - Xian 1980

8. 可视化制图

静态地图（matplotlib）

import matplotlib.pyplot as plt
import rasterio
from rasterio.plot import show

fig, ax = plt.subplots(figsize=(12, 8))

with rasterio.open("ndvi.tif") as src:
    # 显示栅格
    img = show(src, ax=ax, cmap='RdYlGn', title='NDVI Distribution')
    # 添加颜色条
    im = img.get_images()[0]
    fig.colorbar(im, ax=ax, label='NDVI', shrink=0.6)

# 叠加矢量边界
import geopandas as gpd
boundary = gpd.read_file("boundary.shp")
boundary.boundary.plot(ax=ax, color='black', linewidth=1)

plt.savefig("ndvi_map.png", dpi=300, bbox_inches='tight')
plt.close()

交互式地图（folium）

import folium
import rasterio
from matplotlib import cm, colors

# 创建底图
m = folium.Map(location=[37.87, 112.55], zoom_start=10, tiles='CartoDB positron')

# 添加矢量边界
import geopandas as gpd
gdf = gpd.read_file("boundary.shp").to_crs(epsg=4326)
folium.GeoJson(gdf, name='Boundary').add_to(m)

# 添加标记
folium.Marker([37.87, 112.55], popup='太原').add_to(m)

# 保存为 HTML
m.save("interactive_map.html")

9. 常用分析脚本

技能内置以下可执行脚本（位于 scripts/ 目录）：

| 脚本 | 用途 | 用法 | |------|------|------| | raster_info.py | 快速查看栅格文件信息 | python raster_info.py image.tif | | calc_indices.py | 计算遥感指数（NDVI/NDWI/NDBI） | python calc_indices.py image.tif --index ndvi --nir 4 --red 3 | | clip_raster.py | 用矢量裁剪栅格 | python clip_raster.py image.tif boundary.shp output.tif | | zonal_stats.py | 分区统计 | python zonal_stats.py raster.tif zones.shp --stats mean,sum | | download_satellite_imagery.py | 从 GEE 下载卫星影像 | python download_satellite_imagery.py --place 太原市 --satellite sentinel2 --start 2024-06-01 --end 2024-09-30 --cloud 20 |

地图合规须知

生成任何地图时，必须遵守中国地图合规规则：

• 允许的数据源：腾讯地图、高德地图、百度地图、天地图
• 禁止使用：Google Maps、Apple Maps、Mapbox、OpenStreetMap 直接瓦片
• 底图渲染请配合内置技能 geo-map-compliance-guard 使用
• 交互式地图（folium）默认使用 OpenStreetMap 瓦片，在中国境内使用时需替换为合规底图

注意事项

1. 大文件处理: 遥感影像可能很大，使用 rasterio.open() 的上下文管理器确保文件正确关闭；大文件使用窗口化读取 src.read(window=Window(0, 0, 512, 512))
2. NoData 值: 计算前务必处理 NoData，避免影响统计结果
3. 坐标系一致性: 栅格与矢量做空间运算前，必须确保 CRS 一致
4. 数据类型: 遥感影像常用 uint16，计算指数时需转为 float
5. 内存管理: 大栅格用 rasterio.windows.Window 分块处理
6. 中文路径: Windows 下文件路径含中文时，确保 Python 脚本以 UTF-8 编码运行