#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
FY3D/MERSI 哈默投影 NDVI/EVI -> EPSG:4326
修复：Y 轴方向翻转导致的纬度整体偏低问题
"""
import os
import h5py
import numpy as np
from pyproj import CRS
import rasterio
from rasterio.transform import from_bounds
from rasterio.warp import calculate_default_transform, reproject, Resampling
# ---------- 用户配置 ----------
H5_FILE_PATH = "FY3D_MERSI_0000_L3_NVI_MLT_HAM_20251130_AOTD_0250M_MS.HDF"
OUTPUT_TIF_PATH = "FY3D_NDVI_4326_fixed.tif"
DATASET_NAME = "250M_10day_NDVI" # 或 EVI
OUTPUT_NODATA = -9999.0
EARTH_RADIUS_KM = 6363.961
EARTH_RADIUS_M = EARTH_RADIUS_KM * 1000.0
# --------------------------------
def get_fy3d_meta_params(h5_file):
attrs = h5_file.attrs
def get(k, d=0.0): return attrs.get(k, d)
return dict(
center_lon=float(get("Projection_Center_Longitude", 0.0)),
center_lat=float(get("Projection_Center_Latitude", 0.0)),
left_top_x=float(get("Left-Top_X", 0.0)),
left_top_y=float(get("Left-Top_Y", 0.0)),
resolution_x=float(get("Resolution_X", 0.25)),
resolution_y=float(get("Resolution_Y", 0.25)),
data_lines=int(get("Data_Lines", 4000)),
data_pixels=int(get("Data_Pixels", 4000)),
fill_value=int(get("FillValue", -32768)),
slope=float(get("Slope", 1.0e-4)),
intercept=float(get("Intercept", 0.0)),
)
def build_xy_coords_km(pp):
"""返回像素中心坐标（km）"""
W, H = pp["data_pixels"], pp["data_lines"]
dx, dy = pp["resolution_x"], pp["resolution_y"]
x0, y0 = pp["left_top_x"], pp["left_top_y"]
x = x0 + (np.arange(W) + 0.5) * dx
y = y0 + (np.arange(H) + 0.5) * dy # ✅ 关键修复：Y 轴向上
return x, y
def write_reprojected_from_hammer(ndvi_array, x_km, y_km, out_tif,
center_lon=0.0, nodata=OUTPUT_NODATA,
R=EARTH_RADIUS_M):
H, W = ndvi_array.shape
assert W == x_km.size and H == y_km.size
dx_km, dy_km = x_km[1] - x_km[0], y_km[1] - y_km[0]
hx, hy = dx_km/2, dy_km/2
xmin_km = float(x_km.min() - hx)
xmax_km = float(x_km.max() + hx)
ymin_km = float(y_km.min() - hy)
ymax_km = float(y_km.max() + hy)
# 转米
xmin_m, xmax_m = xmin_km*1000, xmax_km*1000
ymin_m, ymax_m = ymin_km*1000, ymax_km*1000
src_trans = from_bounds(xmin_m, ymin_m, xmax_m, ymax_m, W, H)
src_crs = CRS.from_proj4(f"+proj=hammer +lon_0={center_lon} +R={int(R)} +units=m +no_defs")
dst_crs = CRS.from_epsg(4326)
dst_trans, dst_w, dst_h = calculate_default_transform(
src_crs, dst_crs, W, H, left=xmin_m, bottom=ymin_m, right=xmax_m, top=ymax_m)
dst = np.full((dst_h, dst_w), nodata, np.float32)
src = np.where(np.isnan(ndvi_array), nodata, ndvi_array).astype(np.float32)
reproject(source=src, destination=dst,
src_transform=src_trans, src_crs=src_crs, src_nodata=nodata,
dst_transform=dst_trans, dst_crs=dst_crs, dst_nodata=nodata,
resampling=Resampling.bilinear)
meta = dict(driver="GTiff", height=dst_h, width=dst_w, count=1,
dtype=dst.dtype, crs=dst_crs, transform=dst_trans,
nodata=nodata, compress="LZW", tiled=True)
with rasterio.open(out_tif, "w", **meta) as f:
f.write(dst, 1)
print("写出完成 →", out_tif)
ul = dst_trans * (0, 0)
lr = dst_trans * (dst_w, dst_h)
print("输出左上 (lon lat):", ul)
print("输出右下 (lon lat):", lr)
def main(h5_path, ds_name, out_tif):
if not os.path.isfile(h5_path):
print("文件不存在:", h5_path)
return
with h5py.File(h5_path, "r") as h5:
print("HDF keys:", list(h5.keys()))
if ds_name not in h5:
print("数据集不存在:", ds_name); return
raw = h5[ds_name][()]
pp = get_fy3d_meta_params(h5)
# 清洗 NDVI
fill, slope, off = pp["fill_value"], pp["slope"], pp["intercept"]
ndvi = raw.astype(np.float32)
ndvi = np.where(ndvi == fill, np.nan, ndvi) * slope + off
ndvi = np.clip(ndvi, -1.0, 1.0)
print("NDVI 范围:", np.nanmin(ndvi), "~", np.nanmax(ndvi))
x_km, y_km = build_xy_coords_km(pp)
print("X km 范围:", x_km.min(), "~", x_km.max())
print("Y km 范围:", y_km.min(), "~", y_km.max())
write_reprojected_from_hammer(ndvi, x_km, y_km, out_tif,
center_lon=pp["center_lon"])
if __name__ == "__main__":
main(H5_FILE_PATH, DATASET_NAME, OUTPUT_TIF_PATH)