如何用Python的NetCDF4/h5py写入GeoTraj类型数据

GeoTraj本质是一维轨迹序列数据，核心是绑定到轨迹点索引的经度、纬度、时间序列，和你提供的网格点写入逻辑核心区别在于维度设计。以下是两种库的实现方案：

一、NetCDF4库实现 #

import netCDF4 as nc
import numpy as np
from datetime import datetime

# 创建NetCDF文件
with nc.Dataset('geotraj_example.nc', 'w', format='NETCDF4') as f:
    # 定义核心维度：轨迹点数量（示例设为10个点）
    f.createDimension('trajectory_point', 10)
    
    # 写入经度变量及属性
    lon = f.createVariable('lon', np.float64, ('trajectory_point',))
    lon.standard_name = 'longitude'
    lon.long_name = 'Longitude of trajectory point'
    lon.units = 'degrees_east'
    lon[:] = np.linspace(110.0, 120.0, 10)  # 模拟东向移动轨迹
    
    # 写入纬度变量及属性
    lat = f.createVariable('lat', np.float64, ('trajectory_point',))
    lat.standard_name = 'latitude'
    lat.long_name = 'Latitude of trajectory point'
    lat.units = 'degrees_north'
    lat[:] = np.linspace(30.0, 35.0, 10)  # 模拟北向移动轨迹
    
    # 写入时间变量（GeoTraj必填）
    time = f.createVariable('time', np.int64, ('trajectory_point',))
    time.standard_name = 'time'
    time.long_name = 'Time of trajectory observation'
    time.units = 'seconds since 1970-01-01 00:00:00'
    # 生成模拟时间序列：2024年1月1日起每小时一个点
    start_time = datetime(2024, 1, 1)
    time_vals = np.array([(start_time + np.timedelta64(i, 'h')).timestamp() for i in range(10)], dtype=np.int64)
    time[:] = time_vals
    
    # 添加GeoTraj规范的全局元数据
    f.title = 'Example GeoTrajectory Dataset'
    f.trajectory_id = 'TRAJ_001'
    f.geospatial_lat_min = lat[:].min()
    f.geospatial_lat_max = lat[:].max()
    f.geospatial_lon_min = lon[:].min()
    f.geospatial_lon_max = lon[:].max()

二、h5py库实现 #

h5py需手动管理数据集和属性，逻辑和NetCDF4一致：

import h5py
import numpy as np
from datetime import datetime

# 创建HDF5文件
with h5py.File('geotraj_example.h5', 'w') as f:
    # 创建轨迹点维度
    f.create_dataset('trajectory_point', data=np.arange(10))
    
    # 写入经度数据集及属性
    lon_ds = f.create_dataset('lon', data=np.linspace(110.0, 120.0, 10), dtype='f8')
    lon_ds.attrs['standard_name'] = 'longitude'
    lon_ds.attrs['long_name'] = 'Longitude of trajectory point'
    lon_ds.attrs['units'] = 'degrees_east'
    
    # 写入纬度数据集及属性
    lat_ds = f.create_dataset('lat', data=np.linspace(30.0, 35.0, 10), dtype='f8')
    lat_ds.attrs['standard_name'] = 'latitude'
    lat_ds.attrs['long_name'] = 'Latitude of trajectory point'
    lat_ds.attrs['units'] = 'degrees_north'
    
    # 写入时间数据集及属性
    start_time = datetime(2024, 1, 1)
    time_vals = np.array([(start_time + np.timedelta64(i, 'h')).timestamp() for i in range(10)], dtype='i8')
    time_ds = f.create_dataset('time', data=time_vals)
    time_ds.attrs['standard_name'] = 'time'
    time_ds.attrs['long_name'] = 'Time of trajectory observation'
    time_ds.attrs['units'] = 'seconds since 1970-01-01 00:00:00'
    
    # 添加全局元数据
    f.attrs['title'] = 'Example GeoTrajectory Dataset'
    f.attrs['trajectory_id'] = 'TRAJ_001'
    f.attrs['geospatial_lat_min'] = lat_ds[:].min()
    f.attrs['geospatial_lat_max'] = lat_ds[:].max()
    f.attrs['geospatial_lon_min'] = lon_ds[:].min()
    f.attrs['geospatial_lon_max'] = lon_ds[:].max()

核心注意点 #

• 所有变量绑定到trajectory_point一维维度，而非网格的二维/多维维度
• 必须包含lon、lat、time三个核心变量，且需附上符合CF规范的单位和属性
• 全局属性需补充轨迹唯一ID、空间范围等元数据，满足GeoTraj格式要求

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