求问超算里grads画站点降水分布图（插值、非散点图)

wujiayi3946

将站点数据转变的python代码：（站点数据是000格式）#!/usr/bin/env python3
import struct
import sys

def micaps3_to_grads(input_file, output_prefix):
    with open(input_file, 'r', encoding='utf-8', errors='ignore') as f:
        lines = f.readlines()

    lines = [line.rstrip('\n') for line in lines]

    header2 = lines[1].strip().split()
    year, month, day, hour = int(header2[0]), int(header2[1]), int(header2[2]), int(header2[3])
    full_year = 2000 + year if year < 100 else year

    print(f"时间: {full_year}年{month}月{day}日{hour}时")

    # 找站点数据行
    station_start_idx = None
    nstn_expected = 0
    for i in range(2, len(lines)):
        parts = lines[i].strip().split()
        if len(parts) == 2 and parts[0] == '1' and parts[1].isdigit():
            station_start_idx = i
            nstn_expected = int(parts[1])
            print(f"文件声明站点数: {nstn_expected}")
            break

    if station_start_idx is None:
        raise ValueError("找不到站点标记行")

    data_start = station_start_idx + 1

    # 严格解析站点数据，直到文件结束或空行
    station_data = []
    for i in range(data_start, len(lines)):
        line = lines[i].strip()
        if not line:
            continue

        parts = line.split()
        if len(parts) < 5:
            continue

        try:
            stid = parts[0]
            lon = float(parts[1])
            lat = float(parts[2])
            elev = float(parts[3])  # 海拔，不写入但检查格式
            rain = float(parts[4])

            # 确保站号严格8字符
            stid_padded = stid.ljust(8)[:8]

            station_data.append({
                'stid': stid_padded,
                'lat': lat,
                'lon': lon,
                'rain': rain
            })
        except Exception as e:
            print(f"警告: 行{i+1}解析失败: {line[:60]}... 错误: {e}")
            continue

    nstn_actual = len(station_data)
    print(f"实际解析站点数: {nstn_actual}")

    if nstn_actual != nstn_expected:
        print(f"警告: 实际站点数({nstn_actual})与声明({nstn_expected})不符，以实际为准")

    # 写二进制 - 严格32字节/记录 (28头+4数据)
    dat_file = f"{output_prefix}.dat"
    with open(dat_file, 'wb') as f:
        for stn in station_data:
            # 站号8字节ASCII
            stid_bytes = stn['stid'].encode('ascii')

            # 头部: 8s + f + f + f + i + i = 28字节 (小端)
            header = struct.pack('<8s f f f i i',
                stid_bytes,      # 8 bytes
                stn['lat'],      # 4 bytes (float32)
                stn['lon'],      # 4 bytes
                0.0,             # 4 bytes (time)
                1,               # 4 bytes (nlev, int32)
                0                # 4 bytes (flag, int32)
            )

            # 数据: 4字节float
            data = struct.pack('<f', stn['rain'])

            # 写入完整记录
            f.write(header)
            f.write(data)

            # 验证记录大小
            if len(header) != 28 or len(data) != 4:
                print(f"错误: 记录大小异常 header={len(header)}, data={len(data)}")

    # 验证文件大小
    import os
    file_size = os.path.getsize(dat_file)
    expected_size = nstn_actual * 32
    print(f"文件大小: {file_size} 字节 (期望: {expected_size} 字节)")

    if file_size != expected_size:
        print(f"警告: 文件大小不匹配！")

    # 写ctl - 使用实际站点数
    ctl_file = f"{output_prefix}.ctl"
    mon_names = ['jan','feb','mar','apr','may','jun','jul','aug','sep','oct','nov','dec']
    mon_str = mon_names[month-1]

    with open(ctl_file, 'w') as f:
        f.write(f"""dset ./{dat_file}
dtype station
stnmap ./{output_prefix}.map
undef -999.0
title Rainfall_{full_year}{month:02d}{day:02d}{hour:02d}

tdef 1 linear {hour:02d}Z{day:02d}{mon_str}{full_year} 1hr

vars 1
p 0 99 Precipitation(mm)
endvars
""")

    print(f"生成ctl: {ctl_file}")
    print(f"准备运行: stnmap -i {ctl_file}")

if __name__ == "__main__":
    if len(sys.argv) != 3:
        print("用法: python3 convert_fix.py <输入> <输出前缀>")
        sys.exit(1)
    micaps3_to_grads(sys.argv[1], sys.argv[2])


用上面生成了.dat和.ctl
但是stnmap -i rain26011817.ctl时，有以下错误：
  Name of binary data set: rain26011817.dat
  Number of times in the data set: 1
  Number of surface variables: 1
  Number of level dependent variables: 0

Starting scan of station data binary file.
Binary data file open: rain26011817.dat

Processing time step 1
  Low Level I/O Error:  Read error on data file
  Error reading 28 bytes at location 38560
  Possible cause: Premature EOF
求问问题在哪里？grads2.2.1版本，在超算里