This commit is contained in:
Anatoly Antonov 2026-05-18 02:09:07 +09:00
parent c4f355a32e
commit 7a8d5d13fa
72 changed files with 4510 additions and 4104 deletions

155
scripts/build_elevation.py Normal file
View file

@ -0,0 +1,155 @@
#!/usr/bin/env python3
"""
Download ETOPO 2022 30-arc-second elevation data and convert to ruaumoko-compatible
binary format (int16 little-endian, 21601 lat x 43200 lon, south-to-north).
Output: ~1.74 GiB binary file.
Usage: python3 build_elevation.py [output_path]
Default output: /srv/ruaumoko-dataset
"""
import sys
import os
import struct
import tempfile
import numpy as np
CELLS_PER_DEGREE = 120
NUM_LATS = 180 * CELLS_PER_DEGREE + 1 # 21601
NUM_LONS = 360 * CELLS_PER_DEGREE # 43200
EXPECTED_SIZE = NUM_LATS * NUM_LONS * 2 # 1,866,326,400
ETOPO_URL = "https://www.ngdc.noaa.gov/thredds/fileServer/global/ETOPO2022/30s/30s_surface_elev_netcdf/ETOPO_2022_v1_30s_N90W180_surface.nc"
def download_etopo(output_path):
"""Download ETOPO 2022 NetCDF and convert to ruaumoko binary format."""
try:
import xarray as xr
except ImportError:
print("ERROR: xarray is required. Install with: pip install xarray netcdf4")
sys.exit(1)
# Check if we can download directly or need a local file
nc_path = os.environ.get("ETOPO_NC_PATH")
if nc_path and os.path.exists(nc_path):
print(f"Using local ETOPO file: {nc_path}")
else:
print(f"Downloading ETOPO 2022 30-second data (~1.1 GB)...")
print(f" URL: {ETOPO_URL}")
print(f" (Set ETOPO_NC_PATH env var to use a pre-downloaded file)")
import urllib.request
with tempfile.NamedTemporaryFile(suffix=".nc", delete=False) as f:
nc_path = f.name
try:
urllib.request.urlretrieve(ETOPO_URL, nc_path, _progress)
print()
except Exception as e:
os.unlink(nc_path)
print(f"\nDownload failed: {e}")
print("\nAlternative: manually download ETOPO 2022 30s NetCDF from:")
print(" https://www.ncei.noaa.gov/products/etopo-global-relief-model")
print(f" Then set ETOPO_NC_PATH=/path/to/file.nc and re-run")
sys.exit(1)
print(f"Opening NetCDF dataset...")
ds = xr.open_dataset(nc_path)
# ETOPO 2022 30s has:
# - lat: -90 to +90, 21601 points (south to north)
# - lon: -180 to +180, 43201 points
# We need:
# - lat: -90 to +90, 21601 points (south to north) ← same
# - lon: 0 to 360 (exclusive), 43200 points ← need to shift and drop last
z = ds["z"] # elevation variable
print(f" Shape: {z.shape}")
print(f" Lat range: {float(z.lat.min())} to {float(z.lat.max())}")
print(f" Lon range: {float(z.lon.min())} to {float(z.lon.max())}")
# Sort latitude south-to-north (should already be, but ensure)
z = z.sortby("lat")
# Shift longitude from [-180, 180] to [0, 360)
print("Shifting longitude to [0, 360)...")
z = z.assign_coords(lon=(z.lon % 360))
z = z.sortby("lon")
data = z.values
print(f" Raw shape after sort: {data.shape}")
# Handle longitude dimension: drop last col if it wraps (43201 → 43200)
if data.shape[1] == NUM_LONS + 1:
data = data[:, :NUM_LONS]
elif data.shape[1] != NUM_LONS:
print(f"ERROR: unexpected lon dimension: {data.shape[1]}, expected {NUM_LONS} or {NUM_LONS+1}")
sys.exit(1)
# Handle latitude dimension: ETOPO 2022 is cell-centered (21600 rows),
# ruaumoko expects grid-registered (21601 rows including both poles).
# Pad by duplicating edge rows for the poles.
if data.shape[0] == NUM_LATS - 1:
print(f" Padding latitude from {data.shape[0]} to {NUM_LATS} (adding north pole row)")
north_pole = data[-1:, :] # duplicate +89.99... as +90
data = np.concatenate([data, north_pole], axis=0)
elif data.shape[0] != NUM_LATS:
print(f"ERROR: unexpected lat dimension: {data.shape[0]}, expected {NUM_LATS} or {NUM_LATS-1}")
sys.exit(1)
print(f"Final grid shape: {data.shape}")
print(f"Elevation range: {data.min():.1f} to {data.max():.1f} metres")
# Write as int16 little-endian
print(f"Writing to {output_path}...")
elev_int16 = np.clip(data, -32768, 32767).astype(np.dtype("<i2"))
elev_int16.tofile(output_path)
actual_size = os.path.getsize(output_path)
print(f"Written {actual_size:,} bytes (expected {EXPECTED_SIZE:,})")
if actual_size == EXPECTED_SIZE:
print("SUCCESS")
else:
print("WARNING: size mismatch!")
ds.close()
# Spot check
verify(output_path)
def verify(path):
"""Quick spot-check of the elevation dataset."""
data = np.memmap(path, dtype="<i2", mode="r", shape=(NUM_LATS, NUM_LONS))
tests = [
("Mt Everest (~28.0N, 86.9E)", 28.0, 86.9, 8000, 9000),
("Dead Sea (~31.5N, 35.5E)", 31.5, 35.5, -500, 0),
("Pacific Ocean (~0N, 180E)", 0.0, 180.0, -6000, 0),
("Auburn AU (~-34.0S, 138.7E)", -34.03, 138.69, 200, 400),
]
print("\n Spot-check:")
for name, lat, lon, lo, hi in tests:
lat_idx = int((lat + 90) * CELLS_PER_DEGREE)
lon_idx = int(lon * CELLS_PER_DEGREE)
val = int(data[lat_idx, lon_idx])
ok = "OK" if lo <= val <= hi else "FAIL"
print(f" {name}: {val}m [{ok}] (expected {lo}-{hi})")
_last_pct = -1
def _progress(block_num, block_size, total_size):
global _last_pct
if total_size > 0:
pct = int(block_num * block_size * 100 / total_size)
if pct != _last_pct and pct % 5 == 0:
_last_pct = pct
print(f" {pct}%...", end="", flush=True)
if __name__ == "__main__":
output = sys.argv[1] if len(sys.argv) > 1 else "/srv/ruaumoko-dataset"
download_etopo(output)

View file

@ -1,303 +0,0 @@
#!/usr/bin/env python3
import subprocess
import sys
import time
import requests
import json
from typing import Any
import base64
import math
# --- Config ---
REFERENCE_API_URL = "https://fly.stratonautica.ru/api/v2/?profile=standard_profile&pred_type=single&launch_datetime=2025-06-25T20%3A45%3A00Z&launch_latitude=56.6992&launch_longitude=38.8247&launch_altitude=0&ascent_rate=5&burst_altitude=30000&descent_rate=5"
LOCAL_API_URL = "http://localhost:8080/api/v1/prediction?profile=standard_profile&pred_type=single&launch_datetime=2025-06-25T20%3A45%3A00Z&launch_latitude=56.6992&launch_longitude=38.8247&launch_altitude=0&ascent_rate=5&burst_altitude=30000&descent_rate=5"
LOCAL_API_PAYLOAD = {
"launch_latitude": 56.6992,
"launch_longitude": 38.8247,
"launch_datetime": "2025-06-25T20-45-000Z",
"launch_altitude": 0,
"profile": "standard_profile",
"ascent_rate": 5,
"burst_altitude": 30000,
"descent_rate": 5,
"format": "json"
}
READY_URL = "http://localhost:8080/ready"
# --- Utility functions ---
def run_compose_up():
print("[INFO] Running docker-compose down --remove-orphans ...")
result = subprocess.run(["docker-compose", "down", "--remove-orphans"], capture_output=True)
if result.returncode != 0:
print("[ERROR] docker-compose down failed:", result.stderr.decode())
sys.exit(1)
print("[INFO] docker-compose down completed.")
print("[INFO] Running docker-compose up -d ...")
result = subprocess.run(["docker-compose", "up", "-d"], capture_output=True)
if result.returncode != 0:
print("[ERROR] docker-compose up failed:", result.stderr.decode())
sys.exit(1)
print("[INFO] docker-compose up -d completed.")
return True
def wait_for_ready(timeout=900):
print(f"[INFO] Waiting for {READY_URL} to be ready ...")
start = time.time()
while time.time() - start < timeout:
try:
resp = requests.get(READY_URL, timeout=10)
if resp.status_code == 200:
data = resp.json()
if data.get("status") == "ok":
print("[INFO] Service is ready.")
return
else:
print(f"[INFO] Not ready yet: {data}")
else:
print(f"[INFO] /ready returned status {resp.status_code}")
except Exception as e:
print(f"[INFO] Exception while polling /ready: {e}")
time.sleep(10)
print(f"[ERROR] Service did not become ready in {timeout} seconds.")
sys.exit(1)
def fetch_reference():
print(f"[INFO] Fetching reference prediction from {REFERENCE_API_URL}")
resp = requests.get(REFERENCE_API_URL, timeout=60)
if resp.status_code != 200:
print(f"[ERROR] Reference API returned {resp.status_code}: {resp.text}")
sys.exit(1)
return resp.json()
def fetch_local():
print(f"[INFO] Fetching local prediction from {LOCAL_API_URL}")
resp = requests.get(LOCAL_API_URL, timeout=60)
if resp.status_code != 200:
print(f"[ERROR] Local API returned {resp.status_code}: {resp.text}")
sys.exit(1)
return resp.json()
def haversine(lat1, lon1, lat2, lon2):
"""Calculate the great-circle distance between two points on the Earth (specified in decimal degrees). Returns distance in kilometers."""
R = 6371.0 # Earth radius in kilometers
lat1, lon1, lat2, lon2 = map(math.radians, [lat1, lon1, lat2, lon2])
dlat = lat2 - lat1
dlon = lon2 - lon1
a = math.sin(dlat/2)**2 + math.cos(lat1) * math.cos(lat2) * math.sin(dlon/2)**2
c = 2 * math.atan2(math.sqrt(a), math.sqrt(1 - a))
return R * c
def compare_results(reference_data, local_data):
"""Compare prediction results between reference and local APIs."""
print("[INFO] Comparing results ...")
# Extract trajectory data
ref_trajectory = reference_data.get('prediction', [{}])[0].get('trajectory', [])
local_trajectory = local_data.get('prediction', [{}])[0].get('trajectory', [])
print(f"[DEBUG] Reference trajectory length: {len(ref_trajectory)}")
print(f"[DEBUG] Local trajectory length: {len(local_trajectory)}")
# Show first 3 points from both APIs
print("\n[DEBUG] First 3 points - Reference API:")
for i, point in enumerate(ref_trajectory[:3]):
print(f" [{i}] alt={point.get('altitude', 'N/A')}, lat={point.get('latitude', 'N/A')}, lon={point.get('longitude', 'N/A')}, time={point.get('datetime', 'N/A')}")
print("\n[DEBUG] First 3 points - Local API:")
for i, point in enumerate(local_trajectory[:3]):
print(f" [{i}] alt={point.get('altitude', 'N/A')}, lat={point.get('latitude', 'N/A')}, lon={point.get('longitude', 'N/A')}, time={point.get('datetime', 'N/A')}")
# Show last 3 points from both APIs
print("\n[DEBUG] Last 3 points - Reference API:")
for i, point in enumerate(ref_trajectory[-3:]):
idx = len(ref_trajectory) - 3 + i
print(f" [{idx}] alt={point.get('altitude', 'N/A')}, lat={point.get('latitude', 'N/A')}, lon={point.get('longitude', 'N/A')}, time={point.get('datetime', 'N/A')}")
print("\n[DEBUG] Last 3 points - Local API:")
for i, point in enumerate(local_trajectory[-3:]):
idx = len(local_trajectory) - 3 + i
print(f" [{idx}] alt={point.get('altitude', 'N/A')}, lat={point.get('latitude', 'N/A')}, lon={point.get('longitude', 'N/A')}, time={point.get('datetime', 'N/A')}")
# Compare trajectory lengths
if len(ref_trajectory) != len(local_trajectory):
print(f"[DIFF] Trajectory length mismatch: {len(local_trajectory)} vs {len(ref_trajectory)}")
return False
# Compare trajectory points and calculate drift
min_len = min(len(ref_trajectory), len(local_trajectory))
max_drift = 0.0
max_drift_idx = -1
drift_list = []
print("\n[DRIFT] Trajectory point-by-point distance (km):")
for i in range(min_len):
ref_point = ref_trajectory[i]
local_point = local_trajectory[i]
ref_lat = ref_point.get('latitude')
ref_lon = ref_point.get('longitude')
local_lat = local_point.get('latitude')
local_lon = local_point.get('longitude')
drift_km = None
if None not in (ref_lat, ref_lon, local_lat, local_lon):
drift_km = haversine(ref_lat, ref_lon, local_lat, local_lon)
drift_list.append(drift_km)
if drift_km > max_drift:
max_drift = drift_km
max_drift_idx = i
print(f" [{i}] Drift: {drift_km:.3f} km")
else:
print(f" [{i}] Drift: N/A (missing lat/lon)")
if drift_list:
mean_drift = sum(drift_list) / len(drift_list)
print(f"\n[DRIFT] Max drift: {max_drift:.3f} km at idx {max_drift_idx}")
print(f"[DRIFT] Mean drift: {mean_drift:.3f} km over {len(drift_list)} points")
else:
print("[DRIFT] No valid drift data to report.")
# Continue with original comparison for altitude, etc.
for i in range(min_len):
ref_point = ref_trajectory[i]
local_point = local_trajectory[i]
for key in ['altitude', 'latitude', 'longitude']:
ref_val = ref_point.get(key)
local_val = local_point.get(key)
if ref_val is not None and local_val is not None:
if abs(ref_val - local_val) > 0.1:
print(f"[DIFF] At idx {i}, key {key}: {local_val} != {ref_val}")
return False
print("[SUCCESS] Results match!")
return True
def test_custom_profile():
"""Test custom profile with base64 encoded curve."""
print("\n[TEST] Testing custom_profile...")
# Create a simple custom ascent curve (altitude vs time in seconds)
curve_data = {
"altitude": [0, 30000],
"time": [0, 6000]
}
curve_b64 = base64.b64encode(json.dumps(curve_data).encode()).decode()
# Test parameters for custom profile
params = {
"launch_latitude": 56.6992,
"launch_longitude": 38.8247,
"launch_datetime": "2025-06-25T13:28:00Z",
"launch_altitude": 0,
"profile": "custom_profile",
"ascent_curve": curve_b64
}
try:
# Test local API (use GET)
local_resp = requests.get(
"http://localhost:8080/api/v1/prediction",
params=params,
timeout=30
)
local_resp.raise_for_status()
local_data = local_resp.json()
print(f"[INFO] Custom profile test - Local API returned {len(local_data.get('prediction', [{}])[0].get('trajectory', []))} trajectory points")
return True
except Exception as e:
print(f"[ERROR] Custom profile test failed: {e}")
return False
def test_all_profiles():
"""Test all available profiles."""
profiles = [
("standard_profile", "Standard profile test"),
("float_profile", "Float profile test"),
("reverse_profile", "Reverse profile test"),
("custom_profile", "Custom profile test")
]
results = {}
for profile, description in profiles:
print(f"\n[TEST] {description}...")
if profile == "custom_profile":
success = test_custom_profile()
else:
success = test_single_profile(profile)
results[profile] = success
print(f"[RESULT] {profile}: {'PASS' if success else 'FAIL'}")
# Print summary
print("\n" + "="*50)
print("TEST SUMMARY")
print("="*50)
for profile, success in results.items():
status = "PASS" if success else "FAIL"
print(f"{profile:20} : {status}")
total_tests = len(results)
passed_tests = sum(results.values())
print(f"\nTotal tests: {total_tests}, Passed: {passed_tests}, Failed: {total_tests - passed_tests}")
return all(results.values())
def test_single_profile(profile):
"""Test a single profile against reference API."""
# Test parameters
params = {
"launch_latitude": 56.6992,
"launch_longitude": 38.8247,
"launch_datetime": "2025-06-25T13:28:00Z",
"launch_altitude": 0,
"profile": profile,
"ascent_rate": 5,
"burst_altitude": 30000,
"descent_rate": 5
}
# Add float altitude for float profile
if profile == "float_profile":
params["float_altitude"] = 25000
try:
# Test local API (use GET)
local_resp = requests.get(
"http://localhost:8080/api/v1/prediction",
params=params,
timeout=30
)
local_resp.raise_for_status()
local_data = local_resp.json()
print(f"[INFO] {profile} - Local API returned {len(local_data.get('prediction', [{}])[0].get('trajectory', []))} trajectory points")
return True
except Exception as e:
print(f"[ERROR] {profile} test failed: {e}")
return False
def main():
"""Main test function."""
print("[INFO] Starting comprehensive predictor API tests...")
# Run the original standard profile test
print("\n[TEST] Running original standard_profile test...")
run_compose_up()
wait_for_ready()
ref = fetch_reference()
local = fetch_local()
print("[INFO] Comparing results ...")
original_success = compare_results(ref, local)
if original_success:
print("[SUCCESS] Original standard_profile test passed!")
else:
print("[FAIL] Original standard_profile test failed!")
# Test all profiles
print("\n[TEST] Running all profile tests...")
all_profiles_success = test_all_profiles()
# Final result
overall_success = original_success and all_profiles_success
print(f"\n[FINAL RESULT] Overall: {'PASS' if overall_success else 'FAIL'}")
if overall_success:
sys.exit(0)
else:
sys.exit(1)
if __name__ == "__main__":
main()

View file

@ -1,89 +0,0 @@
package main
import (
"context"
"fmt"
"log"
"os"
"time"
"git.intra.yksa.space/gsn/predictor/internal/pkg/grib"
)
func main() {
ctx := context.Background()
// Create S3 downloader
downloader, err := grib.NewS3Downloader(
"/tmp/grib_test",
4, // parallel downloads
"noaa-gfs-bdp-pds",
"us-east-1",
)
if err != nil {
log.Fatalf("Failed to create S3 downloader: %v", err)
}
// Ensure directory exists
if err := os.MkdirAll("/tmp/grib_test", 0o755); err != nil {
log.Fatalf("Failed to create directory: %v", err)
}
// Find nearest run (6-hour intervals: 00, 06, 12, 18 UTC)
now := time.Now().UTC()
hour := now.Hour() - (now.Hour() % 6)
// Use data from 6 hours ago to ensure it's available
run := time.Date(now.Year(), now.Month(), now.Day(), hour, 0, 0, 0, time.UTC).Add(-6 * time.Hour)
fmt.Printf("Testing S3 download for run: %s\n", run.Format("2006-01-02 15:04 MST"))
// List available files first
runStr := run.Format("20060102")
fmt.Printf("Listing available files for %s/%02d...\n", runStr, run.Hour())
files, err := downloader.ListAvailableFiles(ctx, runStr, run.Hour())
if err != nil {
log.Fatalf("Failed to list files: %v", err)
}
fmt.Printf("Found %d files in S3:\n", len(files))
if len(files) > 0 {
// Show first 5 files
for i, file := range files {
if i >= 5 {
fmt.Printf("... and %d more files\n", len(files)-5)
break
}
fmt.Printf(" - %s\n", file)
}
}
// Try downloading just first 3 forecast hours (f000, f001, f002)
fmt.Println("\nTesting download of first 3 forecast hours...")
testRun := run
// Create a timeout context for the download
downloadCtx, cancel := context.WithTimeout(ctx, 5*time.Minute)
defer cancel()
if err := downloader.Run(downloadCtx, testRun); err != nil {
log.Fatalf("Failed to download: %v", err)
}
fmt.Println("\nDownload completed successfully!")
// Check downloaded files
entries, err := os.ReadDir("/tmp/grib_test")
if err != nil {
log.Fatalf("Failed to read directory: %v", err)
}
fmt.Printf("\nDownloaded %d files:\n", len(entries))
for i, entry := range entries {
if i >= 10 {
fmt.Printf("... and %d more files\n", len(entries)-10)
break
}
info, _ := entry.Info()
fmt.Printf(" - %s (%.2f MB)\n", entry.Name(), float64(info.Size())/1024/1024)
}
}

View file

@ -1,68 +0,0 @@
package main
import (
"context"
"fmt"
"io"
"log"
"os"
"github.com/aws/aws-sdk-go-v2/aws"
"github.com/aws/aws-sdk-go-v2/config"
"github.com/aws/aws-sdk-go-v2/service/s3"
)
func main() {
ctx := context.Background()
// Create AWS config with anonymous credentials
cfg, err := config.LoadDefaultConfig(ctx,
config.WithRegion("us-east-1"),
config.WithCredentialsProvider(aws.AnonymousCredentials{}),
)
if err != nil {
log.Fatalf("Failed to load config: %v", err)
}
client := s3.NewFromConfig(cfg)
// Try to download a single file
bucket := "noaa-gfs-bdp-pds"
key := "gfs.20251020/00/atmos/gfs.t00z.pgrb2.0p50.f000"
fmt.Printf("Downloading: s3://%s/%s\n", bucket, key)
input := &s3.GetObjectInput{
Bucket: aws.String(bucket),
Key: aws.String(key),
}
result, err := client.GetObject(ctx, input)
if err != nil {
log.Fatalf("Failed to get object: %v", err)
}
defer result.Body.Close()
// Create output file
outFile := "/tmp/test_grib.part"
f, err := os.Create(outFile)
if err != nil {
log.Fatalf("Failed to create file: %v", err)
}
defer f.Close()
// Copy data
written, err := io.Copy(f, result.Body)
if err != nil {
log.Fatalf("Failed to copy data: %v (wrote %d bytes)", err, written)
}
fmt.Printf("Successfully downloaded %d bytes\n", written)
// Rename
if err := os.Rename(outFile, "/tmp/test_grib"); err != nil {
log.Fatalf("Failed to rename: %v", err)
}
fmt.Println("Download complete!")
}