我经常处理这类海洋遥感HDF数据，给你一套分步可行的方案，分命令行/GDAL+Python和ArcGIS原生工具两种路径，你可以根据自己的熟悉程度选择：

一、命令行+GDAL+Python方案（灵活高效，适合批量处理） #

1. 准备工具 #

先安装GDAL（包含处理HDF4的驱动）：

• Windows：通过OSGeo4W安装，选择GDAL相关包
• Mac：用Homebrew执行 brew install gdal
• Linux：用系统包管理器（比如Ubuntu执行 sudo apt install gdal-bin python3-gdal）

2. 解压.tar文件 #

用命令行批量解压所有年度压缩包：

for tar_file in *.tar; do
  tar -xvf "$tar_file"
done

3. 提取HDF4中的目标数据集 #

首先用gdalinfo查看HDF4文件里的子数据集（找到初级生产力对应的那一项）：

gdalinfo your_file.hdf

输出里会有类似SUBDATASET_1_NAME=HDF4_SDS:UNKNOWN:"your_file.hdf":0的行，这就是你要提取的数据集路径。

然后批量提取所有HDF4为GeoTIFF：

for hdf_file in *.hdf; do
  # 替换下面的SUBDATASET路径为你实际查到的内容
  gdal_translate HDF4_SDS:UNKNOWN:"$hdf_file":0 "${hdf_file%.hdf}_pp_raw.tif"
done

4. 计算单位面积浓度 #

如果原始数据是栅格总碳量（而非单位面积），需要除以每个栅格的实际面积（经纬度投影下，栅格面积随纬度变化）：

步骤4.1：生成栅格面积矩阵（用Python脚本）

import rasterio
import numpy as np

# 读取第一个提取的TIFF获取地理信息
with rasterio.open('first_pp_raw.tif') as src:
    meta = src.meta
    height, width = src.shape
    transform = src.transform
    # 获取所有纬度值（假设左上角是北纬90°）
    latitudes = np.arange(height) * transform[4] + transform[5]
    pixel_width_deg = transform[0]
    pixel_height_deg = abs(transform[4])
    earth_radius = 6371000  # 地球半径（米）

# 计算每个栅格的面积（平方米）
lat_rad = np.deg2rad(latitudes)
pixel_width_m = 2 * np.pi * earth_radius * np.cos(lat_rad) * pixel_width_deg / 360
pixel_height_m = 2 * np.pi * earth_radius * pixel_height_deg / 360
area_matrix = np.tile(pixel_width_m[:, np.newaxis], (1, width)) * pixel_height_m

# 保存面积栅格
meta.update(dtype=rasterio.float32)
with rasterio.open('pixel_area.tif', 'w', **meta) as dst:
    dst.write(area_matrix.astype(rasterio.float32), 1)

步骤4.2：批量计算单位面积浓度

for raw_tif in *_pp_raw.tif; do
  gdal_calc.py -A "$raw_tif" -B pixel_area.tif --outfile="${raw_tif%_pp_raw.tif}_pp_per_m2.tif" --calc="A/B" --NoDataValue=-9999
done

5. 计算多年空间平均值 #

用Python批量读取所有处理后的栅格，计算每个像素的时间序列均值：

import rasterio
import numpy as np
import glob

# 获取所有单位面积浓度栅格
pp_files = glob.glob('*_pp_per_m2.tif')

# 读取元数据并初始化数组
with rasterio.open(pp_files[0]) as src:
    meta = src.meta
    height, width = src.shape
    all_data = np.zeros((len(pp_files), height, width), dtype=np.float32)

# 读取所有栅格数据
for i, file in enumerate(pp_files):
    with rasterio.open(file) as src:
        all_data[i] = src.read(1)

# 计算均值（忽略缺失值）
mean_data = np.nanmean(all_data, axis=0)

# 保存多年平均栅格
meta.update(dtype=rasterio.float32, nodata=-9999)
with rasterio.open('mean_pp_per_m2.tif', 'w', **meta) as dst:
    dst.write(mean_data.astype(rasterio.float32), 1)

6. 转换为Shapefile（可选） #

用GDAL把平均栅格转成矢量Shapefile：

gdal_polygonize.py mean_pp_per_m2.tif -f "ESRI Shapefile" mean_pp.shp mean_pp productivity_value

7. 确保ArcGIS兼容性 #

如果栅格没有地理参考，手动添加WGS84投影（EPSG:4326）：

gdal_edit.py -a_srs EPSG:4326 mean_pp_per_m2.tif

二、ArcGIS原生工具方案（适合熟悉ArcGIS的用户） #

如果你不想用代码，用ArcGIS的工具链也能完成：

1. 提取HDF4数据集：打开HDF To Raster工具，选择HDF文件，在"Subdataset"下拉菜单中找到初级生产力数据集，输出栅格。用ModelBuilder或者ArcPy脚本批量处理所有28个文件。
2. 计算单位面积浓度：如果需要，用Raster Calculator工具，结合栅格面积（可以用Create Constant Raster+公式计算，或者用上面生成的pixel_area.tif）进行除法运算。
3. 计算多年平均值：打开Cell Statistics工具，选择所有处理后的栅格，统计类型选Mean，输出多年平均栅格。
4. 转换为Shapefile：用Raster to Polygon工具，把平均栅格转成Shapefile，注意勾选"Simplify polygons"优化结果。

常见问题排查 #

• 转换失败可能是因为没选对HDF4的子数据集：一定要用gdalinfo或者ArcGIS的HDF预览确认目标数据集路径。
• 地理参考丢失：如果输出栅格没有投影，手动指定EPSG:4326（全球海洋数据默认是经纬度投影）。
• 缺失值处理：计算均值时一定要忽略NaN或无效值，否则结果会失真。

内容的提问来源于stack exchange，提问作者Kristina