Skip to content

initial prep NISAR changes #4

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

Open
ehavazli wants to merge 8 commits into
base: main
Choose a base branch
from
Open

initial prep NISAR changes #4

Changes from all commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
665cc2c
change variable name to follow the latest format
ehavazli Oct 2, 2025
5a3d353
use gdal.WarpOptions instead of kwargs to avoid GDAL Warp TypeError
ehavazli Oct 2, 2025
f13ccff
load ionosphere to MintPy
ehavazli Oct 6, 2025
dde21b2
fix pair datetime extraction
ehavazli Oct 8, 2025
df65357
Merge branch 'insarlab:main' into prep_nisar
ehavazli Oct 16, 2025
99f73e3
extract troposphere layer
ehavazli Jan 29, 2026
756a563
Merge branch 'main' into prep_nisar
ehavazli Jan 29, 2026
3d06fca
change format to reference_secondary
ehavazli Jan 29, 2026
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
269 changes: 248 additions & 21 deletions src/mintpy/prep_nisar.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -28,7 +28,8 @@
'unw' : f"{DATASET_ROOT_UNW}/POL/unwrappedPhase",
'cor' : f"{DATASET_ROOT_UNW}/POL/coherenceMagnitude",
'connComp' : f"{DATASET_ROOT_UNW}/POL/connectedComponents",
#'mask' : f"{DATASET_ROOT_UNW}/mask",
'ion' : f"{DATASET_ROOT_UNW}/POL/ionospherePhaseScreen",
'mask' : f"{DATASET_ROOT_UNW}/mask",
'epsg' : f"{DATASET_ROOT_UNW}/projection",
'xSpacing' : f"{DATASET_ROOT_UNW}/xCoordinateSpacing",
'ySpacing' : f"{DATASET_ROOT_UNW}/yCoordinateSpacing",
Expand All @@ -44,9 +45,11 @@
'end_time' : f"{IDENTIFICATION}/referenceZeroDopplerEndTime",
'rdr_xcoord' : f"{RADARGRID_ROOT}/xCoordinates",
'rdr_ycoord' : f"{RADARGRID_ROOT}/yCoordinates",
'rdr_slant_range': f"{RADARGRID_ROOT}/slantRange",
'rdr_slant_range': f"{RADARGRID_ROOT}/referenceSlantRange",
'rdr_height' : f"{RADARGRID_ROOT}/heightAboveEllipsoid",
'rdr_incidence' : f"{RADARGRID_ROOT}/incidenceAngle",
'rdr_wet_tropo' : f"{RADARGRID_ROOT}/wetTroposphericPhaseScreen",
'rdr_hs_tropo' : f"{RADARGRID_ROOT}/hydrostaticTroposphericPhaseScreen",
'bperp' : f"{RADARGRID_ROOT}/perpendicularBaseline",
}

Expand All @@ -71,6 +74,8 @@ def load_nisar(inps):
# output filename
stack_file = os.path.join(inps.out_dir, "inputs/ifgramStack.h5")
geometry_file = os.path.join(inps.out_dir, "inputs/geometryGeo.h5")
ion_stack_file = os.path.join(inps.out_dir, "inputs/ionStack.h5")
tropo_stack_file = os.path.join(inps.out_dir, "inputs/tropoStack.h5")

# get date info: date12_list
date12_list = _get_date_pairs(input_files)
Expand All @@ -88,10 +93,30 @@ def load_nisar(inps):
outfile=stack_file,
inp_files=input_files,
metadata=metadata,
demFile=inps.dem_file,
bbox=bounds,
date12_list=date12_list,
)

# Load ionosphere layer
prepare_stack(
outfile=ion_stack_file,
inp_files=input_files,
metadata=metadata,
demFile=inps.dem_file,
bbox=bounds,
date12_list=date12_list,
)

# Load troposphere layer
prepare_stack(
outfile=tropo_stack_file,
inp_files=input_files,
metadata=metadata,
demFile=inps.dem_file,
bbox=bounds,
date12_list=date12_list,
)
print("Done.")

return
Expand Down Expand Up @@ -239,6 +264,7 @@ def read_subset(inp_file, bbox, geometry=False):
dataset['cor_data'] = ds[DATASETS['cor']][row1:row2, col1:col2]
dataset['conn_comp'] = (ds[DATASETS['connComp']][row1:row2, col1:col2]).astype(np.float32)
dataset['conn_comp'][dataset['conn_comp'] > 254] = np.nan
dataset['ion_data'] = ds[DATASETS['ion']][row1:row2, col1:col2]
dataset['pbase'] = np.nanmean(ds[PROCESSINFO['bperp']][()])
return dataset

Expand Down Expand Up @@ -275,10 +301,20 @@ def read_and_interpolate_geometry(gunw_file, dem_file, xybbox, mask_file=None):
# Warp DEM to the interferograms grid
input_projection = f"EPSG:{dem_src_epsg}"
output_dem = os.path.join(os.path.dirname(dem_file), 'dem_transformed.tif' )
gdal.Warp(output_dem, dem_file, outputBounds=bounds, format='GTiff',
srcSRS=input_projection, dstSRS=output_projection, resampleAlg=gdal.GRA_Bilinear,
width=subset_cols, height=subset_rows,
options=['COMPRESS=DEFLATE'])

warp_opts = gdal.WarpOptions(
format="GTiff",
outputBounds=bounds,
srcSRS=input_projection,
dstSRS=output_projection,
resampleAlg="bilinear",
width=subset_cols, height=subset_rows,
creationOptions=["COMPRESS=DEFLATE"]
)

gdal.Warp(destNameOrDestDS=output_dem,
srcDSOrSrcDSTab=dem_file,
options=warp_opts)

dem_subset_array = gdal.Open(output_dem, gdal.GA_ReadOnly).ReadAsArray()

Expand Down Expand Up @@ -337,11 +373,122 @@ def interpolate_geometry(X_2d, Y_2d, dem, rdr_coords):
return interpolated_slant_range.reshape(length, width), interpolated_incidence_angle.reshape(length, width)


def read_and_interpolate_troposphere(gunw_file, dem_file, xybbox, mask_file=None):
"""Read the DEM and mask, change projection and interpolate to data grid,
interpolate slant range and incidence angle to data grid"""
dem_dataset = gdal.Open(dem_file, gdal.GA_ReadOnly)
proj = gdal.osr.SpatialReference(wkt=dem_dataset.GetProjection())
dem_src_epsg = int(proj.GetAttrValue('AUTHORITY', 1))
del dem_dataset

rdr_coords = {}

with h5py.File(gunw_file, 'r') as ds:
dst_epsg = ds[DATASETS['epsg']][()]
xcoord = ds[DATASETS['xcoord']][xybbox[0]:xybbox[2]]
ycoord = ds[DATASETS['ycoord']][xybbox[1]:xybbox[3]]

rdr_coords['xcoord_radar_grid'] = ds[PROCESSINFO['rdr_xcoord']][()]
rdr_coords['ycoord_radar_grid'] = ds[PROCESSINFO['rdr_ycoord']][()]
rdr_coords['height_radar_grid'] = ds[PROCESSINFO['rdr_height']][()]
# rdr_coords['perp_baseline'] = ds[PROCESSINFO['bperp']][()]
rdr_coords['wet_tropo'] = ds[PROCESSINFO['rdr_wet_tropo']][()]
rdr_coords['hydrostatic_tropo'] = ds[PROCESSINFO['rdr_hs_tropo']][()]

subset_rows = len(ycoord)
subset_cols = len(xcoord)

Y_2d, X_2d = np.meshgrid(ycoord, xcoord, indexing='ij')
bounds = (min(xcoord), min(ycoord), max(xcoord), max(ycoord))
output_projection = f"EPSG:{dst_epsg}"

# Warp DEM to the interferograms grid
input_projection = f"EPSG:{dem_src_epsg}"
output_dem = os.path.join(os.path.dirname(dem_file), 'dem_transformed.tif' )

warp_opts = gdal.WarpOptions(
format="GTiff",
outputBounds=bounds,
srcSRS=input_projection,
dstSRS=output_projection,
resampleAlg="bilinear",
width=subset_cols, height=subset_rows,
creationOptions=["COMPRESS=DEFLATE"]
)

gdal.Warp(destNameOrDestDS=output_dem,
srcDSOrSrcDSTab=dem_file,
options=warp_opts)

dem_subset_array = gdal.Open(output_dem, gdal.GA_ReadOnly).ReadAsArray()

# Calculate Total troposphere and then interpolate hydrostatic_tropo, wet_tropo and total_tropo
total_tropo = interpolate_troposphere(X_2d, Y_2d, dem_subset_array, rdr_coords)

# Read and warp mask if necessary
if mask_file in ['auto', 'None', None]:
print('*** No mask was found ***')
mask_subset_array = np.ones(dem_subset_array.shape, dtype='byte')
else:
try:
mask_dataset = gdal.Open(mask_file, gdal.GA_ReadOnly)
proj = gdal.osr.SpatialReference(wkt=mask_dataset.GetProjection())
mask_src_epsg = int(proj.GetAttrValue('AUTHORITY', 1))
del mask_dataset

# Warp mask to the interferograms grid
input_projection = f"EPSG:{mask_src_epsg}"
output_mask = os.path.join(os.path.dirname(mask_file), 'mask_transformed.tif' )
gdal.Warp(output_mask, mask_file, outputBounds=bounds, format='GTiff',
srcSRS=input_projection, dstSRS=output_projection, resampleAlg=gdal.GRA_Byte,
width=subset_cols, height=subset_rows,
options=['COMPRESS=DEFLATE'])

mask_subset_array = gdal.Open(output_mask, gdal.GA_ReadOnly).ReadAsArray()

except:
raise OSError('*** Mask is not gdal readable ***')


return total_tropo


def interpolate_troposphere(X_2d, Y_2d, dem, rdr_coords):
"""Interpolate slant range and incidence angle"""
pnts = np.stack((dem.flatten(), Y_2d.flatten(), X_2d.flatten()), axis=-1)
length, width = Y_2d.shape

rdr_coords['total_tropo'] = rdr_coords['hydrostatic_tropo'] + rdr_coords['wet_tropo']

# build the interpolator
# interpolator = RegularGridInterpolator((rdr_coords['height_radar_grid'],
# rdr_coords['ycoord_radar_grid'],
# rdr_coords['xcoord_radar_grid']),
# rdr_coords['hydrostatic_tropo'],
# method='linear')

# interpolated_hydrostatic_tropo = interpolator(pnts)

# interpolator = RegularGridInterpolator((rdr_coords['height_radar_grid'],
# rdr_coords['ycoord_radar_grid'],
# rdr_coords['xcoord_radar_grid']),
# rdr_coords['wet_tropo'],
# method='linear')
# interpolated_wet_tropo = interpolator(pnts)

interpolator = RegularGridInterpolator((rdr_coords['height_radar_grid'],
rdr_coords['ycoord_radar_grid'],
rdr_coords['xcoord_radar_grid']),
rdr_coords['total_tropo'],
method='linear')
interpolated_total_tropo = interpolator(pnts)

return interpolated_total_tropo.reshape(length, width)
######################################

def _get_date_pairs(filenames):
str_list = [Path(f).stem for f in filenames]
return [str(f.split('_')[13]) + '_' + str(f.split('_')[11]) for f in str_list]
return [str(f.split('_')[11].split('T')[0]) + '_' + str(f.split('_')[13].split('T')[0]) for f in str_list]


def prepare_geometry(
Expand Down Expand Up @@ -383,11 +530,50 @@ def prepare_geometry(

return meta

# def prepare_troposphere(
# outfile,
# metaFile,
# metadata,
# bbox,
# demFile,
# ):
# """Prepare the geometry file."""
# print("-" * 50)
# print(f"preparing geometry file: {outfile}")

# # copy metadata to meta
# meta = {key: value for key, value in metadata.items()}

# # Read waterMask, LayoverShadowMask, xybbox:
# geo_ds = read_subset(metaFile, bbox, geometry=True)
# hydrostatic_tropo, wet_tropo, total_tropo = read_and_interpolate_troposphere(metaFile, demFile,
# geo_ds['xybbox'])

# length, width = dem_subset_array.shape

# ds_name_dict = {
# "height": [np.float32, (length, width), dem_subset_array],
# "hydrostaticTroposphere": [np.float32, (length, width), hydrostatic_tropo],
# "wetTroposphere": [np.float32, (length, width), wet_tropo],
# "totalTroposphere": [np.float32, (length, width), total_tropo],
# #"shadowMask": [np.bool_, (length, width), geo_ds['mask']],
# #"waterMask": [np.bool_, (length, width), geo_ds['water_mask']],
# }
# # if maskFile:
# # ds_name_dict['waterMask'] = [np.bool_, (length, width), mask]

# # initiate HDF5 file
# meta["FILE_TYPE"] = "ifgramStack"
# meta['STARTING_RANGE'] = np.nanmin(slant_range)
# writefile.layout_hdf5(outfile, ds_name_dict, metadata=meta)

# return meta

def prepare_stack(
outfile,
inp_files,
metadata,
demFile,
bbox,
date12_list
):
Expand Down Expand Up @@ -426,26 +612,67 @@ def prepare_stack(

# writing data to HDF5 file
print(f"writing data to HDF5 file {outfile} with a mode ...")
if "inputs/ifgramStack.h5" in outfile:
with h5py.File(outfile, "a") as f:
prog_bar = ptime.progressBar(maxValue=num_pair)
for i, file in enumerate(inp_files):
dataset = read_subset(file, bbox)

# read/write *.unw file
f["unwrapPhase"][i] = dataset['unw_data']

# read/write *.cor file
f["coherence"][i] = dataset['cor_data']

# read/write *.unw.conncomp file
f["connectComponent"][i] = dataset['conn_comp']

# read/write perpendicular baseline file
f['bperp'][i] = dataset['pbase']

prog_bar.update(i + 1, suffix=date12_list[i])
elif "inputs/ionStack.h5" in outfile:
print("DOING IONO")
with h5py.File(outfile, "a") as f:
prog_bar = ptime.progressBar(maxValue=num_pair)
for i, file in enumerate(inp_files):
dataset = read_subset(file, bbox)

# read/write *.unw file
f["unwrapPhase"][i] = dataset['ion_data']

# read/write *.cor file
f["coherence"][i] = dataset['cor_data']

# read/write *.unw.conncomp file
f["connectComponent"][i] = dataset['conn_comp']

# read/write perpendicular baseline file
f['bperp'][i] = dataset['pbase']

with h5py.File(outfile, "a") as f:
prog_bar = ptime.progressBar(maxValue=num_pair)
for i, file in enumerate(inp_files):
dataset = read_subset(file, bbox)
prog_bar.update(i + 1, suffix=date12_list[i])
elif "inputs/tropoStack.h5" in outfile:
print("DOING TROPO")
print("NEED TO MASK THE DATA")
with h5py.File(outfile, "a") as f:
prog_bar = ptime.progressBar(maxValue=num_pair)
for i, file in enumerate(inp_files):
geo_ds = read_subset(inp_files[0], bbox, geometry=True)
total_tropo = read_and_interpolate_troposphere(inp_files[0], demFile, geo_ds['xybbox'])

# read/write *.unw file
f["unwrapPhase"][i] = dataset['unw_data']
f["unwrapPhase"][i] = total_tropo

# read/write *.cor file
f["coherence"][i] = dataset['cor_data']
# # read/write *.cor file
# f["coherence"][i] = dataset['cor_data']

# read/write *.unw.conncomp file
f["connectComponent"][i] = dataset['conn_comp']
# # read/write *.unw.conncomp file
# f["connectComponent"][i] = dataset['conn_comp']

# read/write perpendicular baseline file
f['bperp'][i] = dataset['pbase']
# # read/write perpendicular baseline file
# f['bperp'][i] = dataset['pbase']

prog_bar.update(i + 1, suffix=date12_list[i])
prog_bar.close()
prog_bar.update(i + 1, suffix=date12_list[i])
prog_bar.close()

print(f"finished writing to HDF5 file: {outfile}")
return outfile