|
| 1 | +#!/usr/bin/env python3 |
| 2 | +"""Convert an MITgcm binary pickup file to individual init files. |
| 3 | +
|
| 4 | +Reads the pickup.<iter>.data/.meta pair and writes: |
| 5 | + T.init.bin, S.init.bin, U.init.bin, V.init.bin, Eta.init.bin |
| 6 | +
|
| 7 | +The pickup is float64 (MITgcm default for checkpoints); init files are |
| 8 | +written as float32 (matching readBinaryPrec=32 in data PARM01). |
| 9 | +
|
| 10 | +Usage: |
| 11 | + python pickup_to_init.py <pickup_prefix> <output_dir> [--nx 768] [--ny 424] [--nr 50] |
| 12 | +
|
| 13 | +Example: |
| 14 | + python pickup_to_init.py repeat-year-50/001/pickup.0000087600 repeat-year-50/002/ |
| 15 | +""" |
| 16 | + |
| 17 | +import argparse |
| 18 | +import re |
| 19 | +import sys |
| 20 | +from pathlib import Path |
| 21 | + |
| 22 | +import numpy as np |
| 23 | + |
| 24 | + |
| 25 | +def parse_pickup_meta(meta_path: Path) -> dict: |
| 26 | + """Parse a MITgcm .meta file and return dims, precision, and field list.""" |
| 27 | + text = meta_path.read_text() |
| 28 | + |
| 29 | + # Extract dimensions |
| 30 | + dim_match = re.search(r"dimList\s*=\s*\[\s*([\d\s,]+)\]", text) |
| 31 | + if not dim_match: |
| 32 | + raise ValueError(f"Cannot parse dimList from {meta_path}") |
| 33 | + dims = [int(x) for x in dim_match.group(1).replace(",", " ").split()] |
| 34 | + nx, ny = dims[0], dims[3] |
| 35 | + |
| 36 | + # Extract precision |
| 37 | + prec_match = re.search(r"dataprec\s*=\s*\[\s*'(\w+)'\s*\]", text) |
| 38 | + dtype = np.float64 if prec_match and "64" in prec_match.group(1) else np.float32 |
| 39 | + |
| 40 | + # Extract number of records |
| 41 | + nrec_match = re.search(r"nrecords\s*=\s*\[\s*(\d+)\s*\]", text) |
| 42 | + nrecords = int(nrec_match.group(1)) if nrec_match else None |
| 43 | + |
| 44 | + # Extract field list |
| 45 | + fld_match = re.search(r"fldList\s*=\s*\{([^}]+)\}", text) |
| 46 | + if not fld_match: |
| 47 | + raise ValueError(f"Cannot parse fldList from {meta_path}") |
| 48 | + fields = re.findall(r"'(\w+)\s*'", fld_match.group(1)) |
| 49 | + |
| 50 | + return {"nx": nx, "ny": ny, "dtype": dtype, "nrecords": nrecords, "fields": fields} |
| 51 | + |
| 52 | + |
| 53 | +def pickup_to_init(pickup_prefix: str, output_dir: str, nx: int, ny: int, nr: int): |
| 54 | + """Read a pickup file and write individual init .bin files.""" |
| 55 | + meta_path = Path(pickup_prefix + ".meta") |
| 56 | + data_path = Path(pickup_prefix + ".data") |
| 57 | + out = Path(output_dir) |
| 58 | + |
| 59 | + if not meta_path.exists(): |
| 60 | + raise FileNotFoundError(f"Meta file not found: {meta_path}") |
| 61 | + if not data_path.exists(): |
| 62 | + raise FileNotFoundError(f"Data file not found: {data_path}") |
| 63 | + |
| 64 | + meta = parse_pickup_meta(meta_path) |
| 65 | + dtype = meta["dtype"] |
| 66 | + fields = meta["fields"] |
| 67 | + rec_size = nx * ny |
| 68 | + |
| 69 | + print(f"Pickup: {data_path}") |
| 70 | + print(f" Grid: {nx} x {ny} x {nr}") |
| 71 | + print(f" Precision: {dtype}") |
| 72 | + print(f" Fields: {fields}") |
| 73 | + print(f" Total records: {meta['nrecords']}") |
| 74 | + |
| 75 | + # Map pickup field names to init file names and their depth (nr for 3D, 1 for 2D) |
| 76 | + field_map = { |
| 77 | + "Uvel": ("U.init.bin", nr), |
| 78 | + "Vvel": ("V.init.bin", nr), |
| 79 | + "Theta": ("T.init.bin", nr), |
| 80 | + "Salt": ("S.init.bin", nr), |
| 81 | + "EtaN": ("Eta.init.bin", 1), |
| 82 | + } |
| 83 | + |
| 84 | + # Compute byte offsets for each field in the pickup |
| 85 | + bytes_per_val = np.dtype(dtype).itemsize |
| 86 | + rec_bytes = rec_size * bytes_per_val |
| 87 | + |
| 88 | + # Build offset table: walk through fields in order |
| 89 | + offsets = {} |
| 90 | + current_rec = 0 |
| 91 | + for fld in fields: |
| 92 | + # 3D fields have nr levels, 2D fields have 1 level |
| 93 | + if fld in ("EtaN", "dEtaHdt", "EtaH"): |
| 94 | + nlevels = 1 |
| 95 | + else: |
| 96 | + nlevels = nr |
| 97 | + offsets[fld] = (current_rec, nlevels) |
| 98 | + current_rec += nlevels |
| 99 | + |
| 100 | + print(f" Computed record layout: {offsets}") |
| 101 | + |
| 102 | + # Read and write the fields we need |
| 103 | + out.mkdir(parents=True, exist_ok=True) |
| 104 | + with open(data_path, "rb") as f: |
| 105 | + for fld_name, (init_name, nlevels) in field_map.items(): |
| 106 | + if fld_name not in offsets: |
| 107 | + print(f" WARNING: field '{fld_name}' not found in pickup, skipping") |
| 108 | + continue |
| 109 | + |
| 110 | + start_rec, expected_levels = offsets[fld_name] |
| 111 | + assert expected_levels == nlevels, ( |
| 112 | + f"Level mismatch for {fld_name}: expected {nlevels}, got {expected_levels}" |
| 113 | + ) |
| 114 | + |
| 115 | + # Seek to field start and read |
| 116 | + f.seek(start_rec * rec_bytes) |
| 117 | + data = np.fromfile(f, dtype=dtype, count=rec_size * nlevels) |
| 118 | + data = data.reshape((nlevels, ny, nx)) |
| 119 | + |
| 120 | + # Convert to float32 for init files |
| 121 | + init_path = out / init_name |
| 122 | + data.astype(np.float32).tofile(init_path) |
| 123 | + size_mb = init_path.stat().st_size / 1e6 |
| 124 | + print(f" Wrote {init_path} ({size_mb:.1f} MB)") |
| 125 | + |
| 126 | + print("Done.") |
| 127 | + |
| 128 | + |
| 129 | +def main(): |
| 130 | + parser = argparse.ArgumentParser(description="Convert MITgcm pickup to init files") |
| 131 | + parser.add_argument("pickup_prefix", help="Path prefix (without .data/.meta)") |
| 132 | + parser.add_argument("output_dir", help="Directory to write init files") |
| 133 | + parser.add_argument("--nx", type=int, default=768, help="Grid points in X") |
| 134 | + parser.add_argument("--ny", type=int, default=424, help="Grid points in Y") |
| 135 | + parser.add_argument("--nr", type=int, default=50, help="Grid points in Z") |
| 136 | + args = parser.parse_args() |
| 137 | + |
| 138 | + pickup_to_init(args.pickup_prefix, args.output_dir, args.nx, args.ny, args.nr) |
| 139 | + |
| 140 | + |
| 141 | +if __name__ == "__main__": |
| 142 | + main() |
0 commit comments