# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

#

# Author: Simran Sangha, David Bekaert, Alex Fore

# Copyright (c) 2023, by the California Institute of Technology. ALL RIGHTS

# RESERVED. United States Government Sponsorship acknowledged.

#

# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

"""

Extract and organize specified layer(s).

If no layer is specified, extract product bounding box shapefile(s)

"""

import os

import sys

import glob

import time

import copy

import json

import shutil

import logging

import datetime

import tarfile

import subprocess

import threading

import dask

import rioxarray

import rasterio

import osgeo

import osgeo_utils.gdal_calc

import pyproj

import numpy as np

import scipy.interpolate

import shapely.geometry

import ARIAtools.product

import ARIAtools.util.ionosphere

import ARIAtools.util.interp

import ARIAtools.util.vrt

import ARIAtools.util.shp

import ARIAtools.util.misc

import ARIAtools.util.seq_stitch

from ARIAtools.constants import ARIA_PX_SIZES

LOGGER = logging.getLogger(__name__)

# metadata layer quality check, correction applied if necessary

# only apply to geometry layers and prods derived from older ISCE versions

GEOM_LYRS = ['bPerpendicular', 'bParallel', 'incidenceAngle',

'lookAngle', 'azimuthAngle']

class MetadataQualityCheck:

"""

Metadata quality control function.

Artifacts recognized based off of covariance of cross-profiles.

Bug-fix varies based off of layer of interest.

Verbose mode generates a series of quality control plots with

these profiles.

"""

def __init__(self, data_array, prod_key, outname, verbose=None):

# Pass inputs

self.data_array = data_array

self.prod_key = prod_key

self.outname = outname

self.verbose = verbose

self.data_array_band = data_array.GetRasterBand(1).ReadAsArray()

# mask by nodata value

no_data_value = self.data_array.GetRasterBand(1).GetNoDataValue()

self.data_array_band = np.ma.masked_where(

self.data_array_band == no_data_value, self.data_array_band)

# Run class

self.__run__()

def __truncateArray__(self, data_array_band, Xmask, Ymask):

# Mask columns/rows which are entirely made up of 0s

# first must crop all columns with no valid values

nancols = np.all(data_array_band.mask == True, axis=0)

data_array_band = data_array_band[:, ~nancols]

Xmask = Xmask[:, ~nancols]

Ymask = Ymask[:, ~nancols]

# first must crop all rows with no valid values

nanrows = np.all(data_array_band.mask == True, axis=1)

data_array_band = data_array_band[~nanrows]

Xmask = Xmask[~nanrows]

Ymask = Ymask[~nanrows]

return data_array_band, Xmask, Ymask

def __getCovar__(self, prof_direc, profprefix=''):

from scipy.stats import linregress

# Mask columns/rows which are entirely made up of 0s

if (self.data_array_band.mask.size != 1 and

True in self.data_array_band.mask):

Xmask, Ymask = np.meshgrid(

np.arange(0, self.data_array_band.shape[1], 1),

np.arange(0, self.data_array_band.shape[0], 1))

self.data_array_band, Xmask, Ymask = self.__truncateArray__(

self.data_array_band, Xmask, Ymask)

# append prefix for plot names

prof_direc = profprefix + prof_direc

# Cycle between range and azimuth profiles

rsquaredarr = []

std_errarr = []

# iterate through transpose of matrix if looking in azimuth

data_array_band = (

self.data_array_band.T if 'azimuth' in prof_direc else

self.data_array_band)

for i in enumerate(data_array_band):

mid_line = i[1]

xarr = np.array(range(len(mid_line)))

# remove masked values from slice

if mid_line.mask.size != 1:

if True in mid_line.mask:

xarr = xarr[~mid_line.mask]

mid_line = mid_line[~mid_line.mask]

# chunk array to better isolate artifacts

chunk_size = 4

for j in range(0, len(mid_line.tolist()), chunk_size):

chunk = mid_line.tolist()[j:j + chunk_size]

xarr_chunk = xarr[j:j + chunk_size]

# make sure each iteration contains at least minimum number of

# elements

if j == range(0, len(mid_line.tolist()), chunk_size)[-2] and \

len(mid_line.tolist()) % chunk_size != 0:

chunk = mid_line.tolist()[j:]

xarr_chunk = xarr[j:]

# linear regression and get covariance

slope, bias, rsquared, p_value, std_err = linregress(

xarr_chunk, chunk)

rsquaredarr.append(abs(rsquared)**2)

std_errarr.append(std_err)

# terminate early if last iteration would have small chunk size

if len(chunk) > chunk_size:

break

# exit loop/make plots in verbose mode if R^2 and standard error

# anomalous, or if on last iteration

if (min(rsquaredarr) < 0.9 and max(std_errarr) > 0.01) or \

(i[0] == (len(data_array_band) - 1)):

if self.verbose:

# Make quality-control plots

import matplotlib.pyplot as plt

ax0 = plt.figure().add_subplot(111)

ax0.scatter(xarr, mid_line, c='k', s=7)

refline = np.linspace(min(xarr), max(xarr), 100)

ax0.plot(refline, (refline * slope) + bias,

linestyle='solid', color='red')

ax0.set_ylabel('%s array' % (self.prod_key))

ax0.set_xlabel('distance')

ax0.set_title('Profile along %s' % (prof_direc))

ax0.annotate(

'R\u00b2 = %f\nStd error= %f' % (

min(rsquaredarr), max(std_errarr)),

(0, 1), xytext=(4, -4), xycoords='axes fraction',

textcoords='offset points', fontweight='bold',

ha='left', va='top')

if min(rsquaredarr) < 0.9 and max(std_errarr) > 0.01:

ax0.annotate(

'WARNING: R\u00b2 and standard error\nsuggest '

'artifact exists', (1, 1), xytext=(4, -4),

xycoords='axes fraction',

textcoords='offset points', fontweight='bold',

ha='right', va='top')

plt.margins(0)

plt.tight_layout()

plt.savefig(os.path.join(

os.path.dirname(os.path.dirname(self.outname)),

'metadatalyr_plots', self.prod_key,

os.path.basename(self.outname) + '_%s.png' % (

prof_direc)))

plt.close()

break

return rsquaredarr, std_errarr

def __run__(self):

from scipy.linalg import lstsq

# Get R^2/standard error across range

rsquaredarr_rng, std_errarr_rng = self.__getCovar__('range')

# Get R^2/standard error across azimuth

rsquaredarr_az, std_errarr_az = self.__getCovar__('azimuth')

# filter out normal values from arrays

rsquaredarr = [0.97]

std_errarr = [0.0015]

if min(rsquaredarr_rng) < 0.97 and max(std_errarr_rng) > 0.0015:

rsquaredarr.append(min(rsquaredarr_rng))

std_errarr.append(max(std_errarr_rng))

if min(rsquaredarr_az) < 0.97 and max(std_errarr_az) > 0.0015:

rsquaredarr.append(min(rsquaredarr_az))

std_errarr.append(max(std_errarr_az))

# if R^2 and standard error anomalous, fix array

if min(rsquaredarr) < 0.97 and max(std_errarr) > 0.0015:

# Cycle through each band

for i in range(1, 5):

self.data_array_band = self.data_array.GetRasterBand(

i).ReadAsArray()

# mask by nodata value

no_data_value = self.data_array.GetRasterBand(

i).GetNoDataValue()

self.data_array_band = np.ma.masked_where(

self.data_array_band == no_data_value,

self.data_array_band)

negs_percent = ((self.data_array_band < 0).sum()

/ self.data_array_band.size) * 100

# Unique bug-fix for bPerp layers with sign-flips

if ((self.prod_key == 'bPerpendicular' and

min(rsquaredarr) < 0.8 and max(std_errarr) > 0.1) and

(negs_percent != 100 or negs_percent != 0)):

# Circumvent Bperp sign-flip bug by comparing percentage of

# positive and negative values

self.data_array_band = abs(self.data_array_band)

if negs_percent > 50:

self.data_array_band *= -1

else:

# regular grid covering the domain of the data

X, Y = np.meshgrid(

np.arange(0, self.data_array_band.shape[1], 1),

np.arange(0, self.data_array_band.shape[0], 1))

Xmask, Ymask = np.meshgrid(

np.arange(0, self.data_array_band.shape[1], 1),

np.arange(0, self.data_array_band.shape[0], 1))

# best-fit linear plane: for very large artifacts, must

# mask array for outliers to get best fit

if min(rsquaredarr) < 0.85 and max(std_errarr) > 0.0015:

maj_percent = ((self.data_array_band <

self.data_array_band.mean()).sum()

/ self.data_array_band.size) * 100

# mask all values above mean

if maj_percent > 50:

self.data_array_band = np.ma.masked_where(

self.data_array_band > self.data_array_band.mean(),

self.data_array_band)

# mask all values below mean

else:

self.data_array_band = np.ma.masked_where(

self.data_array_band < self.data_array_band.mean(),

self.data_array_band)

# Mask columns/rows which are entirely made up of 0s

if (self.data_array_band.mask.size != 1 and

True in self.data_array_band.mask):

self.data_array_band, Xmask, Ymask = \

self.__truncateArray__(

self.data_array_band, Xmask, Ymask)

# truncated grid covering the domain of the data

Xmask = Xmask[~self.data_array_band.mask]

Ymask = Ymask[~self.data_array_band.mask]

self.data_array_band = self.data_array_band[

~self.data_array_band.mask]

XX = Xmask.flatten()

YY = Ymask.flatten()

A = np.c_[XX, YY, np.ones(len(XX))]

C, _, _, _ = lstsq(A, self.data_array_band.data.flatten())

# evaluate it on grid

self.data_array_band = C[0] * X + C[1] * Y + C[2]

# mask by nodata value

no_data_value = self.data_array.GetRasterBand(

i).GetNoDataValue()

self.data_array_band = np.ma.masked_where(

self.data_array_band == no_data_value,

self.data_array_band)

np.ma.set_fill_value(self.data_array_band, no_data_value)

# update band

self.data_array.GetRasterBand(i).WriteArray(

self.data_array_band.filled())

# Pass warning and get R^2/standard error across range/azimuth

# (only do for first band)

if i == 1:

# make sure appropriate unit is passed to print statement

lyrunit = "\N{DEGREE SIGN}"

if (self.prod_key == 'bPerpendicular' or

self.prod_key == 'bParallel'):

lyrunit = 'm'

LOGGER.warning((

"%s layer for IFG %s has R\u00b2 of %.4f and standard "

"error of %.4f%s, automated correction applied") % (

self.prod_key, os.path.basename(self.outname),

min(rsquaredarr), max(std_errarr), lyrunit))

rsquaredarr_rng, std_errarr_rng = self.__getCovar__(

'range', profprefix='corrected')

rsquaredarr_az, std_errarr_az = self.__getCovar__(

'azimuth', profprefix='corrected')

self.data_array_band = None

# --- CLOSE THE CLASS-LEVEL POINTER ---

# Capture the dataset to return it, then explicitly sever

# the class's internal link to the GDAL memory object.

safe_return_array = self.data_array

self.data_array = None

return safe_return_array

def crop_only_manager(outname, lyrname, ifg_tag, gdal_warp_kwargs):

"""

Manage cropping of existing, extracted layers

"""

LOGGER.debug('Cropping %s - %s', ifg_tag, lyrname)

# Crop

gdal_warp_kwargs['format'] = 'ENVI'

warp_options = osgeo.gdal.WarpOptions(**gdal_warp_kwargs)

ds = osgeo.gdal.Warp(

outname + '_crop', outname + '.vrt', options=warp_options

)

ds = None

for crop_name in glob.glob(outname + '_crop*'):

fname = os.path.basename(crop_name).replace('_crop', '')

fname = os.path.join(os.path.dirname(crop_name), fname)

os.rename(crop_name, fname)

# Update VRT

ds_trans = osgeo.gdal.Translate(outname + '.vrt', outname, format='VRT')

ds_trans = None

return

def merged_productbbox(

metadata_dict, product_dict, workdir='./', bbox_file=None,

croptounion=False, num_threads='2', minimumOverlap=0.0081,

verbose=None, runlog=None):

"""

Extract/merge productBoundingBox layers for each pair.

Also update dict, report common track bbox

(default is to take common intersection, but user may specify union),

report common track union to accurately interpolate metadata fields,

and expected shape for DEM.

"""

# If specified workdir doesn't exist, create it

os.makedirs(workdir, exist_ok=True)

# determine if NISAR GUNW

is_nisar_file = False

track_fileext = product_dict[0]['unwrappedPhase'][0].split('"')[1]

if track_fileext.endswith('.h5'):

is_nisar_file = True

# If specified, check if user's bounding box meets minimum threshold area

lyr_proj = int(metadata_dict[0]['projection'][0])

if bbox_file is not None:

user_bbox = ARIAtools.util.shp.open_shp(bbox_file)

overlap_area = ARIAtools.util.shp.shp_area(user_bbox, lyr_proj)

if overlap_area < minimumOverlap:

raise Exception(f"User bound box {bbox_file} has an area of only "

f"{overlap_area}km\u00b2, below specified "

f"minimum threshold area "

f"{minimumOverlap}km\u00b2")

# Check if product bounding box exists from previous run

prods_TOTbbox = os.path.join(workdir, 'productBoundingBox.json')

prods_TOTbbox_metadatalyr = os.path.join(

workdir, 'productBoundingBox_croptounion_formetadatalyr.json')

if os.path.exists(prods_TOTbbox) \

and os.path.exists(prods_TOTbbox_metadatalyr):

exist_bbox = ARIAtools.util.shp.open_shp(prods_TOTbbox)

exist_metadatalyr = \

ARIAtools.util.shp.open_shp(prods_TOTbbox_metadatalyr)

# Save copy of file to disk

if runlog:

log_data = runlog.load()

run_time = log_data['run_times'][-1]

else:

run_time = datetime.datetime.now().strftime('%Y%m%d-%H%M%S')

copy_ext = f"{run_time}.json"

bbox_copyname = prods_TOTbbox.replace('.json', copy_ext)

shutil.copyfile(prods_TOTbbox, bbox_copyname)

metadatalyr_copyname = \

prods_TOTbbox_metadatalyr.replace('.json', copy_ext)

shutil.copyfile(prods_TOTbbox_metadatalyr, metadatalyr_copyname)

LOGGER.debug(

'Copying existing productBoundingBox to %s', bbox_copyname)

LOGGER.debug(

'Copying existing metadatalyr to %s', metadatalyr_copyname)

else:

exist_bbox = None

# Extract/merge productBoundingBox layers

for scene in product_dict:

# Get pair name, expected in dictionary

pair_name = scene["pair_name"][0]

outname = os.path.join(workdir, pair_name + '.json')

if os.path.exists(outname):

os.remove(outname)

# Create union of productBoundingBox layers

for prods_bbox in scene["productBoundingBox"]:

if os.path.exists(outname):

union_bbox = ARIAtools.util.shp.open_shp(outname)

prods_bbox = prods_bbox.union(union_bbox)

ARIAtools.util.shp.save_shp(

outname, prods_bbox, lyr_proj)

scene["productBoundingBox"] = [outname]

prods_TOTbbox = os.path.join(workdir, 'productBoundingBox.json')

# Need to track bounds of max extent

# to avoid metadata interpolation issues

sceneareas = [

ARIAtools.util.shp.open_shp(i['productBoundingBox'][0]).area

for i in product_dict]

ind_max_area = sceneareas.index(max(sceneareas))

product_bbox = product_dict[ind_max_area]['productBoundingBox'][0]

ARIAtools.util.shp.save_shp(

prods_TOTbbox_metadatalyr, ARIAtools.util.shp.open_shp(product_bbox),

lyr_proj)

# Initiate intersection file with bbox, if bbox specified

if bbox_file is not None:

ARIAtools.util.shp.save_shp(

prods_TOTbbox, ARIAtools.util.shp.open_shp(bbox_file),

lyr_proj)

# Initiate intersection with largest scene, if bbox NOT specified

else:

ARIAtools.util.shp.save_shp(

prods_TOTbbox, ARIAtools.util.shp.open_shp(product_bbox),

lyr_proj)

rejected_scenes = []

for scene in product_dict:

scene_obj = scene['productBoundingBox'][0]

prods_bbox = ARIAtools.util.shp.open_shp(scene_obj)

total_bbox = ARIAtools.util.shp.open_shp(prods_TOTbbox)

total_bbox_metadatalyr = ARIAtools.util.shp.open_shp(

prods_TOTbbox_metadatalyr)

# Generate footprint for the union of all products

if croptounion:

# Get union

total_bbox = total_bbox.union(prods_bbox)

total_bbox_metadatalyr = total_bbox_metadatalyr.union(prods_bbox)

# Save to file

ARIAtools.util.shp.save_shp(

prods_TOTbbox, total_bbox, lyr_proj)

ARIAtools.util.shp.save_shp(

prods_TOTbbox_metadatalyr, total_bbox_metadatalyr,

lyr_proj)

# Generate footprint for the common intersection of all products

else:

# Now pass track intersection for cutline

prods_bbox = prods_bbox.intersection(total_bbox)

# Estimate percentage of overlap with bbox

if prods_bbox.geom_type == 'MultiPolygon':

LOGGER.debug(

f'Rejected scene {scene_obj} is type MultiPolygon')

rejected_scenes.append(product_dict.index(scene))

os.remove(scene_obj)

continue

if prods_bbox.bounds == () or prods_bbox.is_empty:

LOGGER.debug(f'Rejected scene {scene_obj} '

f'has no common overlap with bbox')

rejected_scenes.append(product_dict.index(scene))

os.remove(scene_obj)

else:

overlap_area = ARIAtools.util.shp.shp_area(

prods_bbox, lyr_proj)

# Kick out scenes below specified overlap threshold

if overlap_area < minimumOverlap:

LOGGER.debug(f'Rejected scene {scene_obj} has only '

f'{overlap_area}km\u00b2 overlap with bbox')

rejected_scenes.append(product_dict.index(scene))

os.remove(scene_obj)

else:

ARIAtools.util.shp.save_shp(

prods_TOTbbox, prods_bbox, lyr_proj)

# Need to track bounds of max extent

# to avoid metadata interpolation issues

total_bbox_metadatalyr = total_bbox_metadatalyr.union(

ARIAtools.util.shp.open_shp(

scene['productBoundingBox'][0]))

ARIAtools.util.shp.save_shp(

prods_TOTbbox_metadatalyr, total_bbox_metadatalyr,

lyr_proj)

# Remove scenes with insufficient overlap w.r.t. bbox

if rejected_scenes != []:

LOGGER.warning(("%d out of %d interferograms rejected for not "

"meeting specified spatial thresholds"),

len(rejected_scenes), len(product_dict))

metadata_dict = [

i for j, i in enumerate(metadata_dict) if j not in rejected_scenes]

product_dict = [

i for j, i in enumerate(product_dict) if j not in rejected_scenes]

if product_dict == []:

raise Exception(

'No common track overlap, footprints cannot be generated.')

# If bbox specified, intersect with common track intersection/union

if bbox_file is not None:

user_bbox = ARIAtools.util.shp.open_shp(bbox_file)

total_bbox = ARIAtools.util.shp.open_shp(prods_TOTbbox)

user_bbox = user_bbox.intersection(total_bbox)

ARIAtools.util.shp.save_shp(

prods_TOTbbox, user_bbox, lyr_proj)

else:

bbox_file = prods_TOTbbox

# Compare current bbox to existing bbox

if exist_bbox and runlog:

exist_area = ARIAtools.util.shp.shp_area(exist_bbox, lyr_proj)

# Calculate overlap area

new_bbox = ARIAtools.util.shp.open_shp(prods_TOTbbox)

new_area = ARIAtools.util.shp.shp_area(new_bbox, lyr_proj)

area_ratio = new_area / exist_area

olap_bbox = exist_bbox.intersection(new_bbox)

olap_area = ARIAtools.util.shp.shp_area(olap_bbox, lyr_proj)

# Compare areas

delta_area = np.abs(olap_area - exist_area)

delta_area = np.round(delta_area * 1E7) * 1E-7

olap_ratio = olap_area / exist_area

olap_ratio = np.round(olap_ratio * 1E7) * 1E-7

LOGGER.debug(

'Area difference (|prev - new|): %f km\u00b2', delta_area)

LOGGER.debug('Area ratio (new/prev): %f', area_ratio)

LOGGER.debug('Overlap ratio (new/prev): %f', olap_ratio)

if (delta_area != 0.0) or (olap_ratio != 1.0):

LOGGER.warning('Product bbox changed in size from previous '

'run %f vs %f', new_area, exist_area)

if shapely.equals(new_bbox, exist_bbox):

# Same bbox within machine precision

update_mode = 'skip'

elif olap_ratio < 1.0:

# For smaller bbox, need to crop

update_mode = 'crop_only'

else:

# If no prior products exist, or new AOI is larger

update_mode = 'full_extract'

runlog.update('update_mode', update_mode)

LOGGER.info('Update mode: %s', update_mode)

# Warp the first scene with the output-bounds defined above

# ensure output-bounds are an integer multiple of interferometric grid

# and adjust if necessary

OG_bounds = list(

ARIAtools.util.shp.open_shp(bbox_file).bounds)

gdal_warp_kwargs = {

'format': 'MEM', 'multithread': False, 'dstSRS': f'EPSG:{lyr_proj}'}

ds_vrt = osgeo.gdal.BuildVRT('', product_dict[0]['unwrappedPhase'][0])

with osgeo.gdal.config_options({"GDAL_NUM_THREADS": num_threads}):

warp_options = osgeo.gdal.WarpOptions(**gdal_warp_kwargs)

ds = osgeo.gdal.Warp('', ds_vrt, options=warp_options)

arrres = [abs(ds.GetGeoTransform()[1]),

abs(ds.GetGeoTransform()[-1])]

ds = None

ds_vrt = None

# Adjust arrres to supported resolution

for i, res in enumerate(arrres):

res_ndx = np.argmin([

np.abs(res - px_size) for px_size in ARIA_PX_SIZES

])

arrres[i] = ARIA_PX_SIZES[res_ndx]

# warp again with fixed transform and bounds

gdal_warp_kwargs['outputBounds'] = OG_bounds

gdal_warp_kwargs['xRes'] = arrres[0]

gdal_warp_kwargs['yRes'] = arrres[1]

gdal_warp_kwargs['targetAlignedPixels'] = True

ds_vrt = osgeo.gdal.BuildVRT('', product_dict[0]['unwrappedPhase'][0])

with osgeo.gdal.config_options({"GDAL_NUM_THREADS": num_threads}):

warp_options = osgeo.gdal.WarpOptions(**gdal_warp_kwargs)

ds = osgeo.gdal.Warp('', ds_vrt, options=warp_options)

# Get shape of full res layers

arrshape = [ds.RasterYSize, ds.RasterXSize]

ds_gt = ds.GetGeoTransform()

new_bounds = [ds_gt[0], ds_gt[3] + (ds_gt[-1] * arrshape[0]),

ds_gt[0] + (ds_gt[1] * arrshape[1]), ds_gt[3]]

if OG_bounds != new_bounds:

# Use shapely to make list

user_bbox = shapely.geometry.Polygon(np.column_stack((

np.array([new_bounds[0], new_bounds[2], new_bounds[2],

new_bounds[0], new_bounds[0]]),

np.array([new_bounds[1], new_bounds[1], new_bounds[3],

new_bounds[3], new_bounds[1]]))))

# Save polygon in shapefile

bbox_file = os.path.join(

os.path.dirname(workdir), 'user_bbox.json')

ARIAtools.util.shp.save_shp(

bbox_file, user_bbox, lyr_proj)

total_bbox = ARIAtools.util.shp.open_shp(prods_TOTbbox)

user_bbox = user_bbox.intersection(total_bbox)

ARIAtools.util.shp.save_shp(

prods_TOTbbox, user_bbox, lyr_proj)

# Get projection of full res layers

proj = ds.GetProjection()

ds = None

ds_vrt = None

# Run additional checks and update runlog if provided

if runlog is None:

update_mode = 'full_extract'

else:

log_data = runlog.load()

update_mode = log_data['update_mode']

# Check other parameters

if ('arrres' in log_data.keys()) \

and (arrres != log_data['arrres']):

runlog.update('update_mode', 'full_extract')

LOGGER.warning('arrres has changed. '

'Setting update mode to full_extract.')

if ('lyr_proj' in log_data.keys()) \

and (lyr_proj != log_data['lyr_proj']):

runlog.update('update_mode', 'full_extract')

LOGGER.warning('lyr_proj has changed. '

'Setting update mode to full_extract.')

# Update log

runlog.update('croptounion', croptounion)

runlog.update('prods_TOTbbox', prods_TOTbbox)

runlog.update('prods_TOTbbox_metadatalyr', prods_TOTbbox_metadatalyr)

runlog.update('arrres', arrres)

runlog.update('lyr_proj', lyr_proj)

runlog.update('metadata_dict', metadata_dict)

runlog.update('product_dict', product_dict)

runlog.update('bbox_file', bbox_file)

runlog.update('is_nisar_file', is_nisar_file)

return (metadata_dict, product_dict, bbox_file, prods_TOTbbox,

prods_TOTbbox_metadatalyr, arrres, proj, update_mode,

is_nisar_file)

def create_raster_from_gunw(fname, data_lis, proj, driver, hgt_field=None,

sign_multiplier=1, dem=None):

"""Wrapper to create raster and apply projection using Rioxarray (Safe)"""

# --- Height-based band subsetting optimisation ---

# If a DEM is provided, subset to only the height bands that span the

# DEM elevation range *before* the expensive remote I/O Warp.

subset_vrts = []

effective_data_lis = data_lis

subsetted_heightsMeta = None

if (dem is not None and hgt_field

and not os.environ.get('ARIA_DISABLE_HEIGHT_SUBSET')):

try:

heightsMeta_str = ARIAtools.util.vrt.get_hgt_meta(

data_lis[0], hgt_field)

if heightsMeta_str:

heightsMeta = np.array(

heightsMeta_str[1:-1].split(','), dtype='float32')

dem_min, dem_max = (

ARIAtools.util.interp._compute_dem_range(dem))

band_indices = (

ARIAtools.util.interp._get_height_subset_indices(

heightsMeta, dem_min, dem_max, pad=0))

if len(band_indices) < len(heightsMeta):

band_list = [int(i + 1) for i in band_indices]

translate_opts = osgeo.gdal.TranslateOptions(

format='VRT', bandList=band_list)

subsetted_data = []

for idx, src in enumerate(data_lis):

sub_vrt = fname + f'_src{idx}_hsubset.vrt'

ds_sub = osgeo.gdal.Translate(

sub_vrt, src, options=translate_opts)

ds_sub = None

subset_vrts.append(sub_vrt)

subsetted_data.append(sub_vrt)

effective_data_lis = subsetted_data

subsetted_heightsMeta = heightsMeta[band_indices]

except Exception:

pass # fall back to reading all bands

# 1) Build a lightweight reference warp (VRT)

ref_vrt = fname + "_ref.vrt"

ds = osgeo.gdal.Warp(

ref_vrt,

effective_data_lis[0],

format='VRT',

dstSRS=proj,

dstNodata=np.nan,

multithread=False

)

ds = None # Close immediately

# 2) Read the derived resolution

ds = osgeo.gdal.Open(ref_vrt, osgeo.gdal.GA_ReadOnly)

gt = ds.GetGeoTransform()

xres, yres = gt[1], abs(gt[5])

ds = None # Close immediately

# 3) Warp + mosaic to temp Tiff

mosaic_tif = fname + "_warp.tif"

ds = osgeo.gdal.Warp(

mosaic_tif,

effective_data_lis,

format='GTiff',

xRes=xres, yRes=yres,

dstSRS=proj,

dstNodata=np.nan,

multithread=False,

creationOptions=["TILED=YES", "COMPRESS=LZW", "BIGTIFF=IF_SAFER"]

)

ds = None # Close immediately

# 3b) Fix stale NETCDF dimension metadata after band subsetting.

# gdal.Warp propagates the original height dimension metadata

# (e.g. 20 values) even though the data now has fewer bands.

# rioxarray uses this metadata to build coordinates, causing a

# dimension mismatch. Update it to match the actual band count.

if subsetted_heightsMeta is not None:

ds_fix = osgeo.gdal.Open(mosaic_tif, osgeo.gdal.GA_Update)

if ds_fix is not None:

meta = ds_fix.GetMetadata()

for key in list(meta.keys()):

if key.startswith('NETCDF_DIM_') and key.endswith('_VALUES'):

dim_name = key[len('NETCDF_DIM_'):-len('_VALUES')]

new_vals = '{' + ','.join(

str(h) for h in subsetted_heightsMeta) + '}'

ds_fix.SetMetadataItem(key, new_vals)

def_key = f'NETCDF_DIM_{dim_name}_DEF'

if def_key in meta:

old_def = meta[def_key]

dtype_str = old_def.strip('{}[]').split(',')[-1]

ds_fix.SetMetadataItem(

def_key,

'{' + str(len(subsetted_heightsMeta)) +

',' + dtype_str + '}')

ds_fix.FlushCache()

ds_fix = None

# 4) Open with rioxarray (Context Manager prevents locking)

with rioxarray.open_rasterio(mosaic_tif, masked=True) as da:

da = da.rio.write_nodata(np.nan, encoded=True)

# --- FIX: Remove _FillValue from attrs ---

if "_FillValue" in da.attrs:

del da.attrs["_FillValue"]

# -----------------------------------------

# Flip the sign for NISAR convention

if sign_multiplier == -1:

da = da * -1

# Let GDAL use its safe, default thread pool

da.rio.to_raster(fname, driver=driver, crs=proj)

# 5) Clean up (Now safe because 'da' is closed)

if os.path.exists(mosaic_tif):

os.remove(mosaic_tif)

if os.path.exists(ref_vrt):

os.remove(ref_vrt)

# 6) Create VRT file

buildvrt_options = osgeo.gdal.BuildVRTOptions(outputSRS=proj)

ds_vrt = osgeo.gdal.BuildVRT(

fname + '.vrt', fname, options=buildvrt_options

)

ds_vrt = None # Close immediately

# 7) Add height info

if hgt_field is not None:

if subsetted_heightsMeta is not None:

# Write the subsetted height values

hgt_meta = '{' + ','.join(

str(h) for h in subsetted_heightsMeta) + '}'

else:

hgt_meta = ARIAtools.util.vrt.get_hgt_meta(

data_lis[0], hgt_field)

ds_meta_update = osgeo.gdal.Open(

fname + '.vrt', osgeo.gdal.GA_Update

)

ds_meta_update.SetMetadataItem(hgt_field, hgt_meta)

ds_meta_update = None # Close immediately

# Clean up temporary subset VRTs

for v in subset_vrts:

if os.path.exists(v):

os.remove(v)

return

def prep_metadatalayers(

outname, metadata_arr, dem, layer, layers, is_nisar_file=False,

proj='4326', driver='ENVI', model_name=None, sign_multiplier=1):

"""Wrapper to prep metadata layer for extraction"""

if dem is None:

raise Exception('No DEM input specified. '

'Cannot extract 3D imaging geometry '

'layers without DEM to intersect with.')

ifg = os.path.basename(outname)

out_dir = os.path.dirname(outname)

ref_outname = copy.deepcopy(outname)

# ionosphere layer, heights do not exist to exit

if metadata_arr[0].split('/')[-1] == 'ionosphere':

ds_vrt = osgeo.gdal.BuildVRT(outname + '.vrt', metadata_arr)

ds_vrt= None

return [0], None, outname

# capture model if tropo product

if 'tropo' in layer:

if not is_nisar_file:

out_dir = os.path.join(out_dir, model_name)

outname = os.path.join(out_dir, ifg)

if not os.path.exists(out_dir):

os.mkdir(out_dir)

# Get height values

ds_meta = osgeo.gdal.Open(metadata_arr[0])

zdim = ds_meta.GetMetadataItem('NETCDF_DIM_EXTRA')[1:-1]

ds_meta = None # Close

hgt_field = f'NETCDF_DIM_{zdim}_VALUES'

# Check if height layers are consistent

if (

not os.path.exists(outname + ".vrt")

and len(

{

ARIAtools.util.vrt.get_hgt_meta(i, hgt_field)

for i in metadata_arr

}

)

!= 1

):

heights = [

ARIAtools.util.vrt.get_hgt_meta(i, hgt_field)

for i in metadata_arr

]

raise Exception(

"Inconsistent heights for metadata layer(s) "

f"{metadata_arr}; corresponding heights: {heights}"

)

if 'tropo' in layer or layer == 'solidEarthTide':

# get ref and sec paths

if not is_nisar_file:

date_dir = os.path.join(out_dir, 'dates')

if not os.path.exists(date_dir):

os.mkdir(date_dir)

ref_outname = os.path.join(date_dir, ifg.split('_')[0])

sec_outname = os.path.join(date_dir, ifg.split('_')[1])

ref_str = 'reference/' + layer

sec_str = 'secondary/' + layer

sec_metadata_arr = [

i[:-len(ref_str)] + sec_str for i in metadata_arr]

tup_outputs = [

(ref_outname, metadata_arr), (sec_outname, sec_metadata_arr)]

else:

ref_outname = os.path.join(out_dir, ifg)

sec_outname = ref_outname

tup_outputs = [(ref_outname, metadata_arr)]

# write ref and sec files

for i in tup_outputs:

# delete temporary files to circumvent potential inconsistent dims

for j in glob.glob(i[0] + '*'):

if os.path.isfile(j):

os.remove(j)

create_raster_from_gunw(i[0], i[1], proj, driver, hgt_field,

sign_multiplier, dem=dem)

if not is_nisar_file:

# compute differential

generate_diff(

ref_outname, sec_outname, outname, layer, layer, False,

hgt_field, proj, driver, dem=dem)

# write raster to file if it does not exist

if layer in layers:

for i in [ref_outname, sec_outname]:

if not os.path.exists(i):

create_raster_from_gunw(i, [i], proj, driver, hgt_field)

else:

if not os.path.exists(outname + '.vrt'):

if is_nisar_file:

# Need to compute azimuthAngle from

# losUnitVectorX and losUnitVectorY

if layer == 'azimuthAngle':

losx_arr = copy.deepcopy(metadata_arr)

losy_arr = [

path.replace('losUnitVectorX', 'losUnitVectorY')

for path in metadata_arr

]

# Initiate temp los dimension files

losx_name = os.path.join(out_dir, f'temp_{ifg}_losx_arr')

losy_name = os.path.join(out_dir, f'temp_{ifg}_losy_arr')

create_raster_from_gunw(

losx_name, losx_arr, proj, driver, hgt_field,

dem=dem

)

create_raster_from_gunw(

losy_name, losy_arr, proj, driver, hgt_field,

dem=dem

)

# Get NoData value from input

ds_temp = osgeo.gdal.Open(losx_name + '.vrt')

src_nodata = ds_temp.GetRasterBand(1).GetNoDataValue()

ds_temp = None

# Construct the Calc String safely

# If src_nodata exists and is NOT nan, we must mask it

# manually.

if src_nodata is not None and not np.isnan(src_nodata):

# Logic: If (A == nodata) OR (B == nodata), return

# NaN. Else calculate the angle.

calc_cmd = (

f"numpy.where("

f"(A=={src_nodata})|(B=={src_nodata}), "

f"numpy.nan, "

f"numpy.degrees(numpy.arctan2(-B, -A)))"

)

else:

# If input is already NaN (or None), the math

# handles it naturally.

calc_cmd = "numpy.degrees(numpy.arctan2(-B, -A))"

# Compute azimuthAngle from the outputs above

# We map path_x to 'A' and path_y to 'B'

osgeo_utils.gdal_calc.Calc(

A=losx_name + '.vrt',

B=losy_name + '.vrt',

outfile=outname,

calc=calc_cmd,

format=driver,

allBands="A", # processes every band

quiet=True

)

# Manually enforce NoData = NaN on the output header

ds_update = osgeo.gdal.Open(

outname, osgeo.gdal.GA_Update

)

if ds_update: