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, 蘭伯特等角投影)投影座標系的GeoTIFF文件？ - 如何將
GLASS的LAI產品(hdf4格式,規則格網,全球)輸出為LCC座標系的GeoTIFF文件？ - 如何將
MODIS的土地覆蓋產品(MCD12C1, hdf4格式, 規則格網, 全球, 與GLASSLAI類似)輸出為LCC座標系的GeoTIFF文件？ - 如何將WGS84地理座標系的聚集指數GeoTIFF文件(全球, 規則格網, 與GLASS和MODIS類似但非hdf文件)輸出為LCC座標系的GeoTIFF文件？
- 如何將WGS84地理座標系的土壤類型GeoTIFF文件(全球, 規則格網, 與GLASS和MODIS類似但非hdf文件)輸出為LCC座標系的GeoTIFF文件？
要求：
- 所有結果包含兩種分辨率，9000m和3000m
- 相同分辨率的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 = None4.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產品進行過重投影：
- https://blog.csdn.net/m0_63001937/article/details/133977692?s...
- https://blog.csdn.net/m0_63001937/article/details/134238529?s...
最關鍵的部分即：
;+
; 函數用途:
; 基於經緯度數據集對目標數據集進行幾何校正(重投影-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?s...
但是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 發佈！
