Add iono

tools/ARIAtools/ARIAProduct.py

@@ -629,13 +629,20 @@ def __NISARmappingVersion__(self, fname, version):
 
         # assign latitude to assist with sorting
         # get product bounding box
+        # and other variables
         lyr_pref = '/science/LSAR/GUNW/grids/frequencyA' + \
                         f'/unwrappedInterferogram/{file_pol}/'
-        proj_location = lyr_pref + 'projection'
+        center_freq = '/science/LSAR/GUNW/grids/frequencyA/centerFrequency'
         with h5py.File(fname[8:], 'r') as hdf_gunw:
+            # get bbox
             latlon_file_bbox = hdf_gunw[sdskeys[0]]
             latlon_file_bbox = latlon_file_bbox[()]
             latlon_file_bbox = loads(latlon_file_bbox)
+            # get center frequency
+            center_freq_var = float(hdf_gunw[center_freq][()])
+        rdrmetadata_dict['centerFrequency'] = center_freq_var
+        rdrmetadata_dict['wavelength'] = 299792458 /  \
+            rdrmetadata_dict['centerFrequency']
         rdrmetadata_dict['centerLatitude'] =  int(latlon_file_bbox.centroid.y)
         rdrmetadata_dict['projection'] = self.projection
         pyproj_transformer = Transformer.from_crs('EPSG:4326',
@@ -646,8 +653,6 @@ def __NISARmappingVersion__(self, fname, version):
         # hardcoded keys
         rdrmetadata_dict['missionID'] = 'NISAR'
         rdrmetadata_dict['productType'] = 'UNW GEO IFG'
-        rdrmetadata_dict['wavelength'] = 0.24
-        rdrmetadata_dict['centerFrequency'] = 1.2699997500604727e+09
         rdrmetadata_dict['slantRangeSpacing'] = 4.461197291666666
         rdrmetadata_dict['slantRangeStart'] = 727415.98682514
         rdrmetadata_dict['slantRangeEnd'] = 791384.81843777
@@ -741,6 +746,9 @@ def __NISARmappingData__(self, fname, rdrmetadata_dict, sdskeys, version):
         # indivdual frames preserved here
         datalyr_dict['productBoundingBoxFrames'] = \
             datalyr_dict['productBoundingBox']
+
+        # Rewrite iono key
+        datalyr_dict['ionosphere'] = datalyr_dict.pop('ionospherePhaseScreen')
 
         return [rdrmetadata_dict, datalyr_dict]
 

tools/ARIAtools/extractProduct.py

@@ -501,7 +501,6 @@ def merged_productbbox(metadata_dict, product_dict, workdir='./',
 
     # If specified, check if user's bounding box meets minimum threshold area
     lyr_proj = int(metadata_dict[0]['projection'][0])
-    print('lyr_proj', lyr_proj)
     if bbox_file is not None:
         user_bbox = open_shapefile(bbox_file, 0, 0)
         overlap_area = shapefile_area(user_bbox, lyr_proj)
@@ -1199,20 +1198,19 @@ def export_products(full_product_dict, proj, bbox_file, prods_TOTbbox, layers,
 
     if list(set.intersection(*map(set,
                                   [layers, ['ionosphere']]))) != []:
-        lyr_prefix = '/science/grids/corrections/derived/ionosphere/ionosphere'
+        # set input keys
         key = 'ionosphere'
         product_dict = \
             [[j[key] for j in full_product_dict if key in j.keys()],
              [j["pair_name"] for j in full_product_dict if key in j.keys()]]
-
         workdir = os.path.join(outDir, key)
         prev_outname = deepcopy(workdir)
         prog_bar = progBar.progressBar(maxValue=len(product_dict[0]),
                                        prefix='Generating: '+key+' - ')
 
-
         lyr_input_dict = dict(input_iono_files = None,
                               arrres = arrres,
+                              epsg = epsg_code,
                               output_iono = None,
                               output_format =  outputFormat, 
                               bounds = bounds,

tools/ARIAtools/ionosphere.py

@@ -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
 
 
@@ -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):
@@ -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'})
@@ -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,
@@ -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,
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
 
-
+