当前位置： 首页 > news >正文

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

news 2025/8/15 11:46:01

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文件形如：

T2

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

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

1.2.2 GLASS-LAI产品

LAI产品形如：
GLASS-LAI产品

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

GLASS-LAI的元数据

1.2.3 MCD12C1-土地覆盖产品

MCD12C1产品形如：

MCD12C1产品

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

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

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

聚集指数和土壤类型的基本信息1
聚集指数和土壤类型的基本信息2

02 思路

思路1：

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

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

思路2：

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

LCC投影参数
其中：

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
LCC投影转换示意图

如果你不想去研究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_NearestNeighbourif out_type is None:if np.issubdtype(src_var.dtype, np.integer):out_type = gdal.GDT_Int16elif np.issubdtype(src_var.dtype, np.floating):out_type = gdal.GDT_Float32else:out_type = gdal.GDT_Float32if 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_Float32if resamlpe_alg is None:resamlpe_alg = gdal.GRA_NearestNeighbourrows, cols = src_var.shapevar_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_resul_y_center = cen_x + (rows - 1) / 2 * y_res# 计算左上角像元(第一行第一列)的左上角角点坐标ul_x_corner = ul_x_center - x_res / 2ul_y_corner = ul_y_center + y_res / 2# 构建GeoTransformgeo_transform = [ul_x_corner, x_res, 0, ul_y_corner, 0, -y_res]return geo_transformdef 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.shapedriver = 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.RasterXSizeheight = ref_ds.RasterYSizeref_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 = Nonewarp_options = gdal.WarpOptions(format='GTiff',dstSRS=ref_srs_wkt,  # 目标投影, LCCxRes=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.RasterXSizeheight = ref_ds.RasterYSizeref_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 = Nonewarp_options = gdal.WarpOptions(format='GTiff',dstSRS=ref_srs_wkt,  # 目标投影, LCCxRes=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 CRSgdal.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.RasterXSizeheight = ref_ds.RasterYSizeref_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.nanscaled_data = original_data * 0.01# 创建一个新的栅格文件用于写入结果var_mem_path = '/vsimem/var.tif'driver = gdal.GetDriverByName('GTiff')rows, cols = original_data.shapemem_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,  # 目标投影, LCCxRes=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.RasterXSizeheight = ref_ds.RasterYSizeref_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.shapemem_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,  # 目标投影, LCCxRes=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=5if ~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, *] = ixfor 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_rowtarget_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, $  ; 地理坐标系为WGS84Geogcitationgeokey: '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_NearestNeighbourif out_type is None:if np.issubdtype(src_var.dtype, np.integer):out_type = gdal.GDT_Int16elif np.issubdtype(src_var.dtype, np.floating):out_type = gdal.GDT_Float32else:out_type = gdal.GDT_Float32if 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)

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

http://www.dtcms.com/a/331793.html

相关文章：

  • 深度学习 --- ResNet神经网络
  • 【递归完全搜索】CCC 2008 - 24点游戏Twenty-four
  • 【完整源码+数据集+部署教程】膝关节屈伸运动检测系统源码和数据集：改进yolo11-RFAConv
  • pip和dnf只下载不安装离线包
  • 沈帅波出席茅台红缨子高粱节探讨产业赋能新模式
  • Ansys FreeFlow入门：对搅拌罐进行建模
  • 【159页PPT】机械制造行业数字化转型某著名企业U8系统全解决方案（附下载方式）
  • Avalonia_SukiUI明暗主题切换时部分元素颜色不变
  • jetson orin nx(8G)烧录super系统实录
  • Ubuntu下载、安装、编译指定版本python
  • 机器学习--KNN算法
  • Linux入门指南：基础开发工具---yum/apt
  • 单北斗GNSS变形监测应用解析
  • 读《精益数据分析》：移情（Empathy）—— 验证真实需求，避免伪需求陷阱
  • 大模型工程化落地：从模型选择到性能优化的实战指南
  • C#笔记啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊啊
  • 机器学习学习报告
  • 【博客系统测试报告】---接口自动化测试
  • AI幻觉终结之后：GPT-5开启的“可靠性”新赛道与开发者生存指南
  • JAVA中正则表达式详解
  • 前端八股文-CSS3篇
  • 考研408《计算机组成原理》复习笔记，第四章(2)——指令寻址和数据寻址
  • K8s-kubernetes(二)资源限制-详细介绍
  • 2025 年电赛 C 题 发挥部分 1：多正方形 / 重叠正方形高精度识别与最小边长测量
  • 悲观锁乐观锁与事务注解在项目实战中的应用场景及详细解析
  • 如何解决EMI中传导干扰
  • Spring-解决项目依赖异常问题
  • 【从零开始java学习|第六篇】运算符的使用与注意事项
  • 因果推断在用户流失预警的案例研究
  • 第2节：多模态的核心问题（多模态大模型基础教程）