From 2d3226bef4b1bb6127e65545673fefcd9152a0f8 Mon Sep 17 00:00:00 2001 From: Forrest Williams Date: Thu, 15 May 2025 11:21:47 -0500 Subject: [PATCH 01/11] add absolute location error assessment --- ale/ale.py | 309 ++++++++++++++++++++++++++++++++++++++++++++ ale/environment.yml | 20 +++ 2 files changed, 329 insertions(+) create mode 100644 ale/ale.py create mode 100644 ale/environment.yml diff --git a/ale/ale.py b/ale/ale.py new file mode 100644 index 0000000..c756b4c --- /dev/null +++ b/ale/ale.py @@ -0,0 +1,309 @@ +from argparse import ArgumentParser +from datetime import datetime +from pathlib import Path + +import isce3 +import matplotlib.pyplot as plt +import numpy as np +import pandas as pd +import requests +from lmfit import Model +from osgeo import gdal, osr +from pyproj import Transformer +from shapely import box +from shapely.geometry import Point +from shapely.ops import transform + + +gdal.UseExceptions() + + +def gaussfit(x, y, A, x0, y0, sigma_x, sigma_y, theta): + theta = np.radians(theta) + sigx2 = sigma_x**2 + sigy2 = sigma_y**2 + a = np.cos(theta) ** 2 / (2 * sigx2) + np.sin(theta) ** 2 / (2 * sigy2) + b = np.sin(theta) ** 2 / (2 * sigx2) + np.cos(theta) ** 2 / (2 * sigy2) + c = np.sin(2 * theta) / (4 * sigx2) - np.sin(2 * theta) / (4 * sigy2) + + expo = -a * (x - x0) ** 2 - b * (y - y0) ** 2 - 2 * c * (x - x0) * (y - y0) + return A * np.exp(expo) + + +def get_cr_df(bounds, epsg, date, outdir): + rosamond_bounds = [-124.409591, 32.534156, -114.131211, 42.009518] + rosamond_bounds = box(*rosamond_bounds) + transformer = Transformer.from_crs('EPSG:4326', f'EPSG:{epsg}', always_xy=True) + rosamond_bounds = transform(transformer.transform, rosamond_bounds) + assert bounds.intersects(rosamond_bounds), f'Images does not intersect with Rosamond bounds {rosamond_bounds}.' + date_str = date.strftime('%Y-%m-%d+%H\u0021%M') + + crdata = outdir / f'{date_str.split("+")[0]}_crdata.csv' + if not crdata.exists(): + res = requests.get( + f'https://uavsar.jpl.nasa.gov/cgi-bin/corner-reflectors.pl?date={str(date_str)}&project=rosamond_plate_location' + ) + crdata.write_bytes(res.content) + cr_df = pd.read_csv(crdata) + new_cols = { + ' "Corner ID"': 'ID', + 'Latitude (deg)': 'lat', + 'Longitude (deg)': 'lon', + 'Azimuth (deg)': 'azm', + 'Height Above Ellipsoid (m)': 'hgt', + 'Side Length (m)': 'slen', + } + cr_df.rename(columns=new_cols, inplace=True) + cr_df.drop(columns=cr_df.keys()[-1], inplace=True) + return cr_df + + +def add_image_location(cr_df, epsg, x_start, y_start, x_spacing, y_spacing, bounds): + blank = [np.nan] * cr_df.shape[0] + blank_bool = [False] * cr_df.shape[0] + cr_df = cr_df.assign( + UTMx=blank, + UTMy=blank, + xloc=blank, + yloc=blank, + xloc_floats=blank, + yloc_floats=blank, + inPoly=blank_bool, + ) + transformer = Transformer.from_crs('EPSG:4326', f'EPSG:{epsg}', always_xy=True) + for idx, row in cr_df.iterrows(): + row['UTMx'], row['UTMy'] = transformer.transform(row['lon'], row['lat']) + row['xloc_floats'] = (row['UTMx'] - x_start) / x_spacing + row['xloc'] = int(round(row['xloc_floats'])) + row['yloc_floats'] = (row['UTMy'] - y_start) / y_spacing + row['yloc'] = int(round(row['yloc_floats'])) + row['inPoly'] = bounds.intersects(Point(row['UTMx'], row['UTMy'])) + cr_df.iloc[idx] = row + + cr_df = cr_df[cr_df['inPoly']] + cr_df.drop('inPoly', axis=1, inplace=True) + cr_df = cr_df.reset_index(drop=True) + return cr_df + + +def filter_valid_data(cr_df, data): + cr_df = cr_df.assign(has_data=False) + for idx, row in cr_df.iterrows(): + xloc = int(row['xloc']) + yloc = int(row['yloc']) + local_data = data[yloc - 2 : yloc + 2, xloc - 3 : xloc + 3] + row['has_data'] = bool(~np.all(np.isnan(local_data))) + cr_df.iloc[idx] = row + cr_df = cr_df[cr_df['has_data']] + cr_df.drop('has_data', axis=1, inplace=True) + cr_df = cr_df.reset_index(drop=True) + return cr_df + + +def filter_orientation(cr_df, direction): + # selecting CRs according to orbit direction + if direction == 'DESC': + # only east-looking CRs (for right-looking descending) + cr_df = cr_df[cr_df['azm'] > 340].reset_index(drop=True) + elif direction == 'ASC': + # only west-looking CRs (for right-looking ascending) + cr_df = cr_df[cr_df['azm'] < 200].reset_index(drop=True) + else: + raise ValueError(f'Invalid direction: {direction}. Use "ASC" or "DESC".') + cr_df = cr_df.loc[cr_df['slen'] > 0.8].reset_index(drop=True) # excluding SWOT CRs (0.7 m as a side length) + return cr_df + + +def plot_crs_on_image(cr_df, data, outdir, fileprefix): + buffer = 50 + min_x = cr_df['xloc'].min() - buffer + max_x = cr_df['xloc'].max() + buffer + min_y = cr_df['yloc'].min() - buffer + max_y = cr_df['yloc'].max() + buffer + + fig, ax = plt.subplots(figsize=(15, 7)) + ax.imshow(data, cmap='gray', interpolation='bilinear', vmin=0.3, vmax=1.7, origin='upper') + ax.set_xlim(min_x, max_x) + ax.set_ylim(min_y, max_y) + ax.axis('off') + + for sl in pd.unique(cr_df.slen): + xx = cr_df.loc[cr_df['slen'] == sl]['xloc'] + yy = cr_df.loc[cr_df['slen'] == sl]['yloc'] + id = cr_df.loc[cr_df['slen'] == sl]['ID'] + color = {2.4384: 'blue', 4.8: 'red', 2.8: 'yellow'}.get(sl, 'green') + ax.scatter(xx, yy, color=color, marker='o', facecolor='none', lw=1) + for id_ann, xx_ann, yy_ann in zip(id, xx, yy): + id_ann = f'{int(id_ann)} ({sl:.1f} m)' + ax.annotate(id_ann, (xx_ann, yy_ann), fontsize=10, color='grey') + + ax.set_aspect(1) + ax.set_title('Corner Reflector Locations') + plt.gca().invert_yaxis() + fig.savefig(outdir / f'{fileprefix}_CR_locations.png', dpi=300, bbox_inches='tight') + + +def calculate_ale_for_cr(point, data, outdir, fileprefix, search_window=100, oversample_factor=4): + ybuff = search_window // 2 + xbuff = search_window // 2 + yrange = np.s_[int(point['yloc'] - ybuff) : int(point['yloc'] + ybuff)] + xrange = np.s_[int(point['xloc'] - xbuff) : int(point['xloc'] + xbuff)] + cropped_data = data[yrange, xrange] + yind, xind = np.unravel_index(np.argmax(cropped_data, axis=None), cropped_data.shape) + + center_buff = 32 + yind_full = int(point['yloc'] - ybuff) + yind + xind_full = int(point['xloc'] - xbuff) + xind + ycenter = np.s_[int(yind_full - center_buff) : int(yind_full + center_buff)] + xcenter = np.s_[int(xind_full - center_buff) : int(xind_full + center_buff)] + centered_data = data[ycenter, xcenter] + + data_ovs = isce3.cal.point_target_info.oversample(centered_data, oversample_factor, baseband=True) + + yoff2 = int(data_ovs.shape[0] / 2) + xoff2 = int(data_ovs.shape[1] / 2) + numpix = 8 + zoom_half_size = numpix * oversample_factor + data_ovs_zoom = data_ovs[ + yoff2 - zoom_half_size : yoff2 + zoom_half_size, xoff2 - zoom_half_size : xoff2 + zoom_half_size + ] + + N = numpix * 2 * oversample_factor + x = np.linspace(0, numpix * 2 * oversample_factor - 1, N) + y = np.linspace(0, numpix * 2 * oversample_factor - 1, N) + Xg, Yg = np.meshgrid(x, y) + fmodel = Model(gaussfit, independent_vars=('x', 'y')) + theta = 0.1 # deg + x0 = numpix * oversample_factor + y0 = numpix * oversample_factor + sigx = 2 + sigy = 5 + A = np.max(data_ovs_zoom) + result = fmodel.fit(data_ovs_zoom, x=Xg, y=Yg, A=A, x0=x0, y0=y0, sigma_x=sigx, sigma_y=sigy, theta=theta) + fit = fmodel.func(Xg, Yg, **result.best_values) + + ypeak_ovs = result.best_values['y0'] + yoff2 - zoom_half_size + ypeak_centered = ypeak_ovs / oversample_factor + ypeak = ypeak_centered + yind_full - center_buff + point['yloc_cr'] = ypeak + + xpeak_ovs = result.best_values['x0'] + xoff2 - zoom_half_size + xpeak_centered = xpeak_ovs / oversample_factor + xpeak = xpeak_centered + xind_full - center_buff + point['xloc_cr'] = xpeak + + xreal_centered = point['xloc'] - xind_full + center_buff + xreal_ovs = xreal_centered * oversample_factor + + yreal_centered = point['yloc'] - yind_full + center_buff + yreal_ovs = yreal_centered * oversample_factor + + plt.rcParams.update({'font.size': 14}) + fig, ax = plt.subplots(1, 3, figsize=(15, 7)) + ax[0].imshow(centered_data, cmap='gray', interpolation=None, origin='upper') + ax[0].plot(xpeak_centered, ypeak_centered, 'r+', label='Return Peak') + ax[0].plot(xreal_centered, yreal_centered, 'b+', label='CR Location') + ax[0].legend() + ax[0].set_title(f'Corner Reflector ID: {int(point["ID"])}') + ax[1].imshow(data_ovs, cmap='gray', interpolation=None, origin='upper') + ax[1].plot(xpeak_ovs, ypeak_ovs, 'r+') + ax[1].plot(xreal_ovs, yreal_ovs, 'b+') + ax[1].set_title(f'Oversampled Corner Reflector ID: {point["ID"]}') + ax[2].imshow(fit, cmap='gray', interpolation=None, origin='upper') + ax[2].plot(result.best_values['x0'], result.best_values['y0'], 'r+') + ax[2].set_title(f'Gaussian Fit Corner Reflector ID: {int(point["ID"])}') + [axi.axis('off') for axi in ax] + fig.tight_layout() + fig.savefig(outdir / f'{fileprefix}_CR_{int(point["ID"])}.png', dpi=300, bbox_inches='tight') + + return point + + +def cr_mean(data): + return np.round(np.nanmean(data), 3) + + +def cr_spread(data): + return np.round(np.nanstd(data) / np.sqrt(np.size(data)), 3) + + +def plot_ale(cr_df, outdir, fileprefix): + east_ale = cr_df['easting_ale'] + north_ale = cr_df['northing_ale'] + ale = cr_df['ale'] + fig, ax = plt.subplots(figsize=(8, 8)) + ax.scatter(east_ale, north_ale, s=20, c='k', alpha=0.6, marker='o') + ax.grid(True) + ax.set_xlim(-15.25, 15.25) + ax.set_ylim(-15.25, 15.25) + ax.axhline(0, color='black') + ax.axvline(0, color='black') + east_metric = f'Easting: {cr_mean(east_ale)} +/- {cr_spread(east_ale)} m' + north_metric = f'Northing: {cr_mean(north_ale)} +/- {cr_spread(north_ale)} m' + overall_metric = f'Overall: {cr_mean(ale)} +/- {cr_spread(ale)} m' + ax.set_title(f'{east_metric}, {north_metric}, {overall_metric}', fontsize=10) + ax.set_xlabel('Easting Error (m)') + ax.set_ylabel('Northing Error (m)') + fig.suptitle('Absolute Location Error') + plt.savefig(outdir / f'{fileprefix}_ale.png', dpi=300, bbox_inches='tight', transparent=True) + + +def ale(filepath, date, direction, outdir, fileprefix): + outdir.mkdir(parents=True, exist_ok=True) + ds = gdal.Open(str(filepath)) + data = ds.GetRasterBand(1).ReadAsArray() + geotransform = ds.GetGeoTransform() + x_start = geotransform[0] + 0.5 * geotransform[1] + y_start = geotransform[3] + 0.5 * geotransform[5] + x_end = x_start + geotransform[1] * ds.RasterXSize + y_end = y_start + geotransform[5] * ds.RasterYSize + bounds = (x_start, y_start, x_end, y_end) + bounds = box(*bounds) + x_spacing = geotransform[1] + y_spacing = geotransform[5] + + srs = osr.SpatialReference() + srs.ImportFromWkt(ds.GetProjectionRef()) + epsg = int(srs.GetAuthorityCode(None)) + cr_df = get_cr_df(bounds, epsg, date, outdir) + cr_df = add_image_location(cr_df, epsg, x_start, y_start, x_spacing, y_spacing, bounds) + cr_df = filter_valid_data(cr_df, data) + cr_df = filter_orientation(cr_df, direction) + cr_df = cr_df.assign(yloc_cr=np.nan, xloc_cr=np.nan) + plot_crs_on_image(cr_df, data, outdir, fileprefix) + for idx, cr in cr_df.iterrows(): + cr = calculate_ale_for_cr(cr, data, outdir, fileprefix) + cr_df.iloc[idx] = cr + + cr_df['easting_ale'] = (cr_df['xloc_cr'] - cr_df['xloc_floats']) * x_spacing + cr_df['northing_ale'] = (cr_df['yloc_cr'] - cr_df['yloc_floats']) * y_spacing + cr_df['ale'] = np.sqrt(cr_df['northing_ale'] ** 2 + cr_df['easting_ale'] ** 2) + cr_df.to_csv(outdir / (fileprefix + '_ale.csv'), index=False) + plot_ale(cr_df, outdir, fileprefix) + + +def main(): + """Example: + python ale.py rtc.tif 2024-12-03 DESC out --fileprefix rtc + """ + parser = ArgumentParser(description='Absolute Location Error Estimation.') + parser.add_argument('filepath', type=str, help='Path to the file to be processed') + parser.add_argument('date', type=str, help='Date of the image collection (YYYY-MM-DD)') + parser.add_argument('direction', choices=['ASC', 'DESC'], help='Direction of the image collection') + parser.add_argument('outdir', type=str, help='Directory to save the results') + parser.add_argument( + '--fileprefix', type=str, help='Prefix for the output filenames (default: input filename)', default=None + ) + args = parser.parse_args() + args.filepath = Path(args.filepath) + args.date = datetime.strptime(args.date, '%Y-%m-%d') + args.outdir = Path(args.outdir) + if args.fileprefix is None: + args.fileprefix = Path(args.filepath).stem + assert args.filepath.exists(), f'File {args.filepath} does not exist.' + + ale(args.filepath, args.date, args.direction, args.outdir, args.fileprefix) + + +if __name__ == '__main__': + main() diff --git a/ale/environment.yml b/ale/environment.yml new file mode 100644 index 0000000..9650d48 --- /dev/null +++ b/ale/environment.yml @@ -0,0 +1,20 @@ +name: multirtc-acc +channels: + - conda-forge + - nodefaults +dependencies: + # For running + - python + - pip + - gdal>=3.0 + - numpy>=1.20 + - isce3 + - requests + - lxml + - shapely + - pyproj + - pandas + - matplotlib + - lmfit + # For static analysis + - ruff From 48f2d56712578026f87458b3c9caea8fadab7fae Mon Sep 17 00:00:00 2001 From: Forrest Williams Date: Tue, 20 May 2025 10:04:02 -0500 Subject: [PATCH 02/11] update for LOS --- ale/ale.py | 40 ++++++++++++++++++++++---------------- src/multirtc/prep_umbra.py | 17 ++++++++++++++++ 2 files changed, 40 insertions(+), 17 deletions(-) diff --git a/ale/ale.py b/ale/ale.py index c756b4c..20d4f0f 100644 --- a/ale/ale.py +++ b/ale/ale.py @@ -97,20 +97,17 @@ def filter_valid_data(cr_df, data): cr_df = cr_df[cr_df['has_data']] cr_df.drop('has_data', axis=1, inplace=True) cr_df = cr_df.reset_index(drop=True) + cr_df = cr_df.loc[cr_df['slen'] > 0.8].reset_index(drop=True) # excluding SWOT CRs (0.7 m as a side length) return cr_df -def filter_orientation(cr_df, direction): - # selecting CRs according to orbit direction - if direction == 'DESC': - # only east-looking CRs (for right-looking descending) - cr_df = cr_df[cr_df['azm'] > 340].reset_index(drop=True) - elif direction == 'ASC': - # only west-looking CRs (for right-looking ascending) - cr_df = cr_df[cr_df['azm'] < 200].reset_index(drop=True) - else: - raise ValueError(f'Invalid direction: {direction}. Use "ASC" or "DESC".') - cr_df = cr_df.loc[cr_df['slen'] > 0.8].reset_index(drop=True) # excluding SWOT CRs (0.7 m as a side length) +def filter_orientation(cr_df, azimuth_angle, margin=60): + valid = [] + for i, row in cr_df.iterrows(): + angle_diff = row['azm'] - azimuth_angle + angle_diff = angle_diff + 360 if angle_diff < -180 else np.abs(angle_diff) + valid.append(angle_diff <= margin) + cr_df = cr_df[np.array(valid)].reset_index(drop=True) return cr_df @@ -227,12 +224,20 @@ def cr_spread(data): return np.round(np.nanstd(data) / np.sqrt(np.size(data)), 3) -def plot_ale(cr_df, outdir, fileprefix): +def plot_ale(cr_df, azmangle, outdir, fileprefix): east_ale = cr_df['easting_ale'] north_ale = cr_df['northing_ale'] ale = cr_df['ale'] + los = np.deg2rad(90 - azmangle + 180) fig, ax = plt.subplots(figsize=(8, 8)) ax.scatter(east_ale, north_ale, s=20, c='k', alpha=0.6, marker='o') + ax.annotate( + 'LOS', + xytext=(np.cos(los) * 10, np.sin(los) * 10), + xy=(0, 0), + arrowprops=dict(edgecolor='darkblue', arrowstyle='<-'), + color='darkblue', + ) ax.grid(True) ax.set_xlim(-15.25, 15.25) ax.set_ylim(-15.25, 15.25) @@ -248,7 +253,7 @@ def plot_ale(cr_df, outdir, fileprefix): plt.savefig(outdir / f'{fileprefix}_ale.png', dpi=300, bbox_inches='tight', transparent=True) -def ale(filepath, date, direction, outdir, fileprefix): +def ale(filepath, date, azmangle, outdir, fileprefix): outdir.mkdir(parents=True, exist_ok=True) ds = gdal.Open(str(filepath)) data = ds.GetRasterBand(1).ReadAsArray() @@ -268,7 +273,7 @@ def ale(filepath, date, direction, outdir, fileprefix): cr_df = get_cr_df(bounds, epsg, date, outdir) cr_df = add_image_location(cr_df, epsg, x_start, y_start, x_spacing, y_spacing, bounds) cr_df = filter_valid_data(cr_df, data) - cr_df = filter_orientation(cr_df, direction) + cr_df = filter_orientation(cr_df, azmangle) cr_df = cr_df.assign(yloc_cr=np.nan, xloc_cr=np.nan) plot_crs_on_image(cr_df, data, outdir, fileprefix) for idx, cr in cr_df.iterrows(): @@ -279,7 +284,7 @@ def ale(filepath, date, direction, outdir, fileprefix): cr_df['northing_ale'] = (cr_df['yloc_cr'] - cr_df['yloc_floats']) * y_spacing cr_df['ale'] = np.sqrt(cr_df['northing_ale'] ** 2 + cr_df['easting_ale'] ** 2) cr_df.to_csv(outdir / (fileprefix + '_ale.csv'), index=False) - plot_ale(cr_df, outdir, fileprefix) + plot_ale(cr_df, azmangle, outdir, fileprefix) def main(): @@ -289,7 +294,7 @@ def main(): parser = ArgumentParser(description='Absolute Location Error Estimation.') parser.add_argument('filepath', type=str, help='Path to the file to be processed') parser.add_argument('date', type=str, help='Date of the image collection (YYYY-MM-DD)') - parser.add_argument('direction', choices=['ASC', 'DESC'], help='Direction of the image collection') + parser.add_argument('azmangle', type=int, help='Azimuth angle of the image (clockwise from North in degrees)') parser.add_argument('outdir', type=str, help='Directory to save the results') parser.add_argument( '--fileprefix', type=str, help='Prefix for the output filenames (default: input filename)', default=None @@ -297,12 +302,13 @@ def main(): args = parser.parse_args() args.filepath = Path(args.filepath) args.date = datetime.strptime(args.date, '%Y-%m-%d') + assert 0 <= args.azmangle <= 360, f'Azimuth angle {args.azmangle} is out of range [0, 360].' args.outdir = Path(args.outdir) if args.fileprefix is None: args.fileprefix = Path(args.filepath).stem assert args.filepath.exists(), f'File {args.filepath} does not exist.' - ale(args.filepath, args.date, args.direction, args.outdir, args.fileprefix) + ale(args.filepath, args.date, args.azmangle, args.outdir, args.fileprefix) if __name__ == '__main__': diff --git a/src/multirtc/prep_umbra.py b/src/multirtc/prep_umbra.py index bc0a968..66ef4ba 100644 --- a/src/multirtc/prep_umbra.py +++ b/src/multirtc/prep_umbra.py @@ -32,6 +32,8 @@ class UmbraSICD: wavelength: float polarization: str lookside: str # 'right' or 'left' + look_angle: int + incidence_angle: int shape: tuple[int, int] scp_index: tuple[int, int] scp_time: float @@ -61,6 +63,18 @@ def calculate_orbit(sensing_start, pos_arp, vel_arp): svs.append(isce3.core.StateVector(t, pos, vel_arp)) return isce3.core.Orbit(svs, sensing_start_isce) + @staticmethod + def calculate_look_angle(arp_ecef, scp_ecef): + look_unit = scp_ecef - arp_ecef + look_unit = look_unit / np.linalg.norm(look_unit) # normalize + look_unit_horz = look_unit.copy() + look_unit_horz[2] = 0 + look_unit_horz /= np.linalg.norm(look_unit_horz) + azimuth_rad = np.arctan2(look_unit_horz[1], look_unit_horz[0]) # in radians + azimuth = np.rad2deg(azimuth_rad) + azimuth = (90 - azimuth) % 360 + return int(np.round(azimuth)) + @staticmethod def calculate_range_range_rate_offset(scp_pos, arp_pos, arp_vel, time_coa): arp_minus_scp = arp_pos - scp_pos @@ -127,6 +141,7 @@ def from_sarpy_sicd(cls, sicd, file_path): transform_matrix = cls.calculate_transform_matrix(sicd.PFA, coa_time) beta0_coeff = sicd.Radiometric.BetaZeroSFPoly.Coefs sigma0_coeff = sicd.Radiometric.SigmaZeroSFPoly.Coefs + look_angle = cls.calculate_look_angle(arp_pos, scp_pos) umbra_sicd = cls( id=Path(file_path).with_suffix('').name, file_path=file_path, @@ -135,6 +150,8 @@ def from_sarpy_sicd(cls, sicd, file_path): wavelength=wavelength, polarization=polarization, lookside=lookside, + look_angle=int(look_angle), + incidence_angle=int(np.round(sicd.SCPCOA.IncidenceAng)), scp_index=(sicd.ImageData.SCPPixel.Row, sicd.ImageData.SCPPixel.Col), scp_time=scp_time, scp_pos=scp_pos, From 2fd4d2a761224e65b21a459ef64813564209f06b Mon Sep 17 00:00:00 2001 From: Forrest Williams Date: Wed, 21 May 2025 07:08:46 -0500 Subject: [PATCH 03/11] small changes --- ale/ale.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/ale/ale.py b/ale/ale.py index 20d4f0f..904d719 100644 --- a/ale/ale.py +++ b/ale/ale.py @@ -101,7 +101,7 @@ def filter_valid_data(cr_df, data): return cr_df -def filter_orientation(cr_df, azimuth_angle, margin=60): +def filter_orientation(cr_df, azimuth_angle, margin=135): valid = [] for i, row in cr_df.iterrows(): angle_diff = row['azm'] - azimuth_angle @@ -228,7 +228,7 @@ def plot_ale(cr_df, azmangle, outdir, fileprefix): east_ale = cr_df['easting_ale'] north_ale = cr_df['northing_ale'] ale = cr_df['ale'] - los = np.deg2rad(90 - azmangle + 180) + los = np.deg2rad(90 - azmangle) fig, ax = plt.subplots(figsize=(8, 8)) ax.scatter(east_ale, north_ale, s=20, c='k', alpha=0.6, marker='o') ax.annotate( From 9d9a35778eda35ac05ce15e19c23128346a6311b Mon Sep 17 00:00:00 2001 From: Forrest Williams Date: Wed, 21 May 2025 07:09:33 -0500 Subject: [PATCH 04/11] update ignore --- .gitignore | 3 +++ 1 file changed, 3 insertions(+) diff --git a/.gitignore b/.gitignore index bbca664..026ca37 100644 --- a/.gitignore +++ b/.gitignore @@ -240,3 +240,6 @@ Sessionx.vim tags # Persistent undo [._]*.un~ + +# Data +ale/umbra* From d0ab18ce728139c88c4bc60211f6b7ab4c9cb8a6 Mon Sep 17 00:00:00 2001 From: Forrest Williams Date: Tue, 27 May 2025 15:37:40 -0500 Subject: [PATCH 05/11] switch los definition --- src/multirtc/prep_umbra.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/multirtc/prep_umbra.py b/src/multirtc/prep_umbra.py index 66ef4ba..cfbbc08 100644 --- a/src/multirtc/prep_umbra.py +++ b/src/multirtc/prep_umbra.py @@ -68,11 +68,12 @@ def calculate_look_angle(arp_ecef, scp_ecef): look_unit = scp_ecef - arp_ecef look_unit = look_unit / np.linalg.norm(look_unit) # normalize look_unit_horz = look_unit.copy() - look_unit_horz[2] = 0 + look_unit_horz[2] = 0.0 look_unit_horz /= np.linalg.norm(look_unit_horz) - azimuth_rad = np.arctan2(look_unit_horz[1], look_unit_horz[0]) # in radians + azimuth_rad = np.arctan2(look_unit_horz[0], look_unit_horz[1]) # in radians azimuth = np.rad2deg(azimuth_rad) - azimuth = (90 - azimuth) % 360 + azimuth = 90 - azimuth + azimuth = 360 + azimuth if azimuth < 0 else azimuth return int(np.round(azimuth)) @staticmethod From 4bc68777c38f511ff453cb8ebee6fccbfe0ea174 Mon Sep 17 00:00:00 2001 From: Forrest Williams Date: Fri, 30 May 2025 09:04:14 -0500 Subject: [PATCH 06/11] get correct look angle and do ale --- ale/ale.py | 13 +++++---- src/multirtc/define_geogrid.py | 18 ------------- src/multirtc/prep_umbra.py | 49 +++++++++++++++++++++++----------- 3 files changed, 40 insertions(+), 40 deletions(-) diff --git a/ale/ale.py b/ale/ale.py index 904d719..ccb0292 100644 --- a/ale/ale.py +++ b/ale/ale.py @@ -101,13 +101,12 @@ def filter_valid_data(cr_df, data): return cr_df -def filter_orientation(cr_df, azimuth_angle, margin=135): - valid = [] - for i, row in cr_df.iterrows(): - angle_diff = row['azm'] - azimuth_angle - angle_diff = angle_diff + 360 if angle_diff < -180 else np.abs(angle_diff) - valid.append(angle_diff <= margin) - cr_df = cr_df[np.array(valid)].reset_index(drop=True) +def filter_orientation(cr_df, azimuth_angle): + looking_east = azimuth_angle < 180 + if looking_east: + cr_df = cr_df[(cr_df['azm'] < 200) & (cr_df['azm'] > 20)].reset_index(drop=True) + else: + cr_df = cr_df[cr_df['azm'] > 340].reset_index(drop=True) return cr_df diff --git a/src/multirtc/define_geogrid.py b/src/multirtc/define_geogrid.py index 87c85c5..8c09a50 100644 --- a/src/multirtc/define_geogrid.py +++ b/src/multirtc/define_geogrid.py @@ -10,26 +10,8 @@ def get_point_epsg(lat, lon): - """ - Get EPSG code based on latitude and longitude - coordinates of a point - - Parameters - ---------- - lat: float - Latitude coordinate of the point - lon: float - Longitude coordinate of the point - - Returns - ------- - epsg: int - UTM zone, Polar stereographic (North / South) - """ - # "Warp" longitude value into [-180.0, 180.0] if (lon >= 180.0) or (lon <= -180.0): lon = (lon + 180.0) % 360.0 - 180.0 - if lat >= 75.0: epsg = 3413 elif lat <= -75.0: diff --git a/src/multirtc/prep_umbra.py b/src/multirtc/prep_umbra.py index cfbbc08..6f5dd0b 100644 --- a/src/multirtc/prep_umbra.py +++ b/src/multirtc/prep_umbra.py @@ -63,18 +63,33 @@ def calculate_orbit(sensing_start, pos_arp, vel_arp): svs.append(isce3.core.StateVector(t, pos, vel_arp)) return isce3.core.Orbit(svs, sensing_start_isce) - @staticmethod - def calculate_look_angle(arp_ecef, scp_ecef): - look_unit = scp_ecef - arp_ecef - look_unit = look_unit / np.linalg.norm(look_unit) # normalize - look_unit_horz = look_unit.copy() - look_unit_horz[2] = 0.0 - look_unit_horz /= np.linalg.norm(look_unit_horz) - azimuth_rad = np.arctan2(look_unit_horz[0], look_unit_horz[1]) # in radians - azimuth = np.rad2deg(azimuth_rad) - azimuth = 90 - azimuth - azimuth = 360 + azimuth if azimuth < 0 else azimuth - return int(np.round(azimuth)) + def calculate_look_angles(arp_ecef, scp_ecef): + # Convert observer ECEF to geodetic lat/lon/alt + transformer = pyproj.Transformer.from_crs('EPSG:4978', 'EPSG:4979', always_xy=True) + lon_deg, lat_deg, _ = transformer.transform(*scp_ecef) + + # ECEF to ENU rotation matrix centered at the observer + sin_lat = np.sin(np.deg2rad(lat_deg)) + cos_lat = np.cos(np.deg2rad(lat_deg)) + sin_lon = np.sin(np.deg2rad(lon_deg)) + cos_lon = np.cos(np.deg2rad(lon_deg)) + rotation_matrix = np.array( + [ + [-sin_lon, cos_lon, 0], + [-sin_lat * cos_lon, -sin_lat * sin_lon, cos_lat], + [cos_lat * cos_lon, cos_lat * sin_lon, sin_lat], + ] + ) + + # Transform the LOS vector into ENU frame + los = np.subtract(arp_ecef, scp_ecef) + topocentric = rotation_matrix @ los + east, north, up = topocentric + + # Compute azimuth and elevation + azimuth = np.arctan2(east, north) % (2 * np.pi) + elevation = np.arcsin(up / np.linalg.norm(topocentric)) + return np.rad2deg(azimuth), np.rad2deg(elevation) @staticmethod def calculate_range_range_rate_offset(scp_pos, arp_pos, arp_vel, time_coa): @@ -142,7 +157,11 @@ def from_sarpy_sicd(cls, sicd, file_path): transform_matrix = cls.calculate_transform_matrix(sicd.PFA, coa_time) beta0_coeff = sicd.Radiometric.BetaZeroSFPoly.Coefs sigma0_coeff = sicd.Radiometric.SigmaZeroSFPoly.Coefs - look_angle = cls.calculate_look_angle(arp_pos, scp_pos) + center = Point(sicd.GeoData.SCP.LLH.Lon, sicd.GeoData.SCP.LLH.Lat) + azimuth_angle, elevation_angle = cls.calculate_look_angles(arp_pos, scp_pos) + assert np.isclose(azimuth_angle, sicd.SCPCOA.AzimAng), 'Azimuth angle does not match SCPCOA AzimuthAng' + assert np.isclose(elevation_angle, sicd.SCPCOA.GrazeAng), 'Elevation angle does not match SCPCOA ElevationAng' + look_angle = int(azimuth_angle + 180) % 360 umbra_sicd = cls( id=Path(file_path).with_suffix('').name, file_path=file_path, @@ -151,12 +170,12 @@ def from_sarpy_sicd(cls, sicd, file_path): wavelength=wavelength, polarization=polarization, lookside=lookside, - look_angle=int(look_angle), + look_angle=look_angle, incidence_angle=int(np.round(sicd.SCPCOA.IncidenceAng)), scp_index=(sicd.ImageData.SCPPixel.Row, sicd.ImageData.SCPPixel.Col), scp_time=scp_time, scp_pos=scp_pos, - center=Point(sicd.GeoData.SCP.LLH.Lon, sicd.GeoData.SCP.LLH.Lat), + center=center, scp_hae=scp_hae, coa_time=coa_time, arp_pos=arp_pos, From 272bc0829b8ab5ae69122293ec23675f773cf284 Mon Sep 17 00:00:00 2001 From: Forrest Williams Date: Fri, 30 May 2025 09:05:57 -0500 Subject: [PATCH 07/11] update changelog --- CHANGELOG.md | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 28c9aff..89b4db2 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -8,6 +8,11 @@ and uses [Semantic Versioning](https://semver.org/spec/v2.0.0.html). ## [0.2.0] +### Added +* Cal/Val scripts for performing absolute geolocation error (ALE) assessment + +## [0.2.0] + ### Added * Support for Capella SICD SLC products From b8c6716985c63cd5391712d84cba39613bc0211f Mon Sep 17 00:00:00 2001 From: Forrest Williams Date: Fri, 30 May 2025 10:17:08 -0500 Subject: [PATCH 08/11] rough addition of pfa as sicd subclass --- src/multirtc/prep_umbra.py | 9 ++- src/multirtc/sicd.py | 152 ++++++++++++++++++++++++++++++++++++- 2 files changed, 154 insertions(+), 7 deletions(-) diff --git a/src/multirtc/prep_umbra.py b/src/multirtc/prep_umbra.py index 6f5dd0b..89361dd 100644 --- a/src/multirtc/prep_umbra.py +++ b/src/multirtc/prep_umbra.py @@ -12,6 +12,7 @@ from multirtc import dem from multirtc.define_geogrid import get_point_epsg +from multirtc.sicd import SicdPfaSlc def check_poly_order(poly): @@ -313,9 +314,11 @@ def prep_umbra(granule_path: Path, work_dir: Optional[Path] = None) -> Path: """ if work_dir is None: work_dir = Path.cwd() - reader = SICDReader(str(granule_path)) - sicd = reader.get_sicds_as_tuple()[0] - umbra_sicd = UmbraSICD.from_sarpy_sicd(sicd, granule_path) + umbra_sicd = SicdPfaSlc(granule_path) + # reader = SICDReader(str(granule_path)) + # sicd = reader.get_sicds_as_tuple()[0] + # umbra_sicd = UmbraSICD.from_sarpy_sicd(sicd, granule_path) + breakpoint() dem_path = work_dir / 'dem.tif' dem.download_opera_dem_for_footprint(dem_path, umbra_sicd.footprint) diff --git a/src/multirtc/sicd.py b/src/multirtc/sicd.py index 685e45e..dae6e8c 100644 --- a/src/multirtc/sicd.py +++ b/src/multirtc/sicd.py @@ -3,6 +3,7 @@ import isce3 import numpy as np +import pyproj from numpy.polynomial.polynomial import polyval2d from osgeo import gdal from sarpy.io.complex.sicd import SICDReader @@ -44,12 +45,47 @@ def __init__(self, sicd_path): + sicd.Grid.Row.UVectECF.get_array() * (0 - self.shift[0]) * self.spacing[0] ) self.starting_range = np.linalg.norm(sicd.SCPCOA.ARPPos.get_array() - starting_row_pos) - last_line_time = sicd.Grid.TimeCOAPoly(0, self.shape[1] - self.shift[1]) - first_line_time = sicd.Grid.TimeCOAPoly(0, -self.shift[1]) - self.az_reversed = last_line_time >= first_line_time + self.raw_time_coa_poly = sicd.Grid.TimeCOAPoly + last_line_time = self.raw_time_coa_poly(0, self.shape[1] - self.shift[1]) + first_line_time = self.raw_time_coa_poly(0, -self.shift[1]) + self.az_reversed = last_line_time > first_line_time + self.arp_pos = sicd.SCPCOA.ARPPos.get_array() + self.scp_pos = sicd.GeoData.SCP.ECF.get_array() + azimuth_angle, elevation_angle = self.calculate_look_angles() + assert np.isclose(azimuth_angle, sicd.SCPCOA.AzimAng), 'Azimuth angle does not match SCPCOA AzimuthAng' + assert np.isclose(elevation_angle, sicd.SCPCOA.GrazeAng), 'Elevation angle does not match SCPCOA ElevationAng' + self.look_angle = int(azimuth_angle + 180) % 360 self.beta0 = sicd.Radiometric.BetaZeroSFPoly self.sigma0 = sicd.Radiometric.SigmaZeroSFPoly + def calculate_look_angles(self): + # Convert observer ECEF to geodetic lat/lon/alt + transformer = pyproj.Transformer.from_crs('EPSG:4978', 'EPSG:4979', always_xy=True) + lon_deg, lat_deg, _ = transformer.transform(*self.scp_pos) + + # ECEF to ENU rotation matrix centered at the observer + sin_lat = np.sin(np.deg2rad(lat_deg)) + cos_lat = np.cos(np.deg2rad(lat_deg)) + sin_lon = np.sin(np.deg2rad(lon_deg)) + cos_lon = np.cos(np.deg2rad(lon_deg)) + rotation_matrix = np.array( + [ + [-sin_lon, cos_lon, 0], + [-sin_lat * cos_lon, -sin_lat * sin_lon, cos_lat], + [cos_lat * cos_lon, cos_lat * sin_lon, sin_lat], + ] + ) + + # Transform the LOS vector into ENU frame + los = np.subtract(self.arp_pos, self.scp_pos) + topocentric = rotation_matrix @ los + east, north, up = topocentric + + # Compute azimuth and elevation + azimuth = np.arctan2(east, north) % (2 * np.pi) + elevation = np.arcsin(up / np.linalg.norm(topocentric)) + return np.rad2deg(azimuth), np.rad2deg(elevation) + def get_xrow_ycol(self) -> np.ndarray: """Calculate xrow and ycol SICD.""" irow = np.tile(np.arange(self.shape[0]), (self.shape[1], 1)).T @@ -119,7 +155,7 @@ def create_complex_beta0(self, outpath, isce_format=True): class SicdRzdSlc(SlcTemplate, SicdSlc): def __init__(self, sicd_path): super().__init__(sicd_path) - assert self.source.Grid.Type == 'RGZERO', 'Only range zero doppler grids supported for Capella data' + assert self.source.Grid.Type == 'RGZERO', 'Only range zero doppler grids are supported for by this class' first_col_time = self.source.RMA.INCA.TimeCAPoly(0 - self.shift[1]) last_col_time = self.source.RMA.INCA.TimeCAPoly(self.shape[1] - self.shift[1]) self.sensing_start = min(first_col_time, last_col_time) @@ -154,3 +190,111 @@ def get_radar_grid(self): ref_epoch=to_isce_datetime(self.reference_time), ) return radar_grid + + +class SicdPfaSlc(SlcTemplate, SicdSlc): + def __init__(self, sicd_path): + super().__init__(sicd_path) + assert self.source.ImageFormation.ImageFormAlgo == 'PFA', 'Only PFA-focused data are supported by this class' + assert self.source.Grid.Type == 'RGAZIM', 'Only range azimuth grids are supported by this class' + assert self.raw_time_coa_poly.Coefs.size == 1, 'Only constant COA time is currently supported' + self.coa_time = self.raw_time_coa_poly.Coefs[0][0] + self.arp_vel = self.source.SCPCOA.ARPVel.get_array() + self.scp_time = self.reference_time + timedelta(self.source.SCPCOA.SCPTime) + self.sensing_start = self.coa_time + self.pfa_vars = self.source.PFA + self.orbit = self.get_orbit() + self.rrdot_offset = self.calculate_range_range_rate_offset() + self.transform_matrix = self.calculate_transform_matrix() + self.transform_matrix_inv = np.linalg.inv(self.transform_matrix) + # Without ISCE3 support for PFA grids, these properties are undefined + self.radar_grid = None + self.doppler_centroid_grid = None + self.prf = np.nan + + def get_orbit(self): + svs = [] + sensing_start_isce = to_isce_datetime(self.scp_time) + for offset_sec in range(-10, 10): + t = self.scp_time + timedelta(offset_sec) + t_isce = to_isce_datetime(t) + pos = self.arp_vel * offset_sec + self.arp_pos + svs.append(isce3.core.StateVector(t_isce, pos, self.arp_vel)) + return isce3.core.Orbit(svs, sensing_start_isce) + + def calculate_range_range_rate_offset(self): + arp_minus_scp = self.arp_pos - self.scp_pos + range_scp_to_coa = np.linalg.norm(arp_minus_scp, axis=-1) + range_rate_scp_to_coa = np.sum(self.arp_vel * arp_minus_scp, axis=-1) / range_scp_to_coa + rrdot_offset = np.array([range_scp_to_coa, range_rate_scp_to_coa]) + return rrdot_offset + + def calculate_transform_matrix(self): + polar_ang_poly = self.pfa_vars.PolarAngPoly + spatial_freq_sf_poly = self.pfa_vars.SpatialFreqSFPoly + polar_ang_poly_der = polar_ang_poly.derivative(der_order=1, return_poly=True) + spatial_freq_sf_poly_der = spatial_freq_sf_poly.derivative(der_order=1, return_poly=True) + + polar_ang_poly_der = polar_ang_poly.derivative(der_order=1, return_poly=True) + spatial_freq_sf_poly_der = spatial_freq_sf_poly.derivative(der_order=1, return_poly=True) + + thetaTgtCoa = polar_ang_poly(self.coa_time) + dThetaDtTgtCoa = polar_ang_poly_der(self.coa_time) + + # Compute polar aperture scale factor (KSF) and derivative + # wrt polar angle + ksfTgtCoa = spatial_freq_sf_poly(thetaTgtCoa) + dKsfDThetaTgtCoa = spatial_freq_sf_poly_der(thetaTgtCoa) + + # Compute spatial frequency domain phase slopes in Ka and Kc directions + # NB: sign for the phase may be ignored as it is cancelled + # in a subsequent computation. + dPhiDKaTgtCoa = np.array([np.cos(thetaTgtCoa), np.sin(thetaTgtCoa)]) + dPhiDKcTgtCoa = np.array([-np.sin(thetaTgtCoa), np.cos(thetaTgtCoa)]) + + transform_matrix = np.zeros((2, 2)) + transform_matrix[0, :] = ksfTgtCoa * dPhiDKaTgtCoa + transform_matrix[1, :] = dThetaDtTgtCoa * (dKsfDThetaTgtCoa * dPhiDKaTgtCoa + ksfTgtCoa * dPhiDKcTgtCoa) + return transform_matrix + + def rowcol2geo(self, rc: np.ndarray, hae: float) -> np.ndarray: + """Transform grid (row, col) coordinates to ECEF coordinates. + + Args: + rc: 2D array of (row, col) coordinates + hae: Height above ellipsoid (meters) + + Returns: + np.ndarray: ECEF coordinates + """ + dem = isce3.geometry.DEMInterpolator(hae) + elp = isce3.core.Ellipsoid() + rgaz = (rc - np.array(self.shift)[None, :]) * np.array(self.spacing)[None, :] + rrdot = np.dot(self.transform_matrix, rgaz.T) + self.rrdot_offset[:, None] + side = isce3.core.LookSide(1) if self.lookside == 'left' else isce3.core.LookSide(-1) + pts_ecf = [] + wvl = 1.0 + for pt in rrdot.T: + r = pt[0] + dop = -pt[1] * 2 / wvl + llh = isce3.geometry.rdr2geo(0.0, r, self.orbit, side, dop, wvl, dem, threshold=1.0e-8, maxiter=50) + pts_ecf.append(elp.lon_lat_to_xyz(llh)) + return np.vstack(pts_ecf) + + def geo2rowcol(self, xyz: np.ndarray) -> np.ndarray: + """Transform ECEF xyz to (row, col). + + Args: + xyz: ECEF coordinates + + Returns: + (row, col) coordinates + """ + rrdot = np.zeros((2, xyz.shape[0])) + rrdot[0, :] = np.linalg.norm(xyz - self.arp_pos[None, :], axis=1) + rrdot[1, :] = np.dot(-self.arp_vel, (xyz - self.arp_pos[None, :]).T) / rrdot[0, :] + rgaz = np.dot(self.transform_matrix_inv, (rrdot - self.rrdot_offset[:, None])) + rgaz /= np.array(self.spacing)[:, None] + rgaz += np.array(self.shift)[:, None] + row_col = rgaz.T.copy() + return row_col From d71333b90dd0f64810284cc885873e63b50918c8 Mon Sep 17 00:00:00 2001 From: Forrest Williams Date: Fri, 30 May 2025 14:57:34 -0500 Subject: [PATCH 09/11] successful umbra run --- src/multirtc/create_rtc.py | 22 ++++++++++------------ src/multirtc/multirtc.py | 13 +++++++------ src/multirtc/prep_umbra.py | 5 ----- src/multirtc/sicd.py | 38 +++++++++++++++++++++++++++++++++++--- 4 files changed, 52 insertions(+), 26 deletions(-) diff --git a/src/multirtc/create_rtc.py b/src/multirtc/create_rtc.py index b165c2e..6f9c45e 100755 --- a/src/multirtc/create_rtc.py +++ b/src/multirtc/create_rtc.py @@ -516,44 +516,42 @@ def rtc(slc, geogrid, opts): logger.info(f'elapsed time: {t_end - t_start}') -def umbra_rtc(umbra_sicd, geogrid, dem_path, output_dir): +def pfa_prototype_geocode(sicd, geogrid, dem_path, output_dir): interp_method = isce3.core.DataInterpMethod.BIQUINTIC - sigma0_data = umbra_sicd.load_corrected_data('sigma0') + sigma0_data = sicd.load_scaled_data('sigma0', power=True) slc_lut = isce3.core.LUT2d( np.arange(sigma0_data.shape[1]), np.arange(sigma0_data.shape[0]), sigma0_data, interp_method ) - ecef = pyproj.CRS(4978) # ECEF on WGS84 Ellipsoid - lla = pyproj.CRS(4979) # WGS84 - lla2ecef = pyproj.Transformer.from_crs(lla, ecef, always_xy=True) + local2ecef = pyproj.Transformer.from_crs(f'EPSG:{geogrid.epsg}', 'EPSG:4978', always_xy=True) + local2ll = pyproj.Transformer.from_crs(f'EPSG:{geogrid.epsg}', 'EPSG:4326', always_xy=True) dem_raster = isce3.io.Raster(str(dem_path)) dem = isce3.geometry.DEMInterpolator() dem.load_dem(dem_raster) dem.interp_method = interp_method - output = np.zeros((geogrid.length, geogrid.width), dtype=np.complex64) + output = np.zeros((geogrid.length, geogrid.width), dtype=np.float32) mask = np.zeros((geogrid.length, geogrid.width), dtype=bool) n_iters = geogrid.width * geogrid.length for i, j in tqdm(itertools.product(range(geogrid.width), range(geogrid.length)), total=n_iters): x = geogrid.start_x + (i * geogrid.spacing_x) y = geogrid.start_y + (j * geogrid.spacing_y) + lon, lat = local2ll.transform(x, y) hae = dem.interpolate_lonlat(np.deg2rad(x), np.deg2rad(y)) # ISCE3 expects lat/lon to be in radians! - ecef_x, ecef_y, ecef_z = lla2ecef.transform(x, y, hae) - row, col = umbra_sicd.geo2rowcol(np.array([(ecef_x, ecef_y, ecef_z)]))[0] + ecef_x, ecef_y, ecef_z = local2ecef.transform(x, y, hae) + row, col = sicd.geo2rowcol(np.array([(ecef_x, ecef_y, ecef_z)]))[0] if slc_lut.contains(row, col): output[j, i] = slc_lut.eval(row, col) mask[j, i] = 1 - output = 10 * np.log10(output) output[mask == 0] = np.nan - - output_path = output_dir / f'{umbra_sicd.id}_{umbra_sicd.polarization}.tif' + output_path = output_dir / f'{sicd.id}_{sicd.polarization}.tif' driver = gdal.GetDriverByName('GTiff') out_ds = driver.Create(str(output_path), geogrid.width, geogrid.length, 1, gdal.GDT_Float32) # account for pixel as area start_x = geogrid.start_x - (geogrid.spacing_x / 2) start_y = geogrid.start_y + (geogrid.spacing_y / 2) out_ds.SetGeoTransform([start_x, geogrid.spacing_x, 0, start_y, 0, geogrid.spacing_y]) - out_ds.SetProjection(pyproj.CRS(4326).to_wkt()) + out_ds.SetProjection(pyproj.CRS(geogrid.epsg).to_wkt()) out_ds.GetRasterBand(1).WriteArray(output) out_ds.GetRasterBand(1).SetNoDataValue(np.nan) out_ds.SetMetadata({'AREA_OR_POINT': 'Area'}) diff --git a/src/multirtc/multirtc.py b/src/multirtc/multirtc.py index 4ae8f3d..c6e45d1 100644 --- a/src/multirtc/multirtc.py +++ b/src/multirtc/multirtc.py @@ -4,7 +4,7 @@ import isce3 -from multirtc.create_rtc import rtc, umbra_rtc +from multirtc.create_rtc import pfa_prototype_geocode, rtc from multirtc.define_geogrid import generate_geogrids from multirtc.prep_burst import prep_burst from multirtc.prep_capella import prep_capella @@ -93,18 +93,19 @@ def opera_rtc_umbra_sicd(granule: str, resolution: int = 30, work_dir: Optional[ raise FileNotFoundError(f'Umbra SICD must be present in input dir {input_dir} for processing.') [d.mkdir(parents=True, exist_ok=True) for d in [input_dir, output_dir]] umbra_sicd, dem_path = prep_umbra(granule_path, work_dir=input_dir) - geogrid = generate_geogrids(umbra_sicd, resolution, rda=False) - umbra_rtc(umbra_sicd, geogrid, dem_path, output_dir=output_dir) + geogrid = umbra_sicd.create_geogrid(spacing_meters=resolution) + pfa_prototype_geocode(umbra_sicd, geogrid, dem_path, output_dir=output_dir) def main(): """Create an OPERA RTC for an Umbra SICD SLC granule Example command: - multirtc umbra_image.ntif --resolution 40 + multirtc UMBRA umbra_image.ntif --resolution 40 """ + supported = ['S1', 'UMBRA', 'CAPELLA'] parser = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.ArgumentDefaultsHelpFormatter) - parser.add_argument('platform', choices=['S1', 'UMBRA', 'CAPELLA'], help='Platform to create RTC for') + parser.add_argument('platform', choices=supported, help='Platform to create RTC for') parser.add_argument('granule', help='Data granule to create an RTC for.') parser.add_argument('--resolution', default=30, type=float, help='Resolution of the output RTC (m)') parser.add_argument('--work-dir', type=Path, default=None, help='Working directory for processing') @@ -117,7 +118,7 @@ def main(): elif args.platform == 'CAPELLA': opera_rtc_capella_sicd(args.granule, args.resolution, args.work_dir) else: - raise NotImplementedError('Only Sentinel-1 burst and Umbra processing are supported at this time') + raise ValueError(f'Unsupported platform {args.platform}. Supported platforms are {", ".join(supported)}.') if __name__ == '__main__': diff --git a/src/multirtc/prep_umbra.py b/src/multirtc/prep_umbra.py index 89361dd..a7bca9b 100644 --- a/src/multirtc/prep_umbra.py +++ b/src/multirtc/prep_umbra.py @@ -315,11 +315,6 @@ def prep_umbra(granule_path: Path, work_dir: Optional[Path] = None) -> Path: if work_dir is None: work_dir = Path.cwd() umbra_sicd = SicdPfaSlc(granule_path) - # reader = SICDReader(str(granule_path)) - # sicd = reader.get_sicds_as_tuple()[0] - # umbra_sicd = UmbraSICD.from_sarpy_sicd(sicd, granule_path) - breakpoint() - dem_path = work_dir / 'dem.tif' dem.download_opera_dem_for_footprint(dem_path, umbra_sicd.footprint) return umbra_sicd, dem_path diff --git a/src/multirtc/sicd.py b/src/multirtc/sicd.py index dae6e8c..b8427ef 100644 --- a/src/multirtc/sicd.py +++ b/src/multirtc/sicd.py @@ -9,6 +9,7 @@ from sarpy.io.complex.sicd import SICDReader from shapely.geometry import Point, Polygon +from multirtc import define_geogrid from multirtc.base import SlcTemplate, to_isce_datetime @@ -25,6 +26,7 @@ def __init__(self, sicd_path): self.filepath = Path(sicd_path) self.footprint = Polygon([(ic.Lon, ic.Lat) for ic in sicd.GeoData.ImageCorners]) self.center = Point(sicd.GeoData.SCP.LLH.Lon, sicd.GeoData.SCP.LLH.Lat) + self.scp_hae = sicd.GeoData.SCP.LLH.HAE self.lookside = 'right' if sicd.SCPCOA.SideOfTrack == 'R' else 'left' center_frequency = sicd.RadarCollection.TxFrequency.Min + sicd.RadarCollection.TxFrequency.Max / 2 @@ -90,15 +92,17 @@ def get_xrow_ycol(self) -> np.ndarray: """Calculate xrow and ycol SICD.""" irow = np.tile(np.arange(self.shape[0]), (self.shape[1], 1)).T irow -= self.scp_index[0] - xrow = irow * self.row_mult + xrow = irow * self.spacing[0] icol = np.tile(np.arange(self.shape[1]), (self.shape[0], 1)) icol -= self.scp_index[1] - ycol = icol * self.col_mult + ycol = icol * self.spacing[1] return xrow, ycol def load_data(self): - return self.source[:, :] + reader = SICDReader(str(self.filepath)) + data = reader[:, :] + return data def load_scaled_data(self, scale, power=False): if scale == 'beta0': @@ -298,3 +302,31 @@ def geo2rowcol(self, xyz: np.ndarray) -> np.ndarray: rgaz += np.array(self.shift)[:, None] row_col = rgaz.T.copy() return row_col + + def create_geogrid(self, spacing_meters): + epsg_local = define_geogrid.get_point_epsg(self.center.y, self.center.x) + ecef = pyproj.CRS(4978) # ECEF on WGS84 Ellipsoid + local = pyproj.CRS(epsg_local) + ecef2local = pyproj.Transformer.from_crs(ecef, local, always_xy=True) + x_spacing = spacing_meters + y_spacing = -1 * spacing_meters + + points = np.array([(0, 0), (0, self.shape[1]), self.shape, (self.shape[0], 0)]) + geos = self.rowcol2geo(points, hae=self.scp_hae) + + points = np.vstack(ecef2local.transform(geos[:, 0], geos[:, 1], geos[:, 2])).T + minx, maxx = np.min(points[:, 0]), np.max(points[:, 0]) + miny, maxy = np.min(points[:, 1]), np.max(points[:, 1]) + + width = (maxx - minx) // x_spacing + length = (maxy - miny) // np.abs(y_spacing) + geogrid = isce3.product.GeoGridParameters( + start_x=float(minx), + start_y=float(maxy), + spacing_x=float(x_spacing), + spacing_y=float(y_spacing), + length=int(length), + width=int(width), + epsg=epsg_local, + ) + return geogrid From 881173f580d6c1543d514266234a30c9c9b7bd30 Mon Sep 17 00:00:00 2001 From: Forrest Williams Date: Fri, 30 May 2025 14:58:23 -0500 Subject: [PATCH 10/11] remove old pfa implementation --- src/multirtc/prep_umbra.py | 300 ------------------------------------- 1 file changed, 300 deletions(-) diff --git a/src/multirtc/prep_umbra.py b/src/multirtc/prep_umbra.py index a7bca9b..0d7450b 100644 --- a/src/multirtc/prep_umbra.py +++ b/src/multirtc/prep_umbra.py @@ -1,310 +1,10 @@ -from dataclasses import dataclass -from datetime import datetime from pathlib import Path from typing import Optional -import isce3 -import numpy as np -import pyproj -from numpy.polynomial.polynomial import polyval2d -from sarpy.io.complex.sicd import SICDReader -from shapely.geometry import Polygon - from multirtc import dem -from multirtc.define_geogrid import get_point_epsg from multirtc.sicd import SicdPfaSlc -def check_poly_order(poly): - assert len(poly.Coefs) == poly.order1 + 1, 'Polynomial order does not match number of coefficients' - - -@dataclass -class Point: - x: float - y: float - - -@dataclass -class UmbraSICD: - id: str - file_path: Path - footprint: Polygon - wavelength: float - polarization: str - lookside: str # 'right' or 'left' - look_angle: int - incidence_angle: int - shape: tuple[int, int] - scp_index: tuple[int, int] - scp_time: float - scp_pos: np.ndarray - center: Point - scp_hae: float - coa_time: float - arp_pos: np.ndarray - arp_vel: np.ndarray - grid_shift: np.ndarray - grid_mult: np.ndarray - rrdot_offset: np.ndarray - transform_matrix: np.ndarray - transform_matrix_inv: np.ndarray - orbit: isce3.core.Orbit - beta0_coeff: np.ndarray - sigma0_coeff: np.ndarray - - @staticmethod - def calculate_orbit(sensing_start, pos_arp, vel_arp): - svs = [] - sensing_start_isce = isce3.core.DateTime(datetime.utcfromtimestamp(sensing_start.astype(int) * 1e-9)) - for offset_sec in range(-5, 6): - t = sensing_start + np.timedelta64(offset_sec, 's') - t = isce3.core.DateTime(datetime.utcfromtimestamp(t.astype(int) * 1e-9)) - pos = vel_arp * offset_sec + pos_arp - svs.append(isce3.core.StateVector(t, pos, vel_arp)) - return isce3.core.Orbit(svs, sensing_start_isce) - - def calculate_look_angles(arp_ecef, scp_ecef): - # Convert observer ECEF to geodetic lat/lon/alt - transformer = pyproj.Transformer.from_crs('EPSG:4978', 'EPSG:4979', always_xy=True) - lon_deg, lat_deg, _ = transformer.transform(*scp_ecef) - - # ECEF to ENU rotation matrix centered at the observer - sin_lat = np.sin(np.deg2rad(lat_deg)) - cos_lat = np.cos(np.deg2rad(lat_deg)) - sin_lon = np.sin(np.deg2rad(lon_deg)) - cos_lon = np.cos(np.deg2rad(lon_deg)) - rotation_matrix = np.array( - [ - [-sin_lon, cos_lon, 0], - [-sin_lat * cos_lon, -sin_lat * sin_lon, cos_lat], - [cos_lat * cos_lon, cos_lat * sin_lon, sin_lat], - ] - ) - - # Transform the LOS vector into ENU frame - los = np.subtract(arp_ecef, scp_ecef) - topocentric = rotation_matrix @ los - east, north, up = topocentric - - # Compute azimuth and elevation - azimuth = np.arctan2(east, north) % (2 * np.pi) - elevation = np.arcsin(up / np.linalg.norm(topocentric)) - return np.rad2deg(azimuth), np.rad2deg(elevation) - - @staticmethod - def calculate_range_range_rate_offset(scp_pos, arp_pos, arp_vel, time_coa): - arp_minus_scp = arp_pos - scp_pos - range_scp_to_coa = np.linalg.norm(arp_minus_scp, axis=-1) - range_rate_scp_to_coa = np.sum(arp_vel * arp_minus_scp, axis=-1) / range_scp_to_coa - rrdot_offset = np.array([range_scp_to_coa, range_rate_scp_to_coa]) - return rrdot_offset - - @staticmethod - def calculate_transform_matrix(pfa, time_coa): - polar_ang_poly = pfa.PolarAngPoly - spatial_freq_sf_poly = pfa.SpatialFreqSFPoly - polar_ang_poly_der = polar_ang_poly.derivative(der_order=1, return_poly=True) - spatial_freq_sf_poly_der = spatial_freq_sf_poly.derivative(der_order=1, return_poly=True) - - polar_ang_poly_der = polar_ang_poly.derivative(der_order=1, return_poly=True) - spatial_freq_sf_poly_der = spatial_freq_sf_poly.derivative(der_order=1, return_poly=True) - - thetaTgtCoa = polar_ang_poly(time_coa) - dThetaDtTgtCoa = polar_ang_poly_der(time_coa) - - # Compute polar aperture scale factor (KSF) and derivative - # wrt polar angle - ksfTgtCoa = spatial_freq_sf_poly(thetaTgtCoa) - dKsfDThetaTgtCoa = spatial_freq_sf_poly_der(thetaTgtCoa) - - # Compute spatial frequency domain phase slopes in Ka and Kc directions - # NB: sign for the phase may be ignored as it is cancelled - # in a subsequent computation. - dPhiDKaTgtCoa = np.array([np.cos(thetaTgtCoa), np.sin(thetaTgtCoa)]) - dPhiDKcTgtCoa = np.array([-np.sin(thetaTgtCoa), np.cos(thetaTgtCoa)]) - - transform_matrix = np.zeros((2, 2)) - transform_matrix[0, :] = ksfTgtCoa * dPhiDKaTgtCoa - transform_matrix[1, :] = dThetaDtTgtCoa * (dKsfDThetaTgtCoa * dPhiDKaTgtCoa + ksfTgtCoa * dPhiDKcTgtCoa) - return transform_matrix - - @classmethod - def from_sarpy_sicd(cls, sicd, file_path): - center_frequency = sicd.RadarCollection.TxFrequency.Min + sicd.RadarCollection.TxFrequency.Max / 2 - wavelength = isce3.core.speed_of_light / center_frequency - polarization = sicd.RadarCollection.RcvChannels[0].TxRcvPolarization.replace(':', '') - lookside = 'right' if sicd.SCPCOA.SideOfTrack == 'R' else 'left' - footprint = Polygon([(ic.Lon, ic.Lat) for ic in sicd.GeoData.ImageCorners]) - - coef_time_coa = sicd.Grid.TimeCOAPoly.Coefs - assert coef_time_coa.size == 1, 'Only constant COA time is currently supported' - coa_time = coef_time_coa[0][0] - - arp_pos = np.array([sicd.SCPCOA.ARPPos.X, sicd.SCPCOA.ARPPos.Y, sicd.SCPCOA.ARPPos.Z]) - arp_vel = np.array([sicd.SCPCOA.ARPVel.X, sicd.SCPCOA.ARPVel.Y, sicd.SCPCOA.ARPVel.Z]) - scp_time = sicd.Timeline.CollectStart + np.timedelta64(int(sicd.SCPCOA.SCPTime * 1e9), 'ns') - orbit = cls.calculate_orbit(scp_time, arp_pos, arp_vel) - grid_shift = np.array( - [ - sicd.ImageData.SCPPixel.Row - sicd.ImageData.FirstRow, - sicd.ImageData.SCPPixel.Col - sicd.ImageData.FirstCol, - ] - ) - grid_mult = np.array([sicd.Grid.Row.SS, sicd.Grid.Col.SS]) - scp_pos = np.array([sicd.GeoData.SCP.ECF.X, sicd.GeoData.SCP.ECF.Y, sicd.GeoData.SCP.ECF.Z]) - scp_hae = sicd.GeoData.SCP.LLH.HAE - rrdot_offset = cls.calculate_range_range_rate_offset(scp_pos, arp_pos, arp_vel, coa_time) - transform_matrix = cls.calculate_transform_matrix(sicd.PFA, coa_time) - beta0_coeff = sicd.Radiometric.BetaZeroSFPoly.Coefs - sigma0_coeff = sicd.Radiometric.SigmaZeroSFPoly.Coefs - center = Point(sicd.GeoData.SCP.LLH.Lon, sicd.GeoData.SCP.LLH.Lat) - azimuth_angle, elevation_angle = cls.calculate_look_angles(arp_pos, scp_pos) - assert np.isclose(azimuth_angle, sicd.SCPCOA.AzimAng), 'Azimuth angle does not match SCPCOA AzimuthAng' - assert np.isclose(elevation_angle, sicd.SCPCOA.GrazeAng), 'Elevation angle does not match SCPCOA ElevationAng' - look_angle = int(azimuth_angle + 180) % 360 - umbra_sicd = cls( - id=Path(file_path).with_suffix('').name, - file_path=file_path, - footprint=footprint, - shape=(sicd.ImageData.NumRows, sicd.ImageData.NumCols), - wavelength=wavelength, - polarization=polarization, - lookside=lookside, - look_angle=look_angle, - incidence_angle=int(np.round(sicd.SCPCOA.IncidenceAng)), - scp_index=(sicd.ImageData.SCPPixel.Row, sicd.ImageData.SCPPixel.Col), - scp_time=scp_time, - scp_pos=scp_pos, - center=center, - scp_hae=scp_hae, - coa_time=coa_time, - arp_pos=arp_pos, - arp_vel=arp_vel, - grid_shift=grid_shift, - grid_mult=grid_mult, - rrdot_offset=rrdot_offset, - transform_matrix=transform_matrix, - transform_matrix_inv=np.linalg.inv(transform_matrix), - orbit=orbit, - beta0_coeff=beta0_coeff, - sigma0_coeff=sigma0_coeff, - ) - return umbra_sicd - - def rowcol2geo(self, rc: np.ndarray, hae: float = None) -> np.ndarray: - """Transforma (row, col) coordinates to ECEF coordinates. - - Args: - rc: Tuple of (row, col) coordinates - hae: Height above ellipsoid (meters) - - Returns: - np.ndarray: ECEF coordinates - """ - if hae is None: - hae = self.scp_hae - - dem = isce3.geometry.DEMInterpolator(hae) - elp = isce3.core.Ellipsoid() - rgaz = (rc - self.grid_shift[None, :]) * self.grid_mult[None, :] - rrdot = np.dot(self.transform_matrix, rgaz.T) + self.rrdot_offset[:, None] - side = isce3.core.LookSide(1) if self.lookside == 'left' else isce3.core.LookSide(-1) - pts_ecf = [] - wvl = 1.0 - for pt in rrdot.T: - r = pt[0] - dop = -pt[1] * 2 / wvl - llh = isce3.geometry.rdr2geo(0.0, r, self.orbit, side, dop, wvl, dem, threshold=1.0e-8, maxiter=50) - pts_ecf.append(elp.lon_lat_to_xyz(llh)) - return np.vstack(pts_ecf) - - def geo2rowcol(self, xyz: np.ndarray) -> np.ndarray: - """Transform ECEF xyz to (row, col). - - Args: - xyz: ECEF coordinates - - Returns: - (row, col) coordinates - """ - rrdot = np.zeros((2, xyz.shape[0])) - rrdot[0, :] = np.linalg.norm(xyz - self.arp_pos[None, :], axis=1) - rrdot[1, :] = np.dot(-self.arp_vel, (xyz - self.arp_pos[None, :]).T) / rrdot[0, :] - rgaz = np.dot(self.transform_matrix_inv, (rrdot - self.rrdot_offset[:, None])) - rgaz /= self.grid_mult[:, None] - rgaz += self.grid_shift[:, None] - row_col = rgaz.T.copy() - return row_col - - def get_geogrid(self, x_spacing, y_spacing): - ecef = pyproj.CRS(4978) # ECEF on WGS84 Ellipsoid - lla = pyproj.CRS(4979) # WGS84 lat/lon/ellipsoid height - local_utm = pyproj.CRS(get_point_epsg(self.center.y, self.center.x)) - lla2utm = pyproj.Transformer.from_crs(lla, local_utm, always_xy=True) - utm2lla = pyproj.Transformer.from_crs(local_utm, lla, always_xy=True) - ecef2lla = pyproj.Transformer.from_crs(ecef, lla, always_xy=True) - - lla_point = (self.center.x, self.center.y) - utm_point = lla2utm.transform(*lla_point) - utm_point_shift = (utm_point[0] + x_spacing, utm_point[1]) - lla_point_shift = utm2lla.transform(*utm_point_shift) - x_spacing = lla_point_shift[0] - lla_point[0] - y_spacing = -1 * x_spacing - - points = np.array([(0, 0), (0, self.shape[1]), self.shape, (self.shape[0], 0)]) - geos = self.rowcol2geo(points) - - points = np.vstack(ecef2lla.transform(geos[:, 0], geos[:, 1], geos[:, 2])).T - minx, maxx = np.min(points[:, 0]), np.max(points[:, 0]) - miny, maxy = np.min(points[:, 1]), np.max(points[:, 1]) - - width = (maxx - minx) // x_spacing - length = (maxy - miny) // np.abs(y_spacing) - geogrid = isce3.product.GeoGridParameters( - start_x=float(minx), - start_y=float(maxy), - spacing_x=float(x_spacing), - spacing_y=float(y_spacing), - length=int(length), - width=int(width), - epsg=4326, - ) - return geogrid - - def get_xrow_ycol(self) -> np.ndarray: - """Calculate xrow and ycol for the umbra_sicd.""" - irow = np.tile(np.arange(self.shape[0]), (self.shape[1], 1)).T - irow -= self.scp_index[0] - xrow = irow * self.grid_mult[0] - - icol = np.tile(np.arange(self.shape[1]), (self.shape[0], 1)) - icol -= self.scp_index[1] - ycol = icol * self.grid_mult[1] - return xrow, ycol - - def load_data(self): - """Load data from the UMBRA SICD file.""" - reader = SICDReader(str(self.file_path)) - data = reader[:, :] - return data - - def load_corrected_data(self, correction='sigma0'): - """Load data with a specified correction""" - if correction == 'sigma0': - coeff = self.sigma0_coeff - elif correction == 'beta0': - coeff = self.beta0_coeff - else: - raise ValueError(f'Unknown correction type: {correction}') - xrow, ycol = self.get_xrow_ycol() - scale_factor = polyval2d(xrow, ycol, coeff) - data = self.load_data() - power = data.real**2 + data.imag**2 - corrected = power * scale_factor - return corrected - - def prep_umbra(granule_path: Path, work_dir: Optional[Path] = None) -> Path: """Prepare data for burst-based processing. From 5d35af1e81f60a21dee77a907707258bafc45c43 Mon Sep 17 00:00:00 2001 From: Forrest Williams Date: Fri, 30 May 2025 15:23:50 -0500 Subject: [PATCH 11/11] update changelog --- CHANGELOG.md | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 89b4db2..7b1b361 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -6,11 +6,17 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/), and this project adheres to [PEP 440](https://www.python.org/dev/peps/pep-0440/) and uses [Semantic Versioning](https://semver.org/spec/v2.0.0.html). -## [0.2.0] +## [0.3.0] + +### Changed +* Changed PFA workflow to use a sublclass of the SicdSlc class ### Added * Cal/Val scripts for performing absolute geolocation error (ALE) assessment +### Fixed +* Property assignment issues that caused bugs in the SicdRzdSlc class + ## [0.2.0] ### Added