Python：如何处理WRF投影(LCC, 兰伯特投影)？

2025年8月14日
| 没有评论

01 问题和说明

1.1 问题

目前需要解决的问题是：

  1. 如何将WRF输出的两个nc文件(变量均为T2，分辨率分别为9000m和3000m, 文件名分别为: wrfout_d01_2008-01-01_T2.nc和wrfout_d02_2008-01-01_T2.nc)输出为LCC(Lambert Conformal Conic, 兰伯特等角投影)投影坐标系的GeoTIFF文件？
  2. 如何将GLASS的LAI产品(hdf4格式,规则格网,全球)输出为LCC坐标系的GeoTIFF文件？
  3. 如何将MODIS的土地覆盖产品(MCD12C1, hdf4格式, 规则格网, 全球, 与GLASSLAI类似)输出为LCC坐标系的GeoTIFF文件？
  4. 如何将WGS84地理坐标系的聚集指数GeoTIFF文件(全球, 规则格网, 与GLASS和MODIS类似但非hdf文件)输出为LCC坐标系的GeoTIFF文件？
  5. 如何将WGS84地理坐标系的土壤类型GeoTIFF文件(全球, 规则格网, 与GLASS和MODIS类似但非hdf文件)输出为LCC坐标系的GeoTIFF文件？

要求：

  1. 所有结果包含两种分辨率，9000m和3000m
  2. 相同分辨率的tiff文件的行列数均要保持一致，且与问题1中WRF输出的nc文件中相同分辨率的栅格矩阵的行列数相同。即对于指定分辨率例如9000m所对应的nc文件，其中T2变量的行列数为224行236列，那么对于MODIS、GLASS等产品应当做重投影/GLT校正(为LCC坐标系)、裁剪掩膜、重采样等操作，使所有产品与nc文件进行空间上的统一和对齐.

1.2 说明

1.2.1 WRF输出的nc文件

问题1中的nc文件形如：

可以发现，其中的LAT和LONG经纬度数据集均为二维数据集，查看其中的栅格矩阵中可以发现并非规则格网(其实很显然如果是规则格网大概率不需要提供两个二维的坐标集，只需要提供两个一维坐标集甚至连一维都不需要，通过角点等地理信息定义投影参数和仿射系数即可)

因此对于此处的nc文件，是最麻烦的事情，处理好了这个坐标问题，也就基本上解决了所有问题的核心问题：投影问题。

1.2.2 GLASS-LAI产品

LAI产品形如：

查看上面的LAI的基本信息基本可以简单计算和判断，这是一个规则格网，坐标系为标准的WGS84坐标系，左上角点为(-180°，90°)，分辨率为0.05°×0.05°。为了确保无误，应该查看元数据StructMetadata.0(如下图)。此外，对于读取的栅格矩阵还应进行无效值处理和比例缩放(依据valid_range和scale_factor处理)。

1.2.3 MCD12C1-土地覆盖产品

MCD12C1产品形如：

土地覆盖产品的坐标与GLASS类似，都是一样的处理，按照规则格网来进行。

1.2.4 聚集指数和土壤类型(GeoTIFF文件)

这里给我的是WGS84坐标系的全球范围的GeoTIFF文件，形如：


02 思路

思路1：

由于最早的要求是空间上统一和对齐即可，也就是只要保证相同分辨率下所有结果的行列数一致即可，并没有要求与最初的nc文件的行列数相同。因此最开始的想法是进行重投影/GLT校正(此处的重投影/GLT校正本质就是将两类坐标建立数学关联(例如通过多项式等),然后一一匹配，对于这方面的底层原理参考附录)，但是GLT校正没有办法保证与源栅格矩阵的行列数保持一致。

但是不管怎么说，重投影/GLT校正仍是遥感数据处理中处理难度较大的一个小部分，在附录中会对该思路稍作解释。

思路2：

WRF输出的nc文件，其实不仅仅是本例中的nc文件，只要是WRF的输入输出，其投影参数都是一致的。即：


其中：

CEN_LAT：投影中心的纬度坐标
CEN_LON：投影中心的经度坐标
TRUELAT1：第一条标准纬线的纬度
TRUELAT2：第二条标准纬线的纬度
MAP_PROJ = 1：指代LCC投影
POLE_LAT = 90.0f：北极点的纬度位置
POLE_LON = 0.0f：北极点的经度位置
MOAD_CEN_LAT = 46.50001f：投影中心的纬度坐标(WRF中同一次模拟不同区域的投影中心的纬度坐标均为此)
STAND_LON = 10.0f：标准经线/中央经线
MAP_PROJ_CHAR = “Lambert Conformal”：兰伯特等角投影

使用proj4语法表示为:

from pyproj import CRS
lcc_proj4 = '+proj=lcc +lat_1=30.0 +lat_2=60.0 +lat_0=46.50001 +lon_0=10.0 +a=6370000 +b=6370000 +units=m +towgs84=0,0,0 +no_defs'
lcc_srs = CRS(lcc_proj4)

因此对于WRF输出的nc文件，可以采用类似规则格网的方法去进行，区别在于规则格网是在地理坐标系下进行(坐标是经纬度)，而本例中是在投影坐标系下进行（坐标是投影坐标,单位m），因此处理会存在不同。

后续除附录外，下文没有特别说明都是针对思路2进行说明和论述。

03 WRF的LCC投影中容易混淆的点

参考：

https://fabienmaussion.info/2018/01/06/wrf-projection/
https://www.pkrc.net/wrf-lambert.html

如果你不想去研究WRF中LCC投影(对于标准的LCC无需此处理，仅只针对WRF中的LCC投影)涉及的内容，那么只需要记住上面这幅图：

(仅针对WRF中的LCC投影)如果你需要将LCC投影的GeoTIFF文件转化为其他投影坐标系，应该将LCC投影的GeoTIFF文件输出为r=6370km（即长短半轴a=6370km,b=6370km）WGS84坐标系(非标准datum)的GeoTIFF文件，然后直接将这个非标准的WGS84坐标系直接重新定义投影为标准的WGS84坐标系(即直接删除原来的非标准WGS84坐标系，然后赋值上新的标准WGS84坐标系)，如果目标是转化为WGS84坐标系，那么到这里就结束了。如果还要转化为其他投影坐标系或者地理坐标系，就在此基础上继续进行投影转换。

不能直接将WRF的LCC投影直接转换为WGS84坐标系或者其他投影或者地理坐标系，因为WRF的LCC投影中所基于的椭球体是简单的球体即r=6370km，而并非如WGS84坐标系中的椭球体。因此直接将WRF的LCC投影转化为WGS84坐标系在python例如gdal、pyproj等中由于二者的基准面不一致也就是椭球体不一致，因此会多一步基准面转换(但是这一步是多余且完全不需要的，且会导致巨大的误差)。

本博客中只涉及了LCC投影与WGS84之间的转换。其中，

GLASS产品虽然是规则格网，但是按照上述的转换步骤，应该是先按照原本的规则格网定义为标准的WGS84坐标系，再重新定义为简化的r=6370km的WGS84坐标系，所以不如直接定义为简化的r=6370km的WGS84坐标系(尽管它原本不是标准WGS84坐标系)。

MCD12C1产品、聚集指数和也是类似上方的处理。

对于WRF输出的nc文件(T2)，按照思路2的方法，是直接定义为LCC投影，因此基本不涉及LCC与其他坐标系的转换。

04 方法

按照思路2的方法。

下述代码存在自定义函数，具体定义如下：

# @Author  : ChaoQiezi
# @Time    : 2025/8/1 下午5:33
# @Email   : chaoqiezi.one@qq.com
# @Wechat  : GIS茄子
# @FileName: utils

"""
This script is used to 存放函数
"""

import numpy as np
from osgeo import gdal, osr
from pyproj import Transformer
gdal.UseExceptions()

# out_bound = [-1061649, -1022770, 1062350, 993229]
out_bound_9000m = [-1000000, -1000000, 1000000, 1000000]
out_bound_3000m = [-380000, -340000, 380000, 280000]


def img_glt(src_var, src_lon, src_lat, out_path, out_res, resample_alg=None, out_type=None, src_srs=None, dst_srs=None,
            out_bound=None):
    """
    基于二维的src_lon(纬度数据集)和src_lat（经度数据集）对src_var（目标数据集）进行GLT校正并输出为tiff文件
    :param src_var: 需要重投影/GLT校正的目标数据集, 要求三维np数组且shape=(波段数, 行数, 列数), 或二维数组且shape=(行数, 列数)
    :param src_lon: 经度数据集(二维np数组)
    :param src_lat: 纬度数据集(二维np数组)
    :param out_path: 重投影/GLT校正结果文件的输出路径(.tif)
    :param out_res: 输出分辨率(单位取决于dst_srs)
    :param resample_alg: 重采样算法(使用gdal.GRA_*), 默认为最近邻
    :param out_type: 输出像元值的数据类型, 默认依据src_var判断整型或者浮点
    :param src_srs: 输入栅格数据集的坐标系（lon和lat的坐标系,一般都是WGS84坐标系), 默认WGS84坐标系
    :param dst_srs: 输出栅格数据集的坐标系(重投影之后的坐标系), 默认WGS84坐标系
    :return: None
    """

    if resample_alg is None:
        resample_alg = gdal.GRA_NearestNeighbour

    if out_type is None:
        if np.issubdtype(src_var.dtype, np.integer):
            out_type = gdal.GDT_Int16
        elif np.issubdtype(src_var.dtype, np.floating):
            out_type = gdal.GDT_Float32
        else:
            out_type = gdal.GDT_Float32

    if len(src_var.shape) == 2:
        src_var = np.expand_dims(src_var, axis=0)  # (new_axis, rows, cols)

    # 定义输入栅格数据集的坐标系
    if src_srs is None:
        src_srs = osr.SpatialReference()
        src_srs.ImportFromEPSG(4326)
    # 定义输出栅格数据集的坐标系
    if dst_srs is None:
        dst_srs = osr.SpatialReference()
        dst_srs.ImportFromEPSG(4326)
        dst_srs.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)
    src_srs.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)  # 使用传统的(lon, lat)顺序, 下同
    dst_srs.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)
    # 获取基本信息

    band_count, rows, cols = src_var.shape

    # 为目标、经纬度分别创建临时tiff文件(内存中创建)
    var_mem_path = '/vsimem/var.tif'  # 内存中临时路径, 下同
    lon_mem_path = '/vsimem/lon.tif'
    lat_mem_path = '/vsimem/lat.tif'
    tiff_driver = gdal.GetDriverByName('GTiff')  # 创建TIFF驱动器
    img_var = tiff_driver.Create(var_mem_path, cols, rows, band_count, out_type)  # 创建tiff文件, 下同
    img_lon = tiff_driver.Create(lon_mem_path, cols, rows, 1, out_type)
    img_lat = tiff_driver.Create(lat_mem_path, cols, rows, 1, out_type)
    for band_ix in range(band_count):
        img_var.GetRasterBand(band_ix + 1).WriteArray(src_var[band_ix, :, :])  # 写入波段数据, lon和lat类似
    img_lon.GetRasterBand(1).WriteArray(src_lon)
    img_lat.GetRasterBand(1).WriteArray(src_lat)
    img_var = img_lon = img_lat = None  # 释放资源, 将缓冲区的数据全部写入

    # 更新元数据-关联地理定位
    img_var = gdal.Open(var_mem_path, gdal.GA_Update)
    """
    使用img_var=None后又使用gdal.Open打开似乎是重复冗余的操作, 实际并不是:
    首要原因是前面是在写入数据集, 而后面是在更新元数据 -- 写入和更新是不应该放在一起否则会混淆报错.
    这是因为写入之后本身就会生成元数据出来, 但是后面又有一个setMetadata更新GEOLOCATION会使得信息写入
    混乱最终导致后续重投影失败
    因此这是必要的.
    """
    img_var.SetMetadata({
        'SRS': src_srs.ExportToWkt(),  # X_DATASET和Y_DATASET所使用的坐标系
        'X_DATASET': lon_mem_path,  # 包含X坐标（通常是经度）的栅格数据集的路径
        'X_BAND': '1',  # 指定使用X坐标栅格数据集中的第一个波段
        'Y_DATASET': lat_mem_path,  # 包含Y坐标（通常是经度）的栅格数据集的路径
        'Y_BAND': '1',  # 指定使用Y坐标栅格数据集中的第一个波段
        'PIXEL_OFFSET': '0',  # X方向的偏移量
        'LINE_OFFSET': '0',  # Y方向的偏移量
        'PIXEL_STEP': '1',  # X方向的步长
        'LINE_STEP': '1'  # Y方向的步长
    }, 'GEOLOCATION')
    """
    对于传入字典中的最后四个变量：PIXEL_OFFSET、LINE_OFFSET、PIXEL_STEP和LINE_STEP。
    应当这样理解，它类似于于GLT校正,但是不同于ENVI IDL中的GLT校正工具是通过地面控制点，这里传入的每一个点都类似于控制点。
    这里传入X和Y的坐标类似于ENVI中的经纬度查找表，那我们知道，如果将所有的坐标点参与GLT校正那么在分辨率高地理覆盖范围大
    的情况下, 计算量会非常大运行时间会非常长，因此这里可以供你选择：
    假如从第200行300列开始从右隔4个像元采集一个像元坐标参与GLT校正，向下隔2个像元采集一个像元坐标参与GLT校正。
    那么应该设置为:
    'PIXEL_OFFSET': '300',  # X方向的偏移量
    'LINE_OFFSET': '200',  # Y方向的偏移量
    'PIXEL_STEP': '4',  # X方向的步长
    'LINE_STEP': '2'  # Y方向的步长

    这样的话大部分区域的像元都参与了运算而且计算量也下来了。
    """
    img_var.SetProjection(src_srs.ExportToWkt())  # 设置输入栅格数据集的坐标系
    img_var = None

    # 重投影/GLT校正
    warp_options = gdal.WarpOptions(
        dstSRS=dst_srs,
        outputBounds=out_bound,
        dstNodata=-999,
        xRes=out_res,
        yRes=out_res,
        resampleAlg=resample_alg,
        transformerOptions=['METHOD=GEOLOC_ARRAY'],
        targetAlignedPixels=True
    )
    gdal.Warp(out_path, var_mem_path, options=warp_options)

    # 释放内存文件
    gdal.Unlink(var_mem_path)
    gdal.Unlink(lon_mem_path)
    gdal.Unlink(lat_mem_path)


def img_warp(src_var, src_srs, src_geo_transform, dst_srs, out_res, out_path, out_bound, out_type=None,
             resamlpe_alg=None, nodata_value=None):
    if out_type is None:
        out_type = gdal.GDT_Float32

    if resamlpe_alg is None:
        resamlpe_alg = gdal.GRA_NearestNeighbour

    rows, cols = src_var.shape

    var_mem_path = '/vsimem/var.tif'
    # 创建输入栅格数据集(内存中创建)
    tiff_driver = gdal.GetDriverByName('GTiff')
    src_lai = tiff_driver.Create(var_mem_path, cols, rows, 1, out_type)
    src_lai.SetProjection(src_srs)  # 设置坐标系
    src_lai.SetGeoTransform(src_geo_transform)  # 设置仿射系数
    src_lai.GetRasterBand(1).WriteArray(src_var)
    src_lai = None

    # 重投影
    warp_options = gdal.WarpOptions(
        dstSRS=dst_srs,
        outputBounds=out_bound,
        xRes=out_res,
        yRes=out_res,
        srcNodata=nodata_value,
        dstNodata=nodata_value,
        resampleAlg=resamlpe_alg,
        targetAlignedPixels=True,
        outputType=out_type
    )
    """
    targetAlignedPixels=True的说明
    保证对齐, 因此outputBounds传输的边界限制是基础, 实际会在边界上多或者少一个像元列或者行,
    但是保证了所有的产品使用这一设置都是一致的边界范围和行列数.
    """
    gdal.Warp(out_path, var_mem_path, options=warp_options)

    gdal.Unlink(var_mem_path)  # 释放内存中创建的虚拟tiff文件


def cal_geo_transform(cen_lon, cen_lat, rows, cols, x_res, y_res, src_srs, dst_srs):
    """
    基于给定的投影中心的地理坐标系计算仿射系数矩阵
    :param cen_lon:
    :param cen_lat:
    :param rows:
    :param cols:
    :param x_res:
    :param y_res:
    :return:
    """

    # 投影中心的坐标转换(WGS84 --> LCC)
    wgs2lcc_transformer = Transformer.from_crs(src_srs, dst_srs, always_xy=True)
    cen_x, cen_y = wgs2lcc_transformer.transform(cen_lon, cen_lat)
    # 计算左上角像元(第一行第一列)的中心坐标
    ul_x_center = cen_x - (cols - 1) / 2 * x_res
    ul_y_center = cen_x + (rows - 1) / 2 * y_res
    # 计算左上角像元(第一行第一列)的左上角角点坐标
    ul_x_corner = ul_x_center - x_res / 2
    ul_y_corner = ul_y_center + y_res / 2
    # 构建GeoTransform
    geo_transform = [ul_x_corner, x_res, 0, ul_y_corner, 0, -y_res]

    return geo_transform


def write_tiff(out_path, img_band, proj_wkt, geo_transform):
    """
    将输入的栅格矩阵输出为GeoTIFF文件
    :param out_path:
    :param img_band:
    :param proj_wkt:
    :param geo_transform:
    :return:
    """
    rows, cols = img_band.shape

    driver = gdal.GetDriverByName('GTiff')
    out_img = driver.Create(out_path, cols, rows, 1, gdal.GDT_Float32)
    out_img.GetRasterBand(1).WriteArray(img_band)
    out_img.SetProjection(proj_wkt)
    out_img.SetGeoTransform(geo_transform)
    out_img.FlushCache()
    out_img = None

4.1 处理WRF输出的NC文件(T2)

# @Author  : ChaoQiezi
# @Time    : 2025/8/13 下午7:21
# @Email   : chaoqiezi.one@qq.com
# @Wechat  : GIS茄子
# @FileName: process_T2

"""
This script is used to 处理T2数据集
"""

import os
import numpy as np
import xarray as xr
from pyproj import CRS
from osgeo import gdal
gdal.UseExceptions()

from utils import cal_geo_transform, write_tiff

# 准备
d01_path = r"E:\Datasets\Objects\reproject_europe\wrfout_d01_2008-01-01_T2.nc"
d02_path = r"E:\Datasets\Objects\reproject_europe\wrfout_d02_2008-01-01_T2.nc"
out_dir = r'E:\Datasets\Objects\reproject_europe\Result2'
rows_d01, cols_d01 = 224, 236
rows_d02, cols_d02 = 201, 252
out_d01_res = 9000  # 输出分辨率9000m
out_d02_res = 3000
cen_lon, cen_lat = 10.0, 46.50001
# 定义LCC(兰伯特投影坐标系)
lcc_proj4 = '+proj=lcc +lat_1=30.0 +lat_2=60.0 +lat_0=46.50001 +lon_0=10.0 +a=6370000 +b=6370000 +units=m +towgs84=0,0,0 +no_defs'
lcc_srs = CRS(lcc_proj4)
# 定义等经纬度格网投影(WGS84坐标系)
wgs84_srs = CRS('+proj=latlong +a=6370000 +b=6370000')  # 一定要添加a和b

# 计算d01的仿射系数
geo_transform = cal_geo_transform(cen_lon, cen_lat, rows_d01, cols_d01, out_d01_res, out_d01_res,
                                  wgs84_srs, lcc_srs)
# 读取T2变量并输出为tiff文件
da = xr.open_dataset(d01_path)
img = np.flipud(da['T2'].values[0, :, :])
out_path = os.path.join(out_dir, os.path.basename(d01_path).replace('.nc', '.tif'))
write_tiff(out_path, img, lcc_srs.to_wkt(), geo_transform)

# 计算d02的仿射系数
geo_transform = cal_geo_transform(cen_lon, cen_lat, rows_d02, cols_d02, out_d02_res, out_d02_res,
                                  wgs84_srs, lcc_srs)
# 读取T2变量并输出为tiff文件
da = xr.open_dataset(d02_path)
img = np.flipud(da['T2'].values[0, :, :])
out_path = os.path.join(out_dir, os.path.basename(d02_path).replace('.nc', '.tif'))
write_tiff(out_path, img, lcc_srs.to_wkt(), geo_transform)

4.2 处理GLASS-LAI产品

# @Author  : ChaoQiezi
# @Time    : 2025/8/13 下午9:55
# @Email   : chaoqiezi.one@qq.com
# @Wechat  : GIS茄子
# @FileName: process_glass.py

"""
This script is used to 
"""


import os
from os import makedirs
import numpy as np
from pyproj import CRS
from pyhdf.SD import SD, SDC
from glob import glob
from osgeo import gdal
gdal.UseExceptions()


# 准备
in_dir = r"E:\Datasets\Objects\reproject_europe\GLASSLAI"
out_dir = r'E:\Datasets\Objects\reproject_europe\Result2\GLASSLAI'
ref_paths = [r"E:\Datasets\Objects\reproject_europe\Result2\wrfout_d01_2008-01-01_T2.tif",
             r"E:\Datasets\Objects\reproject_europe\Result2\wrfout_d02_2008-01-01_T2.tif"]
start_year, end_year = 2008, 2012
var_name = 'LAI'
postfix_dict = {9000: 'd01', 3000: 'd02'}
src_gt = [-180, 0.05, 0, 90, 0, -0.05]
# 定义等经纬度格网投影(WGS84坐标系)
wgs84_srs = CRS('+proj=latlong +a=6370000 +b=6370000')  # 一定要添加a和b
# 定义LCC(兰伯特投影坐标系)
lcc_proj4 = '+proj=lcc +lat_1=30.0 +lat_2=60.0 +lat_0=46.50001 +lon_0=10.0 +a=6370000 +b=6370000 +units=m +towgs84=0,0,0 +no_defs'
lcc_srs = CRS(lcc_proj4)

# 检索HDF文件
for cur_ref_path in ref_paths:
    # 获取参考文件的地理信息
    ref_ds = gdal.Open(cur_ref_path)
    ref_gt = ref_ds.GetGeoTransform()
    width = ref_ds.RasterXSize
    height = ref_ds.RasterYSize
    ref_srs_wkt = ref_ds.GetProjection()
    # 计算输出边界
    xmin = ref_gt[0]
    ymax = ref_gt[3]
    xmax = xmin + width * ref_gt[1]
    ymin = ymax + height * ref_gt[5]
    out_bound = (xmin, ymin, xmax, ymax)
    # 获取目标分辨率
    ref_res = ref_gt[1]

    for cur_year in range(start_year, end_year + 1):
        # 检索当前年份的LAI-HDF文件
        cur_year_paths = glob(os.path.join(in_dir, str(cur_year), 'GLASS*.hdf'))
        cur_out_dir = os.path.join(out_dir, str(cur_year))
        makedirs(cur_out_dir, exist_ok=True)

        for cur_path in cur_year_paths:
            # 读取变量
            h4 = SD(cur_path, SDC.READ)
            img_var = np.float32(h4.select(var_name)[:])
            h4.end()

            # 无效值处理
            img_var[(img_var < 0) | (img_var > 100)] = np.nan
            # 比例缩放
            img_var *= 0.1

            # 创建源tiff文件
            var_mem_path = '/vsimem/var.tif'
            rows, cols = img_var.shape
            # 创建输入栅格数据集(内存中创建)
            tiff_driver = gdal.GetDriverByName('GTiff')
            src_var = tiff_driver.Create(var_mem_path, cols, rows, 1, gdal.GDT_Float32)
            src_var.SetProjection(wgs84_srs.to_wkt())  # 设置坐标系
            src_var.SetGeoTransform(src_gt)  # 设置仿射系数
            src_var.GetRasterBand(1).WriteArray(img_var)
            src_var = None

            warp_options = gdal.WarpOptions(
                format='GTiff',
                dstSRS=ref_srs_wkt,  # 目标投影, LCC
                xRes=ref_res,  # 输出X分辨率
                yRes=ref_res,  # 输出Y分辨率
                outputBounds=out_bound,  # 输出边界范围
                resampleAlg=gdal.GRA_Cubic,  # 重采样算法
                dstNodata=-999
            )
            cur_out_filename = os.path.basename(cur_path).replace('.hdf', '_{}.tif'.format(postfix_dict[ref_res]))
            cur_out_path = os.path.join(cur_out_dir, cur_out_filename)
            ds = gdal.Warp(cur_out_path,
                           var_mem_path,
                           options=warp_options)

            print('处理: {}'.format(cur_out_filename))

    print(f'参考标准({os.path.join(cur_ref_path)})下文件输出完成.')

4.3 处理MCD12C1产品

# @Author  : ChaoQiezi
# @Time    : 2025/8/13 下午8:45
# @Email   : chaoqiezi.one@qq.com
# @Wechat  : GIS茄子
# @FileName: process_mcd12c1.py

"""
This script is used to 
"""

import os
from pyproj import CRS
from pyhdf.SD import SD, SDC
from glob import glob
from osgeo import gdal
gdal.UseExceptions()


# 准备
in_dir = r"E:\Datasets\Objects\reproject_europe\MCD12C1"
out_dir = r'E:\Datasets\Objects\reproject_europe\Result2\MCD12C1'
ref_paths = [r"E:\Datasets\Objects\reproject_europe\Result2\wrfout_d01_2008-01-01_T2.tif",
             r"E:\Datasets\Objects\reproject_europe\Result2\wrfout_d02_2008-01-01_T2.tif"]
land_sover_name = 'Majority_Land_Cover_Type_1'
postfix_dict = {9000: 'd01', 3000: 'd02'}
src_gt = [-180, 0.05, 0, 90, 0, -0.05]
# 定义等经纬度格网投影(WGS84坐标系)
wgs84_srs = CRS('+proj=latlong +a=6370000 +b=6370000')  # 一定要添加a和b
# 定义LCC(兰伯特投影坐标系)
lcc_proj4 = '+proj=lcc +lat_1=30.0 +lat_2=60.0 +lat_0=46.50001 +lon_0=10.0 +a=6370000 +b=6370000 +units=m +towgs84=0,0,0 +no_defs'
lcc_srs = CRS(lcc_proj4)

# 检索HDF文件
hdf_paths = glob(os.path.join(in_dir, 'MCD12C1*.hdf'))
for cur_ref_path in ref_paths:
    # 获取参考文件的地理信息
    ref_ds = gdal.Open(cur_ref_path)
    ref_gt = ref_ds.GetGeoTransform()
    width = ref_ds.RasterXSize
    height = ref_ds.RasterYSize
    ref_srs_wkt = ref_ds.GetProjection()
    # 计算输出边界
    xmin = ref_gt[0]
    ymax = ref_gt[3]
    xmax = xmin + width * ref_gt[1]
    ymin = ymax + height * ref_gt[5]
    out_bound = (xmin, ymin, xmax, ymax)
    # 获取目标分辨率
    ref_res = ref_gt[1]

    for cur_path in hdf_paths:
        h4 = SD(cur_path, SDC.READ)
        land_cover = h4.select(land_sover_name)[:]
        h4.end()

        # 创建源tiff文件
        var_mem_path = '/vsimem/var.tif'
        rows, cols = land_cover.shape
        # 创建输入栅格数据集(内存中创建)
        tiff_driver = gdal.GetDriverByName('GTiff')
        src_var = tiff_driver.Create(var_mem_path, cols, rows, 1, gdal.GDT_Byte)
        src_var.SetProjection(wgs84_srs.to_wkt())  # 设置坐标系
        src_var.SetGeoTransform(src_gt)  # 设置仿射系数
        src_var.GetRasterBand(1).WriteArray(land_cover)
        src_var = None

        warp_options = gdal.WarpOptions(
            format='GTiff',
            dstSRS=ref_srs_wkt,  # 目标投影, LCC
            xRes=ref_res,  # 输出X分辨率
            yRes=ref_res,  # 输出Y分辨率
            outputBounds=out_bound,  # 输出边界范围
            resampleAlg=gdal.GRA_NearestNeighbour,  # 重采样算法
            dstNodata=255
        )
        cur_out_filename = os.path.basename(cur_path).replace('.hdf', '_{}.tif'.format(postfix_dict[ref_res]))
        cur_out_path = os.path.join(out_dir, cur_out_filename)
        ds = gdal.Warp(cur_out_path,
                       var_mem_path,
                       options=warp_options)

        print('处理: {}'.format(cur_out_filename))

    print(f'参考标准({os.path.join(cur_ref_path)})下文件输出完成.')

4.4 处理聚集指数和土壤类型GeoTIFF文件

# @Author  : ChaoQiezi
# @Time    : 2025/8/13 下午10:14
# @Email   : chaoqiezi.one@qq.com
# @Wechat  : GIS茄子
# @FileName: process_tiff.py

"""
This script is used to 
"""

import os
import numpy as np
from osgeo import gdal, osr
from pyproj import CRS

gdal.UseExceptions()

# 准备
clumping_index_path = r"E:\Datasets\Objects\reproject_europe\global_clumping_index_Reprojec_005D.tif"
soil_type_path = r"E:\Datasets\Objects\reproject_europe\USDA_soil_type_005D.tif"
out_dir = r'E:\Datasets\Objects\reproject_europe\Result2'
postfix_dict = {9000: 'd01', 3000: 'd02'}
ref_paths = [r"E:\Datasets\Objects\reproject_europe\Result2\wrfout_d01_2008-01-01_T2.tif",
             r"E:\Datasets\Objects\reproject_europe\Result2\wrfout_d02_2008-01-01_T2.tif"]
wgs84_srs = CRS('+proj=latlong +a=6370000 +b=6370000')  # 一定要添加a和b

# 处理聚集指数
for cur_ref_path in ref_paths:
    # 获取参考文件的地理信息
    ref_ds = gdal.Open(cur_ref_path)
    ref_gt = ref_ds.GetGeoTransform()
    width = ref_ds.RasterXSize
    height = ref_ds.RasterYSize
    ref_srs_wkt = ref_ds.GetProjection()
    # 计算输出边界
    xmin = ref_gt[0]
    ymax = ref_gt[3]
    xmax = xmin + width * ref_gt[1]
    ymin = ymax + height * ref_gt[5]
    out_bound = (xmin, ymin, xmax, ymax)
    # 获取目标分辨率
    ref_res = ref_gt[1]

    # 读取源文件的投影、地理转换信息和栅格矩阵
    src_ds = gdal.Open(clumping_index_path)
    # projection = src_ds.GetProjection()
    # geotransform = src_ds.GetGeoTransform()
    band = src_ds.GetRasterBand(1)
    original_data = np.float32(band.ReadAsArray())
    nodata_value = band.GetNoDataValue()
    # 比例缩放和无效值处理
    original_data[original_data == nodata_value] = np.nan
    scaled_data = original_data * 0.01
    # 创建一个新的栅格文件用于写入结果
    var_mem_path = '/vsimem/var.tif'
    driver = gdal.GetDriverByName('GTiff')
    rows, cols = original_data.shape
    mem_ds = driver.Create(var_mem_path, cols, rows, 1, gdal.GDT_Float32)
    mem_ds.SetProjection(wgs84_srs.to_wkt())
    mem_ds.SetGeoTransform([-180, 0.05, 0, 90, 0, -0.05])
    mem_ds.GetRasterBand(1).WriteArray(scaled_data)
    mem_ds.GetRasterBand(1).SetNoDataValue(np.nan)
    mem_ds = None

    # 重投影等
    warp_options = gdal.WarpOptions(
        format='GTiff',
        dstSRS=ref_srs_wkt,  # 目标投影, LCC
        xRes=ref_res,  # 输出X分辨率
        yRes=ref_res,  # 输出Y分辨率
        outputBounds=out_bound,  # 输出边界范围
        resampleAlg=gdal.GRA_Cubic,  # 重采样算法
        dstNodata=-999
    )
    cur_out_filename = os.path.basename(clumping_index_path).replace('.tif', '_{}.tif'.format(postfix_dict[ref_res]))
    cur_out_path = os.path.join(out_dir, cur_out_filename)
    gdal.Warp(cur_out_path,
              var_mem_path,
              options=warp_options)

    print('处理: {} ({})'.format(cur_out_filename, postfix_dict[ref_res]))

# 处理土壤类型
for cur_ref_path in ref_paths:
    # 获取参考文件的地理信息
    ref_ds = gdal.Open(cur_ref_path)
    ref_gt = ref_ds.GetGeoTransform()
    width = ref_ds.RasterXSize
    height = ref_ds.RasterYSize
    ref_srs_wkt = ref_ds.GetProjection()
    # 计算输出边界
    xmin = ref_gt[0]
    ymax = ref_gt[3]
    xmax = xmin + width * ref_gt[1]
    ymin = ymax + height * ref_gt[5]
    out_bound = (xmin, ymin, xmax, ymax)
    # 获取目标分辨率
    ref_res = ref_gt[1]

    # 读取源文件的投影、地理转换信息和栅格矩阵
    src_ds = gdal.Open(soil_type_path)
    # projection = src_ds.GetProjection()
    # geotransform = src_ds.GetGeoTransform()
    band = src_ds.GetRasterBand(1)
    soil_type = np.float32(band.ReadAsArray())
    nodata_value = band.GetNoDataValue()
    # 创建一个新的栅格文件用于写入结果
    var_mem_path = '/vsimem/var.tif'
    driver = gdal.GetDriverByName('GTiff')
    rows, cols = soil_type.shape
    mem_ds = driver.Create(var_mem_path, cols, rows, 1, gdal.GDT_Int8)
    mem_ds.SetProjection(wgs84_srs.to_wkt())
    mem_ds.SetGeoTransform([-180, 0.05, 0, 90, 0, -0.05])
    mem_ds.GetRasterBand(1).WriteArray(soil_type)
    mem_ds.GetRasterBand(1).SetNoDataValue(nodata_value)
    mem_ds = None

    # 重投影等
    warp_options = gdal.WarpOptions(
        format='GTiff',
        dstSRS=ref_srs_wkt,  # 目标投影, LCC
        xRes=ref_res,  # 输出X分辨率
        yRes=ref_res,  # 输出Y分辨率
        outputBounds=out_bound,  # 输出边界范围
        resampleAlg=gdal.GRA_NearestNeighbour,  # 重采样算法
    )
    cur_out_filename = os.path.basename(soil_type_path).replace('.tif', '_{}.tif'.format(postfix_dict[ref_res]))
    cur_out_path = os.path.join(out_dir, cur_out_filename)
    gdal.Warp(cur_out_path,
              var_mem_path,
              options=warp_options)

    print('处理: {} ({})'.format(cur_out_filename, postfix_dict[ref_res]))

05 附录

5.1 重投影/GLT校正

之前，我曾用IDL对MODIS GRID产品进行过重投影：

  1. https://blog.csdn.net/m0_63001937/article/details/133977692?spm=1011.2415.3001.5331
  2. https://blog.csdn.net/m0_63001937/article/details/134238529?spm=1011.2415.3001.5331

最关键的部分即：

;+
;   函数用途:
;       基于经纬度数据集对目标数据集进行几何校正(重投影-WGS84)
;   函数参数:
;       target: 待校正的目标数据集
;       lon: 对应目标数据集的经度数据集
;       lat: 对应目标数据集的纬度数据集
;       out_res: 输出的分辨率(°)
;       target_warped: 校正后的目标数据集
;       geo_info: 校正后的目标数据集对应的地理结构体
;       degree<关键字参数: 5>: 多项式的次数
;       interp<关键字参数: nearest>: 插值算法(包括: nearest, linear, cublic)
;       sub_percent<关键字参数: 0.1>: 默认使用10%的点位进行几何校正
;-
pro img_warp, target, lon, lat, out_res, target_warped, geo_info, degree=degree, interp=interp, $
    sub_percent=sub_percent

    ; 获取基本信息
    ds_size = size(target)
    lon_lat_corner = hash($  ; 考虑到min()max()为角点像元的中心处而非四个角点的边界
        'min_lon', min(lon) - out_res / 2.0d, $
        'max_lon', max(lon) + out_res / 2.0d, $
        'min_lat', min(lat) - out_res / 2.0d, $
        'max_lat', max(lat) + out_res / 2.0d)
    col_count_out = ceil((lon_lat_corner['max_lon'] - lon_lat_corner['min_lon']) / out_res)
    row_count_out = ceil((lon_lat_corner['max_lat'] - lon_lat_corner['min_lat']) / out_res)
    interp_code = hash($
        'nearest', 0, $
        'linear', 1, $
        'cublic', 2)
    if ~keyword_set(interp) then interp='nearest'
    if ~keyword_set(degree) then degree=5
    if ~keyword_set(sub_percent) then sub_percent = 0.1

    ; 原始的行列号矩阵
    row_ori = make_array(ds_size[1:2], /integer)
    col_ori = make_array(ds_size[1:2], /integer)
    for ix=0, ds_size[1]-1 do col_ori[ix, *] = ix
    for ix=0, ds_size[2]-1 do row_ori[*, ix] = ix

    ; 校正后的行列号矩阵
    col_warp = floor((lon - lon_lat_corner['min_lon']) / out_res)
    row_warp = floor((lon_lat_corner['max_lat'] - lat) / out_res)

    ; 获取sub_percent的均匀采样点索引
    all_count = ds_size[1] * ds_size[2]
    sample_count = floor(all_count * sub_percent)
    sample_ix = floor((findgen(sample_count) / double(sample_count)) * all_count)

    polywarp, col_ori[sample_ix], row_ori[sample_ix], col_warp[sample_ix], row_warp[sample_ix], $
        degree, k_col, k_row
    target_warped = poly_2d(target, k_col, k_row, interp_code[interp], $
        col_count_out, row_count_out, missing=!VALUES.F_nan)

    geo_info={$
        Modelpixelscaletag: [out_res, out_res, 0.0], $  ; 分辨率
        Modeltiepointtag: [0.0, 0.0, 0.0, lon_lat_corner['min_lon'], lon_lat_corner['max_lat'], 0.0], $  ; 角点信息
        Gtmodeltypegeokey: 2, $  ; 设置为地理坐标系
        Gtrastertypegeokey: 1, $  ; 像素的表示类型, 北上图像(North-Up)
        Geographictypegeokey: 4326, $  ; 地理坐标系为WGS84
        Geogcitationgeokey: 'GCS_WGS_1984', $
        Geogangularunitsgeokey: 9102}  ; 单位为度
end

但是上述是利用IDL从底层实现GLT校正，对于Python，我之前对FY卫星进行处理时使用与IDL类似的思路进行操作：

https://blog.csdn.net/m0_63001937/article/details/131626368?spm=1011.2415.3001.5331

但是gdal中有更封装好的方法针对提供二维的lon和lat，如何将指定变量var重投影为GeoTIFF文件，这是使用gdal实现的方法：

def img_glt(src_var, src_lon, src_lat, out_path, out_res, resample_alg=None, out_type=None, src_srs=None, dst_srs=None,
            out_bound=None):
    """
    基于二维的src_lon(纬度数据集)和src_lat（经度数据集）对src_var（目标数据集）进行GLT校正并输出为tiff文件
    :param src_var: 需要重投影/GLT校正的目标数据集, 要求三维np数组且shape=(波段数, 行数, 列数), 或二维数组且shape=(行数, 列数)
    :param src_lon: 经度数据集(二维np数组)
    :param src_lat: 纬度数据集(二维np数组)
    :param out_path: 重投影/GLT校正结果文件的输出路径(.tif)
    :param out_res: 输出分辨率(单位取决于dst_srs)
    :param resample_alg: 重采样算法(使用gdal.GRA_*), 默认为最近邻
    :param out_type: 输出像元值的数据类型, 默认依据src_var判断整型或者浮点
    :param src_srs: 输入栅格数据集的坐标系（lon和lat的坐标系,一般都是WGS84坐标系), 默认WGS84坐标系
    :param dst_srs: 输出栅格数据集的坐标系(重投影之后的坐标系), 默认WGS84坐标系
    :return: None
    """

    if resample_alg is None:
        resample_alg = gdal.GRA_NearestNeighbour

    if out_type is None:
        if np.issubdtype(src_var.dtype, np.integer):
            out_type = gdal.GDT_Int16
        elif np.issubdtype(src_var.dtype, np.floating):
            out_type = gdal.GDT_Float32
        else:
            out_type = gdal.GDT_Float32

    if len(src_var.shape) == 2:
        src_var = np.expand_dims(src_var, axis=0)  # (new_axis, rows, cols)

    # 定义输入栅格数据集的坐标系
    if src_srs is None:
        src_srs = osr.SpatialReference()
        src_srs.ImportFromEPSG(4326)
    # 定义输出栅格数据集的坐标系
    if dst_srs is None:
        dst_srs = osr.SpatialReference()
        dst_srs.ImportFromEPSG(4326)
        dst_srs.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)
    src_srs.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)  # 使用传统的(lon, lat)顺序, 下同
    dst_srs.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)
    # 获取基本信息

    band_count, rows, cols = src_var.shape

    # 为目标、经纬度分别创建临时tiff文件(内存中创建)
    var_mem_path = '/vsimem/var.tif'  # 内存中临时路径, 下同
    lon_mem_path = '/vsimem/lon.tif'
    lat_mem_path = '/vsimem/lat.tif'
    tiff_driver = gdal.GetDriverByName('GTiff')  # 创建TIFF驱动器
    img_var = tiff_driver.Create(var_mem_path, cols, rows, band_count, out_type)  # 创建tiff文件, 下同
    img_lon = tiff_driver.Create(lon_mem_path, cols, rows, 1, out_type)
    img_lat = tiff_driver.Create(lat_mem_path, cols, rows, 1, out_type)
    for band_ix in range(band_count):
        img_var.GetRasterBand(band_ix + 1).WriteArray(src_var[band_ix, :, :])  # 写入波段数据, lon和lat类似
    img_lon.GetRasterBand(1).WriteArray(src_lon)
    img_lat.GetRasterBand(1).WriteArray(src_lat)
    img_var = img_lon = img_lat = None  # 释放资源, 将缓冲区的数据全部写入

    # 更新元数据-关联地理定位
    img_var = gdal.Open(var_mem_path, gdal.GA_Update)
    """
    使用img_var=None后又使用gdal.Open打开似乎是重复冗余的操作, 实际并不是:
    首要原因是前面是在写入数据集, 而后面是在更新元数据 -- 写入和更新是不应该放在一起否则会混淆报错.
    这是因为写入之后本身就会生成元数据出来, 但是后面又有一个setMetadata更新GEOLOCATION会使得信息写入
    混乱最终导致后续重投影失败
    因此这是必要的.
    """
    img_var.SetMetadata({
        'SRS': src_srs.ExportToWkt(),  # X_DATASET和Y_DATASET所使用的坐标系
        'X_DATASET': lon_mem_path,  # 包含X坐标（通常是经度）的栅格数据集的路径
        'X_BAND': '1',  # 指定使用X坐标栅格数据集中的第一个波段
        'Y_DATASET': lat_mem_path,  # 包含Y坐标（通常是经度）的栅格数据集的路径
        'Y_BAND': '1',  # 指定使用Y坐标栅格数据集中的第一个波段
        'PIXEL_OFFSET': '0',  # X方向的偏移量
        'LINE_OFFSET': '0',  # Y方向的偏移量
        'PIXEL_STEP': '1',  # X方向的步长
        'LINE_STEP': '1'  # Y方向的步长
    }, 'GEOLOCATION')
    """
    对于传入字典中的最后四个变量：PIXEL_OFFSET、LINE_OFFSET、PIXEL_STEP和LINE_STEP。
    应当这样理解，它类似于于GLT校正,但是不同于ENVI IDL中的GLT校正工具是通过地面控制点，这里传入的每一个点都类似于控制点。
    这里传入X和Y的坐标类似于ENVI中的经纬度查找表，那我们知道，如果将所有的坐标点参与GLT校正那么在分辨率高地理覆盖范围大
    的情况下, 计算量会非常大运行时间会非常长，因此这里可以供你选择：
    假如从第200行300列开始从右隔4个像元采集一个像元坐标参与GLT校正，向下隔2个像元采集一个像元坐标参与GLT校正。
    那么应该设置为:
    'PIXEL_OFFSET': '300',  # X方向的偏移量
    'LINE_OFFSET': '200',  # Y方向的偏移量
    'PIXEL_STEP': '4',  # X方向的步长
    'LINE_STEP': '2'  # Y方向的步长

    这样的话大部分区域的像元都参与了运算而且计算量也下来了。
    """
    img_var.SetProjection(src_srs.ExportToWkt())  # 设置输入栅格数据集的坐标系
    img_var = None

    # 重投影/GLT校正
    warp_options = gdal.WarpOptions(
        dstSRS=dst_srs,
        outputBounds=out_bound,
        dstNodata=-999,
        xRes=out_res,
        yRes=out_res,
        resampleAlg=resample_alg,
        transformerOptions=['METHOD=GEOLOC_ARRAY'],
        targetAlignedPixels=True
    )
    gdal.Warp(out_path, var_mem_path, options=warp_options)

    # 释放内存文件
    gdal.Unlink(var_mem_path)
    gdal.Unlink(lon_mem_path)
    gdal.Unlink(lat_mem_path)

关于代码和投影相关的问题仍还有很多东西没有在博客中进行详细说明和探讨，但并不意味着他们不重要，如果你在这方面遇到了困惑，请联系我。

本文由博客一文多发平台 OpenWrite 发布！

flask, python, pytorch, tornado