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

wordpress4.3.1宁波正规优化seo公司

news 2025/9/26 6:44:01
wordpress4.3.1,宁波正规优化seo公司,文职培训机构前十名,做静态网站的参考文献01 问题和说明 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, 兰伯特等角投影)投影坐…

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/405689.html

相关文章：

  • 网站建设是设计师吗响应式网站布局
  • 做游戏开箱网站的法律风险个人网站怎样申请icp
  • 做视频特效的网站湖南常德论坛
  • 网站分析报告范文2000网站建设项目中标通知
  • 公司网站是否有必要销售产品网站开发类合同范本
  • 网站开发语言选择百度收录哪些网站吗
  • 广州建站方法wordpress搭建多个购物网站
  • 营销网站费用你理解的网络营销是什么
  • 做网站项目需要多少钱网站经营性备案流程
  • 域名备案 没有网站吗环境设计案例网站
  • 如何创建一个免费的网站免费做字体的网站好
  • 网站建设实战集团网站建设多少钱
  • wordpress 网站播放器插件珠海移动网站设计
  • 网站建设需要的技术wordpress数据表文档
  • 如何做网站页面免费的服装网站建设的规模和类别
  • 怎样才能做网站百度竞价排名系统
  • 东莞做网站有哪些北京网站设计精选刻
  • 自助建设网站平台中山市企业网站seo哪家好
  • 章丘做网站优化阿里云服务器免费一年
  • 玄圭做网站怎么样网站建设招聘信息
  • 网站开发新闻管理系统的背景河北做it的网站
  • 网站vpsseo是做什么的
  • 宁波网站推广排名永州市规划建设局网站
  • 网站上面的彩票快3怎么做下载的asp网站怎么打开
  • 网站信任 用户转化html指令代码大全
  • 云谷系统网站开发seo的概念是什么
  • 网站开发工具报告网站空间控制面板
  • 教师在哪些网站可以做兼职wordpress移除仪表盘
  • 青岛建设监理协会网站亚马逊雨林地图
  • 英语可以做推广的亲子类网站seo专员是什么职业