在Ubuntu/WSL下用Python+GDAL实现WRF土地利用数据全流程处理

当我们需要提升WRF模拟精度时，使用更精确的土地利用数据是关键。本文将详细介绍如何在Linux环境下，完全基于开源工具链完成从原始数据到WRF可读格式的完整转换流程。

1. 环境准备与数据获取

在开始处理前，我们需要确保系统环境配置正确。对于Ubuntu或WSL用户，建议使用Python 3.8+版本，并安装必要的依赖库。

首先安装基础工具和GDAL：

sudo apt update
sudo apt install -y python3-pip gdal-bin libgdal-dev

然后通过pip安装Python依赖：

pip install numpy gdal

验证GDAL安装是否成功：

gdalinfo --version

获取中国土地利用数据时，可以从权威科研机构网站下载GRID或GeoTIFF格式的数据文件。常见的数据源包括：

• 中国科学院资源环境科学数据中心
• 国家地球系统科学数据中心
• 全球土地覆盖数据集

提示：下载数据时注意检查数据的投影信息和时间戳，确保与你的研究时段匹配。

2. 数据格式转换与投影处理

原始数据可能是各种格式，我们需要将其转换为GDAL能够直接处理的格式。对于GRID格式数据，可以使用以下命令转换为GeoTIFF：

gdal_translate -of GTiff input.grid output.tif

投影转换是数据处理的关键步骤。中国土地利用数据通常使用Krasovsky_1940_Albers投影，而WRF需要WGS84地理坐标系。使用GDAL进行投影转换：

from osgeo import gdal

def reproject_tif(input_file, output_file):
    options = gdal.WarpOptions(
        dstSRS='EPSG:4326',  # WGS84坐标系
        resampleAlg=gdal.GRA_NearestNeighbour
    )
    gdal.Warp(output_file, input_file, options=options)

投影转换后，建议检查结果：

gdalinfo reprojected.tif

输出应包含类似以下信息：

Coordinate System is:
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433]]

3. 土地利用数据重分类

WRF模型支持的土地利用分类体系主要有USGS-24和IGBP-20/21。我们需要将中国土地利用分类映射到这些体系中。

创建重分类映射字典：

classification_map = {
    # 中国分类代码: USGS分类代码
    51: 1,   # 城市用地 → 城市和建筑区
    11: 3,   # 水田 → 灌溉农田和牧场
    31: 7,   # 林地 → 落叶阔叶林
    # 其他分类映射...
}

完整的重分类函数实现：

import numpy as np

def reclassify_landuse(input_array):
    output = np.zeros_like(input_array)
    for src, dst in classification_map.items():
        output[input_array == src] = dst
    return output

处理完成后保存结果：

def save_reclassified(input_file, output_file):
    dataset = gdal.Open(input_file)
    band = dataset.GetRasterBand(1)
    data = band.ReadAsArray()
    
    reclassified = reclassify_landuse(data)
    
    driver = gdal.GetDriverByName('GTiff')
    out_dataset = driver.Create(
        output_file,
        dataset.RasterXSize,
        dataset.RasterYSize,
        1,
        gdal.GDT_Int16
    )
    out_dataset.SetGeoTransform(dataset.GetGeoTransform())
    out_dataset.SetProjection(dataset.GetProjection())
    out_band = out_dataset.GetRasterBand(1)
    out_band.WriteArray(reclassified)
    out_band.SetNoDataValue(255)
    out_dataset.FlushCache()

4. 转换为WRF二进制格式

WRF需要特定的二进制格式(.bil)作为输入。使用GDAL进行格式转换：

gdal_translate -of ENVI -co INTERLEAVE=BSQ reclassified.tif data.bil

这会生成三个文件：

• data.bil (实际数据)
• data.hdr (头文件)
• data.bil.aux.xml (元数据)

检查生成的头文件内容：

ENVI
description = { data.bil }
samples = 6157
lines = 3393
bands = 1
header offset = 0
file type = ENVI Standard
data type = 3
interleave = bsq
byte order = 0
map info = {Geographic Lat/Lon, 1, 1, 66.2922997247769, 54.9785484716575, 0.0117836266203959, 0.0117836266203959, WGS-84}

创建index文件：

type = categorical
category_min = 1
category_max = 24
projection = regular_ll
dx = 0.0117836266203959
dy = 0.0117836266203959
known_x = 1.0
known_y = 1.0
known_lat = 14.9967033487
known_lon = 66.2922997248
wordsize = 1
tile_x = 6157
tile_y = 3393
tile_z = 1
units = "category"
description = "USGS 24-category land use categories"
mminlu = "USGS"
missing_value = 128
iswater = 16
islake = -1
isice = 24
isurban = 1
row_order = top_bottom

最后重命名二进制文件：

mv data.bil 00001-06157.00001-03393

5. 配置WRF使用新数据

将处理好的数据放入WPS_GEOG目录下的适当位置，然后修改GEOGRID.TBL文件：

name = LANDUSEF
priority = 1
dest_type = categorical
abs_path = /path/to/your/data/lucc2015
landmask_water = lucc2015
interp_option = lucc2015:nearest_neighbor
rel_path = lucc2015:lucc2015/

在namelist.wps中指定使用新数据：

&geogrid
 geog_data_res = 'usgs_30s+default','lucc2015+default','lucc2015+default',
/

6. 常见问题排查

Q: GDAL报错"Unable to open EPSG support file" A: 安装proj-data包：

sudo apt install proj-data

Q: 重分类后数据值不正确 A: 检查分类映射表是否完整，可以使用以下代码验证：

unique_values = np.unique(reclassified)
print("结果中的唯一值:", unique_values)

Q: 生成的二进制文件WRF无法识别 A: 确保：

1. 文件命名格式正确
2. index文件参数与数据一致
3. 文件权限设置正确

对于大规模数据处理，可以考虑使用Python多进程加速：

from multiprocessing import Pool

def process_chunk(args):
    # 分块处理函数
    pass

with Pool(processes=4) as pool:
    results = pool.map(process_chunk, chunks)

通过这套完整的开源工具链，我们实现了从原始数据到WRF可用格式的全流程处理，避免了商业软件的依赖，特别适合在Linux服务器和高性能计算环境中使用。