688 lines
26 KiB
Python
688 lines
26 KiB
Python
#!/usr/bin/env python3
|
|
from __future__ import annotations
|
|
|
|
import argparse
|
|
import csv
|
|
import gzip
|
|
import json
|
|
import math
|
|
import re
|
|
import struct
|
|
from datetime import datetime
|
|
from pathlib import Path
|
|
from typing import Any
|
|
|
|
import numpy as np
|
|
|
|
|
|
def read_gamma_value(path: Path, key: str) -> str:
|
|
for line in path.read_text(encoding="utf-8", errors="replace").splitlines():
|
|
parts = line.split()
|
|
if parts and parts[0].rstrip(":") == key.rstrip(":"):
|
|
return parts[1]
|
|
raise KeyError(f"{key} not found in {path}")
|
|
|
|
|
|
def read_float32(path: Path, shape: tuple[int, int] | None = None) -> np.ndarray:
|
|
data = np.fromfile(path, dtype=">f4")
|
|
if shape is not None:
|
|
data = data.reshape(shape)
|
|
return data
|
|
|
|
|
|
def read_float32_pixel(path: Path, width: int, x: int, y: int) -> float:
|
|
with path.open("rb") as handle:
|
|
handle.seek((y * width + x) * 4)
|
|
chunk = handle.read(4)
|
|
if len(chunk) != 4:
|
|
return float("nan")
|
|
return float(struct.unpack(">f", chunk)[0])
|
|
|
|
|
|
def write_scaled_float32(input_path: Path, output_path: Path, scale: float) -> None:
|
|
data = np.fromfile(input_path, dtype=">f4")
|
|
output_path.parent.mkdir(parents=True, exist_ok=True)
|
|
(data * float(scale)).astype(">f4", copy=False).tofile(output_path)
|
|
|
|
|
|
def monitor_valid_mask(rate: np.ndarray, sigma: np.ndarray) -> np.ndarray:
|
|
lines, width = rate.shape
|
|
yy, xx = np.indices(rate.shape)
|
|
edge_mask = (
|
|
(xx > width * 0.1)
|
|
& (xx < width * 0.9)
|
|
& (yy > lines * 0.1)
|
|
& (yy < lines * 0.9)
|
|
)
|
|
finite = np.isfinite(rate) & np.isfinite(sigma)
|
|
valid = finite & edge_mask & (rate != 0.0) & (sigma > 0.0)
|
|
if not valid.any():
|
|
raise RuntimeError("No valid pixels available for monitor point selection")
|
|
return valid
|
|
|
|
|
|
def remove_near_selected(candidate: np.ndarray, selected: list[tuple[int, int]], min_distance: int) -> np.ndarray:
|
|
if not selected:
|
|
return candidate
|
|
yy, xx = np.indices(candidate.shape)
|
|
filtered = candidate.copy()
|
|
min_distance_sq = float(min_distance * min_distance)
|
|
for x, y in selected:
|
|
filtered &= ((xx - float(x)) ** 2 + (yy - float(y)) ** 2) >= min_distance_sq
|
|
return filtered
|
|
|
|
|
|
def pick_scored_point(
|
|
score: np.ndarray,
|
|
candidate: np.ndarray,
|
|
selected: list[tuple[int, int]],
|
|
*,
|
|
min_distance: int,
|
|
) -> tuple[int, int] | None:
|
|
filtered = remove_near_selected(candidate, selected, min_distance)
|
|
if not filtered.any():
|
|
filtered = candidate
|
|
if not filtered.any():
|
|
return None
|
|
safe_score = np.full(score.shape, -np.inf, dtype=np.float64)
|
|
safe_score[filtered] = score[filtered]
|
|
y, x = np.unravel_index(int(np.nanargmax(safe_score)), score.shape)
|
|
if not np.isfinite(safe_score[y, x]):
|
|
return None
|
|
return int(x), int(y)
|
|
|
|
|
|
def pick_auto_point(rate: np.ndarray, sigma: np.ndarray) -> tuple[int, int]:
|
|
valid = monitor_valid_mask(rate, sigma)
|
|
|
|
abs_rate = np.abs(rate[valid])
|
|
sig = sigma[valid]
|
|
rate_min = np.percentile(abs_rate, 85)
|
|
rate_max = np.percentile(abs_rate, 99)
|
|
sigma_max = np.percentile(sig, 40)
|
|
candidate = valid & (np.abs(rate) >= rate_min) & (np.abs(rate) <= rate_max) & (sigma <= sigma_max)
|
|
if not candidate.any():
|
|
candidate = valid
|
|
|
|
score = np.zeros(rate.shape, dtype=np.float32)
|
|
score[candidate] = np.abs(rate[candidate]) / (sigma[candidate] + 1.0e-6)
|
|
y, x = np.unravel_index(int(np.argmax(score)), rate.shape)
|
|
return int(x), int(y)
|
|
|
|
|
|
def pick_auto_points(rate: np.ndarray, sigma: np.ndarray, *, count: int = 5) -> list[dict[str, Any]]:
|
|
valid = monitor_valid_mask(rate, sigma)
|
|
lines, width = rate.shape
|
|
yy, xx = np.indices(rate.shape)
|
|
min_distance = max(24, int(min(width, lines) * 0.08))
|
|
sigma_max = float(np.percentile(sigma[valid], 40))
|
|
low_sigma = valid & (sigma <= sigma_max)
|
|
if not low_sigma.any():
|
|
low_sigma = valid
|
|
|
|
abs_rate = np.abs(rate)
|
|
abs_valid = abs_rate[valid]
|
|
high_abs_min = float(np.percentile(abs_valid, 85))
|
|
high_abs_max = float(np.percentile(abs_valid, 99))
|
|
low_abs_max = float(np.percentile(abs_valid, 25))
|
|
cx = (width - 1) / 2.0
|
|
cy = (lines - 1) / 2.0
|
|
|
|
definitions = [
|
|
{
|
|
"point_id": "auto_away_high_rate_low_sigma",
|
|
"selection": "automatic_away_from_radar_high_rate_low_sigma_non_edge",
|
|
"candidate": low_sigma & (rate < 0.0) & (abs_rate >= high_abs_min) & (abs_rate <= high_abs_max),
|
|
"score": (-rate) / (sigma + 1.0e-6),
|
|
},
|
|
{
|
|
"point_id": "auto_toward_high_rate_low_sigma",
|
|
"selection": "automatic_toward_radar_high_rate_low_sigma_non_edge",
|
|
"candidate": low_sigma & (rate > 0.0) & (abs_rate >= high_abs_min) & (abs_rate <= high_abs_max),
|
|
"score": rate / (sigma + 1.0e-6),
|
|
},
|
|
{
|
|
"point_id": "auto_low_sigma_high_rate",
|
|
"selection": "automatic_low_sigma_high_abs_rate_non_edge",
|
|
"candidate": low_sigma & (abs_rate >= high_abs_min) & (abs_rate <= high_abs_max),
|
|
"score": abs_rate / (sigma + 1.0e-6),
|
|
},
|
|
{
|
|
"point_id": "auto_stable_low_sigma",
|
|
"selection": "automatic_near_zero_rate_low_sigma_non_edge",
|
|
"candidate": low_sigma & (abs_rate <= low_abs_max),
|
|
"score": 1.0 / ((abs_rate + 1.0) * (sigma + 1.0e-6)),
|
|
},
|
|
{
|
|
"point_id": "auto_center_valid",
|
|
"selection": "automatic_valid_pixel_nearest_stack_center",
|
|
"candidate": valid,
|
|
"score": -((xx - cx) ** 2 + (yy - cy) ** 2),
|
|
},
|
|
]
|
|
|
|
selected_xy: list[tuple[int, int]] = []
|
|
selected_points: list[dict[str, Any]] = []
|
|
for definition in definitions:
|
|
if len(selected_points) >= count:
|
|
break
|
|
candidate = definition["candidate"]
|
|
if not candidate.any():
|
|
candidate = low_sigma if low_sigma.any() else valid
|
|
picked = pick_scored_point(
|
|
np.asarray(definition["score"], dtype=np.float64),
|
|
candidate,
|
|
selected_xy,
|
|
min_distance=min_distance,
|
|
)
|
|
if picked is None:
|
|
continue
|
|
x, y = picked
|
|
selected_xy.append((x, y))
|
|
selected_points.append(
|
|
{
|
|
"point_id": definition["point_id"],
|
|
"selection": definition["selection"],
|
|
"range_pixel": x,
|
|
"azimuth_line": y,
|
|
}
|
|
)
|
|
|
|
if not selected_points:
|
|
x, y = pick_auto_point(rate, sigma)
|
|
selected_points.append(
|
|
{
|
|
"point_id": "auto_low_sigma_high_rate",
|
|
"selection": "automatic_low_sigma_high_rate_non_edge",
|
|
"range_pixel": x,
|
|
"azimuth_line": y,
|
|
}
|
|
)
|
|
return selected_points[:count]
|
|
|
|
|
|
def dem_grid(dem_par: Path) -> dict[str, float | int]:
|
|
return {
|
|
"width": int(read_gamma_value(dem_par, "width")),
|
|
"nlines": int(read_gamma_value(dem_par, "nlines")),
|
|
"corner_lon": float(read_gamma_value(dem_par, "corner_lon")),
|
|
"corner_lat": float(read_gamma_value(dem_par, "corner_lat")),
|
|
"post_lon": float(read_gamma_value(dem_par, "post_lon")),
|
|
"post_lat": float(read_gamma_value(dem_par, "post_lat")),
|
|
}
|
|
|
|
|
|
def radar_to_lonlat(x: int, y: int, dem_par: Path, lookup: Path) -> tuple[float | None, float | None]:
|
|
grid = dem_grid(dem_par)
|
|
width = int(grid["width"])
|
|
lines = int(grid["nlines"])
|
|
lut = np.fromfile(lookup, dtype=">c8").reshape((lines, width))
|
|
rng = lut.real
|
|
az = lut.imag
|
|
valid = np.isfinite(rng) & np.isfinite(az) & (rng > 0.0) & (az > 0.0)
|
|
if not valid.any():
|
|
return None, None
|
|
|
|
distance = np.full(rng.shape, np.inf, dtype=np.float32)
|
|
distance[valid] = (rng[valid] - float(x)) ** 2 + (az[valid] - float(y)) ** 2
|
|
gy, gx = np.unravel_index(int(np.argmin(distance)), distance.shape)
|
|
lon = float(grid["corner_lon"]) + (gx + 0.5) * float(grid["post_lon"])
|
|
lat = float(grid["corner_lat"]) + (gy + 0.5) * float(grid["post_lat"])
|
|
return float(lon), float(lat)
|
|
|
|
|
|
def lonlat_to_radar(lon: float, lat: float, dem_par: Path, lookup: Path) -> tuple[int, int]:
|
|
grid = dem_grid(dem_par)
|
|
width = int(grid["width"])
|
|
lines = int(grid["nlines"])
|
|
gx = int(round((lon - float(grid["corner_lon"])) / float(grid["post_lon"]) - 0.5))
|
|
gy = int(round((lat - float(grid["corner_lat"])) / float(grid["post_lat"]) - 0.5))
|
|
gx = max(0, min(width - 1, gx))
|
|
gy = max(0, min(lines - 1, gy))
|
|
lut = np.fromfile(lookup, dtype=">c8").reshape((lines, width))
|
|
value = lut[gy, gx]
|
|
if not (np.isfinite(value.real) and np.isfinite(value.imag) and value.real > 0 and value.imag > 0):
|
|
raise RuntimeError(f"manual lon/lat maps to invalid lookup pixel: lon={lon}, lat={lat}")
|
|
return int(round(float(value.real))), int(round(float(value.imag)))
|
|
|
|
|
|
def safe_point_id(value: str, fallback: str) -> str:
|
|
text = str(value or "").strip() or fallback
|
|
text = re.sub(r"[^A-Za-z0-9_.-]+", "_", text)[:64]
|
|
return text or fallback
|
|
|
|
|
|
def load_diff_files(timeseries_dir: Path) -> list[Path]:
|
|
tab = timeseries_dir / "diff_ts.tab"
|
|
if tab.is_file():
|
|
rows = [line.strip() for line in tab.read_text(encoding="utf-8", errors="replace").splitlines() if line.strip()]
|
|
files = [Path(row.split()[0]) for row in rows if row.split()]
|
|
files = [path for path in files if path.is_file()]
|
|
if files:
|
|
return files
|
|
return sorted(timeseries_dir.glob("diff_ts_*.diff"))
|
|
|
|
|
|
def point_records(
|
|
diff_files: list[Path],
|
|
*,
|
|
dates: list[str],
|
|
width: int,
|
|
x: int,
|
|
y: int,
|
|
scale_mm: float,
|
|
) -> list[dict[str, Any]]:
|
|
records: list[dict[str, Any]] = []
|
|
for index, path in enumerate(diff_files):
|
|
phase = read_float32_pixel(path, width, x, y)
|
|
away_mm = float(phase * scale_mm) if math.isfinite(phase) else float("nan")
|
|
records.append(
|
|
{
|
|
"date": dates[index] if index < len(dates) else f"step_{index + 1:03d}",
|
|
"phase_rad": float(phase),
|
|
"los_away_mm": away_mm,
|
|
"los_toward_mm": -away_mm if math.isfinite(away_mm) else float("nan"),
|
|
}
|
|
)
|
|
return records
|
|
|
|
|
|
def write_point_outputs(
|
|
point_dir: Path,
|
|
point: dict[str, Any],
|
|
*,
|
|
records: list[dict[str, Any]],
|
|
rate_value: float,
|
|
sigma_value: float,
|
|
wavelength: float,
|
|
reference_date: str,
|
|
) -> dict[str, str]:
|
|
import matplotlib
|
|
|
|
matplotlib.use("Agg")
|
|
import matplotlib.pyplot as plt
|
|
|
|
point_id = str(point["point_id"])
|
|
csv_path = point_dir / f"{point_id}_timeseries.csv"
|
|
json_path = point_dir / f"{point_id}_metadata.json"
|
|
png_path = point_dir / f"{point_id}_timeseries.png"
|
|
|
|
with csv_path.open("w", newline="", encoding="utf-8") as handle:
|
|
writer = csv.DictWriter(handle, fieldnames=["date", "phase_rad", "los_away_mm", "los_toward_mm"])
|
|
writer.writeheader()
|
|
writer.writerows(records)
|
|
|
|
metadata = {
|
|
"schema": "insar.sbas-monitor-point/v1",
|
|
"point_id": point_id,
|
|
"selection": point.get("selection"),
|
|
"radar_pixel": {"range": int(point["range_pixel"]), "azimuth": int(point["azimuth_line"])},
|
|
"approx_lonlat": {"lon": point.get("lon"), "lat": point.get("lat")},
|
|
"reference_date": reference_date,
|
|
"los_convention": "toward radar positive; away from radar negative",
|
|
"los_rate_toward_mm_per_year": rate_value,
|
|
"los_sigma_mm_per_year": sigma_value,
|
|
"wavelength_m": wavelength,
|
|
"records": records,
|
|
}
|
|
json_path.write_text(json.dumps(metadata, indent=2, ensure_ascii=False), encoding="utf-8")
|
|
|
|
dates = [record["date"] for record in records]
|
|
disp = [record["los_toward_mm"] for record in records]
|
|
plt.figure(figsize=(8.0, 4.6), dpi=160)
|
|
plt.plot(dates, disp, marker="o", linewidth=2.0, color="#1f77b4")
|
|
plt.axhline(0, color="#666666", linewidth=0.8)
|
|
plt.grid(True, color="#dddddd", linewidth=0.7)
|
|
plt.title(
|
|
f"LOS displacement time series ({point_id})\n"
|
|
f"toward radar positive, rate={rate_value:.2f} mm/yr, sigma={sigma_value:.2f} mm/yr",
|
|
fontsize=10,
|
|
)
|
|
plt.xlabel("Date")
|
|
plt.ylabel("LOS displacement (mm)")
|
|
plt.tight_layout()
|
|
plt.savefig(png_path)
|
|
plt.close()
|
|
return {"png": str(png_path), "csv": str(csv_path), "metadata": str(json_path)}
|
|
|
|
|
|
def read_geotiff_float32(path: Path) -> dict[str, Any]:
|
|
try:
|
|
import rasterio
|
|
|
|
with rasterio.open(path) as src:
|
|
transform = src.transform
|
|
return {
|
|
"array": src.read(1).astype(np.float32, copy=False),
|
|
"width": src.width,
|
|
"height": src.height,
|
|
"nodata": src.nodata,
|
|
"crs": src.crs.to_string() if src.crs else None,
|
|
"transform": (transform.a, transform.b, transform.c, transform.d, transform.e, transform.f),
|
|
}
|
|
except Exception as rasterio_exc:
|
|
try:
|
|
from osgeo import gdal
|
|
except Exception as gdal_exc:
|
|
raise RuntimeError("rasterio or osgeo.gdal is required to read GeoTIFF files") from gdal_exc
|
|
dataset = gdal.Open(str(path), gdal.GA_ReadOnly)
|
|
if dataset is None:
|
|
raise RuntimeError(f"Unable to open GeoTIFF: {path}") from rasterio_exc
|
|
band = dataset.GetRasterBand(1)
|
|
array = band.ReadAsArray().astype(np.float32, copy=False)
|
|
geotransform = dataset.GetGeoTransform()
|
|
return {
|
|
"array": array,
|
|
"width": int(dataset.RasterXSize),
|
|
"height": int(dataset.RasterYSize),
|
|
"nodata": band.GetNoDataValue(),
|
|
"crs": dataset.GetProjection() or None,
|
|
"transform": (
|
|
float(geotransform[1]),
|
|
float(geotransform[2]),
|
|
float(geotransform[0]),
|
|
float(geotransform[4]),
|
|
float(geotransform[5]),
|
|
float(geotransform[3]),
|
|
),
|
|
}
|
|
|
|
|
|
def normalize_crs_label(value: Any) -> str | None:
|
|
text = str(value or "").strip()
|
|
if not text:
|
|
return None
|
|
upper = text.upper()
|
|
if "EPSG" in upper and "4326" in upper:
|
|
return "EPSG:4326"
|
|
if "WGS 84" in upper or "WGS_1984" in upper:
|
|
return "EPSG:4326"
|
|
return text[:240]
|
|
|
|
|
|
def pixel_center(transform: tuple[float, float, float, float, float, float], row: int, col: int) -> tuple[float, float]:
|
|
a, b, c, d, e, f = transform
|
|
x = c + (col + 0.5) * a + (row + 0.5) * b
|
|
y = f + (col + 0.5) * d + (row + 0.5) * e
|
|
return float(x), float(y)
|
|
|
|
|
|
def run_export_points_geojson(args: argparse.Namespace) -> int:
|
|
|
|
toward_path = Path(args.toward_tif)
|
|
away_path = Path(args.away_tif)
|
|
sigma_path = Path(args.sigma_tif)
|
|
output_path = Path(args.output)
|
|
summary_path = Path(args.summary_path)
|
|
output_path.parent.mkdir(parents=True, exist_ok=True)
|
|
summary_path.parent.mkdir(parents=True, exist_ok=True)
|
|
|
|
run_id = str(args.run_id or "").strip()
|
|
date_start = str(args.date_start or "").strip()
|
|
date_end = str(args.date_end or "").strip()
|
|
reference_date = str(args.reference_date or "").strip()
|
|
admin_province = str(args.admin_province or "").strip()
|
|
admin_city = str(args.admin_city or "").strip()
|
|
|
|
fields = [
|
|
"run_id",
|
|
"row",
|
|
"col",
|
|
"lon",
|
|
"lat",
|
|
"los_rate_toward_mm_per_year",
|
|
"los_rate_away_mm_per_year",
|
|
"los_sigma_mm_per_year",
|
|
"date_start",
|
|
"date_end",
|
|
"reference_date",
|
|
"admin_province",
|
|
"admin_city",
|
|
]
|
|
|
|
toward_meta = read_geotiff_float32(toward_path)
|
|
away_meta = read_geotiff_float32(away_path)
|
|
sigma_meta = read_geotiff_float32(sigma_path)
|
|
width = int(toward_meta["width"])
|
|
height = int(toward_meta["height"])
|
|
if (width, height) != (int(away_meta["width"]), int(away_meta["height"])):
|
|
raise RuntimeError("toward and away GeoTIFF dimensions do not match")
|
|
if (width, height) != (int(sigma_meta["width"]), int(sigma_meta["height"])):
|
|
raise RuntimeError("toward and sigma GeoTIFF dimensions do not match")
|
|
|
|
toward = np.asarray(toward_meta["array"], dtype=np.float32)
|
|
away = np.asarray(away_meta["array"], dtype=np.float32)
|
|
sigma = np.asarray(sigma_meta["array"], dtype=np.float32)
|
|
valid = np.isfinite(toward) & np.isfinite(away) & np.isfinite(sigma) & (sigma > 0.0)
|
|
if toward_meta.get("nodata") is not None:
|
|
valid &= toward != float(toward_meta["nodata"])
|
|
if away_meta.get("nodata") is not None:
|
|
valid &= away != float(away_meta["nodata"])
|
|
if sigma_meta.get("nodata") is not None:
|
|
valid &= sigma != float(sigma_meta["nodata"])
|
|
|
|
transform = tuple(float(value) for value in toward_meta["transform"])
|
|
feature_count = 0
|
|
with gzip.open(output_path, "wt", encoding="utf-8", compresslevel=6) as handle:
|
|
handle.write('{"type":"FeatureCollection","features":[\n')
|
|
first = True
|
|
for row in range(height):
|
|
cols = np.where(valid[row])[0]
|
|
for col in cols.tolist():
|
|
lon, lat = pixel_center(transform, row, int(col))
|
|
properties = {
|
|
"run_id": run_id,
|
|
"row": int(row),
|
|
"col": int(col),
|
|
"lon": lon,
|
|
"lat": lat,
|
|
"los_rate_toward_mm_per_year": float(toward[row, col]),
|
|
"los_rate_away_mm_per_year": float(away[row, col]),
|
|
"los_sigma_mm_per_year": float(sigma[row, col]),
|
|
"date_start": date_start,
|
|
"date_end": date_end,
|
|
"reference_date": reference_date,
|
|
"admin_province": admin_province,
|
|
"admin_city": admin_city,
|
|
}
|
|
feature = {
|
|
"type": "Feature",
|
|
"geometry": {"type": "Point", "coordinates": [lon, lat]},
|
|
"properties": properties,
|
|
}
|
|
if not first:
|
|
handle.write(",\n")
|
|
handle.write(json.dumps(feature, ensure_ascii=False, separators=(",", ":")))
|
|
first = False
|
|
feature_count += 1
|
|
handle.write("\n]}\n")
|
|
|
|
crs = normalize_crs_label(toward_meta.get("crs"))
|
|
|
|
summary = {
|
|
"schema": "insar.gamma-sbas-point-vector-summary/v1",
|
|
"generated_at": datetime.utcnow().isoformat(timespec="seconds") + "Z",
|
|
"ready": output_path.is_file() and output_path.stat().st_size > 0,
|
|
"feature_count": feature_count,
|
|
"output_geojson_gz": str(output_path),
|
|
"output_size_bytes": output_path.stat().st_size if output_path.is_file() else 0,
|
|
"fields": fields,
|
|
"source_geotiffs": {
|
|
"los_rate_toward_mm_per_year": str(toward_path),
|
|
"los_rate_away_mm_per_year": str(away_path),
|
|
"los_sigma_mm_per_year": str(sigma_path),
|
|
},
|
|
"width": width,
|
|
"height": height,
|
|
"crs": crs,
|
|
"date_start": date_start,
|
|
"date_end": date_end,
|
|
"reference_date": reference_date,
|
|
"admin_region": {
|
|
"province": admin_province or None,
|
|
"city": admin_city or None,
|
|
},
|
|
"los_convention": "toward radar positive; away from radar negative",
|
|
"frontend_policy": "download_only; do not render full point GeoJSON in browser",
|
|
}
|
|
summary_path.write_text(json.dumps(summary, indent=2, ensure_ascii=False), encoding="utf-8")
|
|
print(json.dumps(summary, indent=2, ensure_ascii=False))
|
|
return 0
|
|
|
|
|
|
def run_phase_to_los(args: argparse.Namespace) -> int:
|
|
write_scaled_float32(Path(args.input), Path(args.output), float(args.scale))
|
|
return 0
|
|
|
|
|
|
def run_monitor_points(args: argparse.Namespace) -> int:
|
|
timeseries_dir = Path(args.timeseries_dir)
|
|
export_dir = Path(args.export_dir)
|
|
point_dir = Path(args.point_dir)
|
|
mli_par = Path(args.mli_par)
|
|
slc_par = Path(args.slc_par) if args.slc_par else mli_par
|
|
dem_par = Path(args.dem_par)
|
|
lookup = Path(args.lookup)
|
|
monitor_config_path = Path(args.monitor_config)
|
|
summary_path = Path(args.summary_path)
|
|
point_dir.mkdir(parents=True, exist_ok=True)
|
|
|
|
width = int(read_gamma_value(mli_par, "range_samples"))
|
|
lines = int(read_gamma_value(mli_par, "azimuth_lines"))
|
|
shape = (lines, width)
|
|
dates = [item.strip() for item in str(args.dates or "").split(",") if item.strip()]
|
|
reference_date = str(args.reference_date or "").strip()
|
|
|
|
radar_freq = float(read_gamma_value(slc_par, "radar_frequency"))
|
|
wavelength = 299792458.0 / radar_freq
|
|
scale_mm = wavelength / (4.0 * math.pi) * 1000.0
|
|
|
|
rate_toward = read_float32(export_dir / "los_rate_toward_mm_per_year.rdc", shape)
|
|
sigma = read_float32(export_dir / "los_sigma_mm_per_year.rdc", shape)
|
|
diff_files = load_diff_files(timeseries_dir)
|
|
if not diff_files:
|
|
raise RuntimeError(f"No diff_ts files found in {timeseries_dir}")
|
|
|
|
config = {}
|
|
if monitor_config_path.is_file():
|
|
config = json.loads(monitor_config_path.read_text(encoding="utf-8"))
|
|
mode = str(config.get("mode") or "auto_low_sigma_high_rate")
|
|
selected_points: list[dict[str, Any]] = []
|
|
|
|
if mode == "manual_lonlat" and config.get("points"):
|
|
for index, raw in enumerate(config.get("points") or []):
|
|
lon = float(raw["lon"])
|
|
lat = float(raw["lat"])
|
|
x, y = lonlat_to_radar(lon, lat, dem_par, lookup)
|
|
selected_points.append(
|
|
{
|
|
"point_id": safe_point_id(raw.get("point_id"), f"manual_{index + 1:03d}"),
|
|
"selection": "manual_lonlat_nearest_lookup_pixel",
|
|
"range_pixel": x,
|
|
"azimuth_line": y,
|
|
"lon": lon,
|
|
"lat": lat,
|
|
}
|
|
)
|
|
else:
|
|
auto_count = int(config.get("auto_count") or 5)
|
|
for point in pick_auto_points(rate_toward, sigma, count=max(1, min(auto_count, 12))):
|
|
lon, lat = radar_to_lonlat(int(point["range_pixel"]), int(point["azimuth_line"]), dem_par, lookup)
|
|
selected_points.append({**point, "lon": lon, "lat": lat})
|
|
|
|
outputs: list[dict[str, Any]] = []
|
|
for point in selected_points:
|
|
x = int(point["range_pixel"])
|
|
y = int(point["azimuth_line"])
|
|
if not (0 <= x < width and 0 <= y < lines):
|
|
raise ValueError(f"pixel out of bounds: x={x}, y={y}, width={width}, lines={lines}")
|
|
records = point_records(diff_files, dates=dates, width=width, x=x, y=y, scale_mm=scale_mm)
|
|
rate_value = float(rate_toward[y, x])
|
|
sigma_value = float(sigma[y, x])
|
|
files = write_point_outputs(
|
|
point_dir,
|
|
point,
|
|
records=records,
|
|
rate_value=rate_value,
|
|
sigma_value=sigma_value,
|
|
wavelength=wavelength,
|
|
reference_date=reference_date,
|
|
)
|
|
outputs.append(
|
|
{
|
|
**point,
|
|
"los_rate_toward_mm_per_year": rate_value,
|
|
"los_sigma_mm_per_year": sigma_value,
|
|
"record_count": len(records),
|
|
"files": files,
|
|
}
|
|
)
|
|
|
|
summary = {
|
|
"schema": "insar.gamma-sbas-monitor-points-summary/v1",
|
|
"generated_at": datetime.utcnow().isoformat(timespec="seconds") + "Z",
|
|
"ready": bool(outputs),
|
|
"mode": mode,
|
|
"reference_date": reference_date,
|
|
"width": width,
|
|
"lines": lines,
|
|
"wavelength_m": wavelength,
|
|
"diff_ts_count": len(diff_files),
|
|
"date_count": len(dates),
|
|
"monitor_points": outputs,
|
|
}
|
|
summary_path.parent.mkdir(parents=True, exist_ok=True)
|
|
summary_path.write_text(json.dumps(summary, indent=2, ensure_ascii=False), encoding="utf-8")
|
|
print(json.dumps(summary, indent=2, ensure_ascii=False))
|
|
return 0
|
|
|
|
|
|
def build_parser() -> argparse.ArgumentParser:
|
|
parser = argparse.ArgumentParser(description="Gamma SBAS product helper tools")
|
|
subparsers = parser.add_subparsers(dest="command", required=True)
|
|
|
|
phase = subparsers.add_parser("phase-to-los")
|
|
phase.add_argument("input")
|
|
phase.add_argument("output")
|
|
phase.add_argument("scale", type=float)
|
|
phase.set_defaults(func=run_phase_to_los)
|
|
|
|
vector = subparsers.add_parser("export-points-geojson")
|
|
vector.add_argument("--toward-tif", required=True)
|
|
vector.add_argument("--away-tif", required=True)
|
|
vector.add_argument("--sigma-tif", required=True)
|
|
vector.add_argument("--output", required=True)
|
|
vector.add_argument("--summary-path", required=True)
|
|
vector.add_argument("--run-id", default="")
|
|
vector.add_argument("--date-start", default="")
|
|
vector.add_argument("--date-end", default="")
|
|
vector.add_argument("--reference-date", default="")
|
|
vector.add_argument("--admin-province", default="")
|
|
vector.add_argument("--admin-city", default="")
|
|
vector.set_defaults(func=run_export_points_geojson)
|
|
|
|
monitor = subparsers.add_parser("monitor-points")
|
|
monitor.add_argument("--monitor-config", required=True)
|
|
monitor.add_argument("--timeseries-dir", required=True)
|
|
monitor.add_argument("--export-dir", required=True)
|
|
monitor.add_argument("--point-dir", required=True)
|
|
monitor.add_argument("--mli-par", required=True)
|
|
monitor.add_argument("--slc-par", required=True)
|
|
monitor.add_argument("--dem-par", required=True)
|
|
monitor.add_argument("--lookup", required=True)
|
|
monitor.add_argument("--dates", default="")
|
|
monitor.add_argument("--reference-date", default="")
|
|
monitor.add_argument("--summary-path", required=True)
|
|
monitor.set_defaults(func=run_monitor_points)
|
|
return parser
|
|
|
|
|
|
def main() -> int:
|
|
parser = build_parser()
|
|
args = parser.parse_args()
|
|
return int(args.func(args))
|
|
|
|
|
|
if __name__ == "__main__":
|
|
raise SystemExit(main())
|