Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add NISAR product support #390

Closed
wants to merge 5 commits into from
Closed

Add NISAR product support #390

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
22e3017
Add NISAR product support
sssangha Mar 20, 2024
1653a29
Add iono
sssangha Mar 20, 2024
a746aa3
Delete unw stitcher
sssangha Mar 20, 2024
78fe05a
Clean tropo key logic
sssangha Mar 21, 2024
8674066
Path bug for reading h5 file
sssangha Mar 22, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
331 changes: 272 additions & 59 deletions tools/ARIAtools/ARIAProduct.py
View file
Open in desktop

Large diffs are not rendered by default.

417 changes: 247 additions & 170 deletions tools/ARIAtools/extractProduct.py
View file
Open in desktop

Large diffs are not rendered by default.

77 changes: 60 additions & 17 deletions tools/ARIAtools/ionosphere.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -7,12 +7,13 @@
#
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

import os
import numpy as np
import xarray as xr
from numpy.typing import NDArray

from typing import List, Optional, Tuple
from osgeo import gdal
from osgeo import gdal, osr
from pathlib import Path


Expand All @@ -23,10 +24,18 @@
_nan_filled_array)


GUNW_LAYERS = {'unwrappedPhase': 'NETCDF:"%s":/science/grids/data/unwrappedPhase',
'coherence': 'NETCDF:"%s":/science/grids/data/coherence',
'connectedComponents': 'NETCDF:"%s":/science/grids/data/connectedComponents',
'ionosphere': 'NETCDF:"%s":/science/grids/corrections/derived/ionosphere/ionosphere'}
GUNW_LAYERS = {
'unwrappedPhase': 'NETCDF:"%s":/science/grids/data/unwrappedPhase',
'coherence': 'NETCDF:"%s":/science/grids/data/coherence',
'connectedComponents': 'NETCDF:"%s":/science/grids/data/connectedComponents',
'ionosphere': 'NETCDF:"%s":/science/grids/corrections/derived/ionosphere/ionosphere'
}

NISAR_GUNW_LAYERS = {
'unwrappedPhase': 'NETCDF:"%s":/science/LSAR/GUNW/grids/frequencyA/unwrappedInterferogram/%s/unwrappedPhase',
'coherence': 'NETCDF:"%s":/science/LSAR/GUNW/grids/frequencyA/unwrappedInterferogram/%s/coherenceMagnitude',
'connectedComponents': 'NETCDF:"%s":/science/LSAR/GUNW/grids/frequencyA/unwrappedInterferogram/%s/connectedComponents',
'ionosphere': 'NETCDF:"%s":/science/LSAR/GUNW/grids/frequencyA/unwrappedInterferogram/%s/ionospherePhaseScreen'}


def fit_surface(data, order=2):
Expand Down Expand Up @@ -75,33 +84,57 @@ def _get_median_offsets2frames(xr_data_list, xr_mask_list, ix1, ix2):
S, N, W, E = _get_overlay(xr_data_list[ix1], xr_data_list[ix2])

# Get overlap
cropped_ds1 = xr_data_list[ix1].ionosphere.sel(y=slice(N, S), x=slice(W, E)).copy()
get_key = list(xr_data_list[ix1].data_vars.keys())[0]
cropped_ds1 = xr_data_list[ix1][get_key].sel(y=slice(N, S), x=slice(W, E)).copy()
cropped_mask1 = xr_mask_list[ix1].mask.sel(y=slice(N, S), x=slice(W, E)).copy()
ds1 = np.ma.masked_array(cropped_ds1.values, mask=~cropped_mask1.values)

cropped_ds2 = xr_data_list[ix2].ionosphere.sel(y=slice(N, S), x=slice(W, E)).copy()
get_key = list(xr_data_list[ix2].data_vars.keys())[0]
cropped_ds2 = xr_data_list[ix2][get_key].sel(y=slice(N, S), x=slice(W, E)).copy()
cropped_mask2 = xr_mask_list[ix2].mask.sel(y=slice(N, S), x=slice(W, E)).copy()
ds2 = np.ma.masked_array(cropped_ds2.values, mask=~cropped_mask2.values)

return np.nanmedian((ds1 - ds2).filled(fill_value=np.nan))

def stitch_ionosphere_frames(input_iono_files: List[str],
direction_N_S: Optional[bool] = True):
proj: Optional[str] = 'EPSG:4326',
direction_N_S: Optional[bool] = True):

# Initalize variables for raster attributes
iono_attr_list = [] # ionosphere raster metadata
iono_xr_list = []
mask_xr_list =[]

# track if product stack is NISAR GUNW or not
nisar_file = False
track_fileext = input_iono_files[0].split(':')[1]
if len(track_fileext.split('.h5')) > 1:
nisar_file = True
# Get polarization
pol_dict = {}
pol_dict['SV'] = 'VV'
pol_dict['SH'] = 'HH'
pol_dict['HHNA'] = 'HH'
basename = os.path.basename(track_fileext)
file_pol = pol_dict[basename.split('_')[10]]

# Loop through files
for iono_file in input_iono_files:
filename = iono_file.split(':')[1]
iono_attr_list.append(get_GUNW_attr(iono_file))
iono_attr_list.append(get_GUNW_attr(iono_file, proj=proj))
iono_xr = xr.open_dataset(iono_file, engine='rasterio').squeeze()

# Generate mask using unwrapPhase connectedComponents
mask_xr = xr.open_dataset(GUNW_LAYERS['connectedComponents'] % filename,
engine='rasterio').squeeze()
if nisar_file:
mask_xr = xr.open_dataset( \
NISAR_GUNW_LAYERS['connectedComponents'] \
%(filename, file_pol),
engine='rasterio').squeeze()
else:
mask_xr = xr.open_dataset( \
GUNW_LAYERS['connectedComponents'] \
% filename,
engine='rasterio').squeeze()
mask = np.bool_(mask_xr.connectedComponents.data != 0)
mask_xr['connectedComponents'].values = mask
mask_xr = mask_xr.rename_vars({'connectedComponents':'mask'})
Expand All @@ -128,10 +161,11 @@ def stitch_ionosphere_frames(input_iono_files: List[str],
# by using only reliable areas (connctedComponents != 0)
for ix1, ix2 in zip(sorted_ix[:-1], sorted_ix[1:]):
diff = _get_median_offsets2frames(iono_xr_list, mask_xr_list, ix1, ix2)
iono_xr_list[ix2]['ionosphere'] += diff
get_key = list(iono_xr_list[ix2].data_vars.keys())[0]
iono_xr_list[ix2][get_key] += diff

# Step 2: Merged ionosphere and mask datasets
data_list = [d.ionosphere.data for d in iono_xr_list]
data_list = [d[list(d.data_vars.keys())[0]].data for d in iono_xr_list]
mask_list = [d.mask.data for d in mask_xr_list]

combined_iono = combine_data_to_single(data_list,
Expand Down Expand Up @@ -164,6 +198,7 @@ def stitch_ionosphere_frames(input_iono_files: List[str],

def export_ionosphere(input_iono_files: List[str],
arrres: List[float],
epsg: Optional[str] = '4326',
output_iono: Optional[str] = './ionosphere',
output_format: Optional[str] = 'ISCE',
bounds: Optional[tuple] = None,
Expand All @@ -180,6 +215,12 @@ def export_ionosphere(input_iono_files: List[str],

# create temp files
temp_iono_out = output_iono.parent / ('temp_' + output_iono.name)

# obtain reference epsg code to assign to intermediate outputs
ref_proj_str = get_GUNW_attr(input_iono_files[0])['PROJECTION']
srs = osr.SpatialReference(wkt=ref_proj_str)
srs.AutoIdentifyEPSG()
ref_proj = srs.GetAuthorityCode(None)

# Create VRT and exit early if only one frame passed,
# and therefore no stitching needed
Expand All @@ -190,7 +231,7 @@ def export_ionosphere(input_iono_files: List[str],
else:
(combined_iono,
snwe, latlon_spacing) = stitch_ionosphere_frames(input_iono_files,
direction_N_S=True)
proj=f'EPSG:{ref_proj}', direction_N_S=True)

# write stitched ionosphere
# outputformat
Expand All @@ -200,7 +241,7 @@ def export_ionosphere(input_iono_files: List[str],

write_GUNW_array(
temp_iono_out, combined_iono, snwe,
format=output_format, verbose=verbose,
format=output_format, epsg=int(ref_proj), verbose=verbose,
update_mode=overwrite, add_vrt=True, nodata=0.0)

# Crop
Expand All @@ -218,6 +259,7 @@ def export_ionosphere(input_iono_files: List[str],
xRes=arrres[0], yRes=arrres[1],
targetAlignedPixels=True,
# cropToCutline = True,
dstSRS=f'EPSG:{epsg}',
outputBounds=bounds
)
ds = None
Expand All @@ -241,11 +283,12 @@ def export_ionosphere(input_iono_files: List[str],
mask = mask_file

mask_array = mask.ReadAsArray()
array = get_GUNW_array(str(output_iono.with_suffix('.vrt')))
array = get_GUNW_array(str(output_iono.with_suffix('.vrt')),
proj=f'EPSG:{epsg}')
update_array = mask_array * array

update_file = gdal.Open(str(output_iono), gdal.GA_Update)
update_file = update_file.GetRasterBand(1).WriteArray(update_array)
update_file = None



Loading
Loading