📌 核心数据特性速览

💡 批量下载脚本片段（Python + requests）

import requests
import os

# 根据研究区经纬度范围计算所需带号
def calc_tile_ids(min_lon, max_lon, min_lat, max_lat):
    # 每6度一个带，从-180开始编号
    tiles = []
    for lon in range(int(min_lon//6)*6, int(max_lon//6)*6+1, 6):
        for lat in range(int(min_lat//6)*6, int(max_lat//6)*6+1, 6):
            ns = 'N' if lat >= 0 else 'S'
            zone = int((lon + 180) / 6) + 1
            tiles.append(f"{ns}{abs(lat)}_{lon}_2020LC030")
    return tiles

# 示例：下载中国东部区域（东经110-122，北纬20-42）
tiles = calc_tile_ids(110, 122, 20, 42)
base_url = "http://www.globallandcover.com/data/2020/"
for tile in tiles:
    url = f"{base_url}{tile}.zip"
    r = requests.get(url, stream=True)
    with open(f"{tile}.zip", 'wb') as f:
        for chunk in r.iter_content(chunk_size=8192):
            f.write(chunk)
    print(f"✅ 已下载: {tile}")

🔧 使用GDAL进行批量解压与镶嵌

# 1. 批量解压所有ZIP文件
for f in *.zip; do unzip -o "$f" -d ./extracted/; done

# 2. 使用gdalbuildvrt构建虚拟栅格（避免直接镶嵌的内存问题）
gdalbuildvrt -overwrite -srcnodata 0 -vrtnodata 0 mosaic.vrt ./extracted/*.tif

# 3. 转换为GeoTIFF并压缩（减少存储空间）
gdal_translate -of GTiff -co COMPRESS=LZW -co BIGTIFF=YES mosaic.vrt final_mosaic.tif

# 4. 使用研究区矢量进行裁剪
gdalwarp -cutline study_area.shp -crop_to_cutline -dstnodata 0 -co COMPRESS=LZW study_area_glc.tif final_mosaic.tif

📊 常见重映射方案示例（以2020版分类为例）

import numpy as np
from osgeo import gdal

# 定义重映射字典（原始值: 新值）
remap_dict = {10: 1, 20: 2, 30: 2, 40: 2, 50: 3, 60: 4, 80: 5, 90: 6, 100: 6}

# 读取栅格数据
ds = gdal.Open('study_area_glc.tif')
band = ds.GetRasterBand(1)
arr = band.ReadAsArray()

# 向量化重映射
keys = np.array(list(remap_dict.keys()))
vals = np.array(list(remap_dict.values()))
mapping = np.arange(0, 256)  # 假设像素值范围0-255
for k, v in remap_dict.items():
    mapping[k] = v
reclassified = mapping[arr]

# 写入新文件
driver = gdal.GetDriverByName('GTiff')
out_ds = driver.Create('reclassified.tif', ds.RasterXSize, ds.RasterYSize, 1, gdal.GDT_Byte)
out_ds.SetProjection(ds.GetProjection())
out_ds.SetGeoTransform(ds.GetGeoTransform())
out_band = out_ds.GetRasterBand(1)
out_band.WriteArray(reclassified)
out_band.SetNoDataValue(0)
out_ds.FlushCache()

#!/usr/bin/env python3
"""
GlobeLand30全自动处理流水线
适用场景：批量处理多个年份、多个研究区
"""
import os, glob, subprocess
import numpy as np
from osgeo import gdal

# ===== 配置参数 =====
YEARS = ['2000', '2010', '2020']
STUDY_AREA = 'my_study_area.shp'
OUTPUT_DIR = './output/'
REMAP_DICT = {10: 1, 20: 2, 30: 2, 40: 2, 50: 3, 60: 4, 80: 5, 90: 6, 100: 6}

def download_tiles(year, lat_range, lon_range):
    # 实现下载逻辑（略，参考2.2节）
    pass

def mosaic_and_clip(year):
    # 镶嵌并裁剪
    subprocess.run(f'gdalbuildvrt mosaic_{year}.vrt ./{year}/*.tif', shell=True)
    subprocess.run(f'gdalwarp -cutline {STUDY_AREA} -crop_to_cutline -co COMPRESS=LZW {year}_clipped.tif mosaic_{year}.vrt', shell=True)

def reclassify(year):
    ds = gdal.Open(f'{year}_clipped.tif')
    arr = ds.GetRasterBand(1).ReadAsArray()
    mapping = np.arange(256)
    for k, v in REMAP_DICT.items():
        mapping[k] = v
    reclass = mapping[arr]
    # 写入文件
    driver = gdal.GetDriverByName('GTiff')
    out = driver.Create(f'{OUTPUT_DIR}{year}_reclass.tif', ds.RasterXSize, ds.RasterYSize, 1, gdal.GDT_Byte)
    out.SetProjection(ds.GetProjection())
    out.SetGeoTransform(ds.GetGeoTransform())
    out.GetRasterBand(1).WriteArray(reclass)
    out.FlushCache()

# 主循环
for year in YEARS:
    print(f"🚀 正在处理 {year} 年数据...")
    download_tiles(year, lat_range=(20, 42), lon_range=(110