Skip to content

Commit

Permalink
prep_nisar: update for changes in NISAR GUNW products (#1158)
Browse files Browse the repository at this point in the history 
* add mask to inputs

* remove output dir from prep_nisar

* update for mask and correct lat/lon subset

* add condition for other grib file pattern

* correct order of bounds

* use gdalwarp to transform projection

* change to new GUNW structure

---------

Co-authored-by: mirzaees <[email protected]>
  • Loading branch information
@mirzaees
mirzaees and mirzaees authored Mar 25, 2024
1 parent e7480d6 commit ccc946c
Show file tree
Hide file tree
Showing 4 changed files with 97 additions and 40 deletions.
11 changes: 10 additions & 1 deletion src/mintpy/cli/prep_nisar.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -40,12 +40,21 @@ def create_parser(subparsers=None):
help="path to the DEM (default: %(default)s).",
)

parser.add_argument(
"-m",
"--mask",
dest="mask_file",
type=str,
default=None,
help="path to the mask (default: %(default)s).",
)

parser.add_argument(
"-o",
"--out-dir",
dest="out_dir",
type=str,
default="./mintpy",
default=".",
help="output directory (default: %(default)s).",
)

Expand Down
6 changes: 5 additions & 1 deletion src/mintpy/load_data.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -634,9 +634,13 @@ def prepare_metadata(iDict):
elif processor == 'nisar':
dem_file = iDict['mintpy.load.demFile']
gunw_files = iDict['mintpy.load.unwFile']
water_mask = iDict['mintpy.load.waterMaskFile']

# run prep_*.py
iargs = ['-i', gunw_files, '-d', dem_file, '-o', '../mintpy']
iargs = ['-i', gunw_files, '-d', dem_file]

if os.path.exists(water_mask):
iargs = iargs + ['--mask', water_mask]

if iDict['mintpy.subset.yx']:
warnings.warn('Subset in Y/X is not implemented for NISAR. \n'
Expand Down
112 changes: 75 additions & 37 deletions src/mintpy/prep_nisar.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -18,23 +18,23 @@
from mintpy.constants import EARTH_RADIUS, SPEED_OF_LIGHT
from mintpy.utils import ptime, writefile

DATASET_ROOT_UNW = '/science/LSAR/GUNW/grids/frequencyA/interferogram/unwrapped'
DATASET_ROOT_UNW = '/science/LSAR/GUNW/grids/frequencyA/unwrappedInterferogram'
PARAMETERS = '/science/LSAR/GUNW/metadata/processingInformation/parameters/unwrappedInterferogram/frequencyA'
IDENTIFICATION = '/science/LSAR/identification'
RADARGRID_ROOT = 'science/LSAR/GUNW/metadata/radarGrid'
DATASETS = {
'xcoord' : f"{DATASET_ROOT_UNW}/xCoordinates",
'ycoord' : f"{DATASET_ROOT_UNW}/yCoordinates",
'xcoord' : f"{DATASET_ROOT_UNW}/POL/xCoordinates",
'ycoord' : f"{DATASET_ROOT_UNW}/POL/yCoordinates",
'unw' : f"{DATASET_ROOT_UNW}/POL/unwrappedPhase",
'cor' : f"{DATASET_ROOT_UNW}/POL/coherenceMagnitude",
'connComp' : f"{DATASET_ROOT_UNW}/POL/connectedComponents",
'layoverShadowMask': f"{DATASET_ROOT_UNW}/layoverShadowMask",
'waterMask' : f"{DATASET_ROOT_UNW}/waterMask",
#'mask' : f"{DATASET_ROOT_UNW}/mask",
'epsg' : f"{DATASET_ROOT_UNW}/projection",
'xSpacing' : f"{DATASET_ROOT_UNW}/xCoordinateSpacing",
'ySpacing' : f"{DATASET_ROOT_UNW}/yCoordinateSpacing",
'polarization' : f"{DATASET_ROOT_UNW}/listOfPolarizations",
'range_look' : f"{DATASET_ROOT_UNW}/numberOfRangeLooks",
'azimuth_look' : f"{DATASET_ROOT_UNW}/numberOfAzimuthLooks",
'polarization' : "/science/LSAR/GUNW/grids/frequencyA/listOfPolarizations",
'range_look' : f"{PARAMETERS}/numberOfRangeLooks",
'azimuth_look' : f"{PARAMETERS}/numberOfAzimuthLooks",
}
PROCESSINFO = {
'centerFrequency': "/science/LSAR/GUNW/grids/frequencyA/centerFrequency",
Expand All @@ -61,9 +61,9 @@ def load_nisar(inps):
print(f"Found {len(input_files)} unwrapped files")

if inps.subset_lat:
bbox = (inps.subset_lat[0], inps.subset_lon[0], inps.subset_lat[1], inps.subset_lon[1])
bbox = (inps.subset_lon[0], inps.subset_lat[0], inps.subset_lon[1], inps.subset_lat[1])
else:
bbox=None
bbox = None

# extract metadata
metadata, bounds = extract_metadata(input_files, bbox)
Expand All @@ -80,7 +80,8 @@ def load_nisar(inps):
metaFile=input_files[0],
bbox=bounds,
metadata=metadata,
demFile=inps.dem_file
demFile=inps.dem_file,
maskFile=inps.mask_file
)

prepare_stack(
Expand Down Expand Up @@ -136,7 +137,15 @@ def extract_metadata(input_files, bbox=None, polarization='HH'):
meta["Y_FIRST"] = y_origin - pixel_height//2
meta["X_STEP"] = pixel_width
meta["Y_STEP"] = pixel_height
meta["X_UNIT"] = meta["Y_UNIT"] = "meters"

if meta["EPSG"] == 4326:
meta["X_UNIT"] = meta["Y_UNIT"] = 'degree'
else:
meta["X_UNIT"] = meta["Y_UNIT"] = "meters"
if str(meta["EPSG"]).startswith('326'):
meta["UTM_ZONE"] = str(meta["EPSG"])[3:] + 'N'
else:
meta["UTM_ZONE"] = str(meta["EPSG"])[3:] + 'S'
meta["EARTH_RADIUS"] = EARTH_RADIUS

# NISAR Altitude
Expand All @@ -155,7 +164,6 @@ def extract_metadata(input_files, bbox=None, polarization='HH'):
utm_bbox = None
bounds = common_raster_bound(input_files, utm_bbox)
meta['bbox'] = ",".join([str(b) for b in bounds])

col1, row1, col2, row2 = get_rows_cols(xcoord, ycoord, bounds)
length = row2 - row1
width = col2 - col1
Expand Down Expand Up @@ -224,9 +232,7 @@ def read_subset(inp_file, bbox, geometry=False):
col1, row1, col2, row2 = get_rows_cols(xcoord, ycoord, bbox)

if geometry:
# Set all values to 1 temporarily because water mask is zero
dataset['water_mask'] = ds[DATASETS['waterMask']][row1:row2, col1:col2] * 0 + 1
dataset['layover_shadow_mask'] = ds[DATASETS['layoverShadowMask']][row1:row2, col1:col2]
# dataset['mask'] = ds[DATASETS['mask']][row1:row2, col1:col2]
dataset['xybbox'] = (col1, row1, col2, row2)
else:
dataset['unw_data'] = ds[DATASETS['unw']][row1:row2, col1:col2]
Expand All @@ -237,13 +243,14 @@ def read_subset(inp_file, bbox, geometry=False):
return dataset


def read_and_interpolate_geometry(gunw_file, dem_file, xybbox):
"""Read the DEM, change projection and interpolate to data grid, interpolate slant range and incidence angle to data grid"""
dataset = gdal.Open(dem_file, gdal.GA_ReadOnly)
geotransform = dataset.GetGeoTransform()
proj = gdal.osr.SpatialReference(wkt=dataset.GetProjection())
src_epsg = int(proj.GetAttrValue('AUTHORITY', 1))
raster_array = dataset.ReadAsArray()
def read_and_interpolate_geometry(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:
Expand All @@ -260,23 +267,50 @@ def read_and_interpolate_geometry(gunw_file, dem_file, xybbox):

subset_rows = len(ycoord)
subset_cols = len(xcoord)
# dem_subset_array = np.zeros((subset_rows, subset_cols), dtype=raster_array.dtype)

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}"

if not src_epsg == dst_epsg:
coord_transform = Transformer.from_crs(dst_epsg, src_epsg, always_xy=True)
x_dem, y_dem = coord_transform.transform(X_2d.flatten(), Y_2d.flatten())
# 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'])

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

# Interpolate slant_range and incidence_angle
slant_range, incidence_angle = interpolate_geometry(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:
x_dem, y_dem = X_2d.flatten(), Y_2d.flatten()
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

cols = ((y_dem - geotransform[3]) / geotransform[5]).astype(int)
rows = ((x_dem - geotransform[0]) / geotransform[1]).astype(int)
# 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'])

dem_subset_array = raster_array[cols.reshape(subset_rows, subset_cols), rows.reshape(subset_rows, subset_cols)]
mask_subset_array = gdal.Open(output_mask, gdal.GA_ReadOnly).ReadAsArray()

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

slant_range, incidence_angle = interpolate_geometry(X_2d, Y_2d, dem_subset_array, rdr_coords)

return dem_subset_array, slant_range, incidence_angle
return dem_subset_array, slant_range, incidence_angle, mask_subset_array


def interpolate_geometry(X_2d, Y_2d, dem, rdr_coords):
Expand Down Expand Up @@ -315,7 +349,8 @@ def prepare_geometry(
metaFile,
metadata,
bbox,
demFile
demFile,
maskFile
):
"""Prepare the geometry file."""
print("-" * 50)
Expand All @@ -326,17 +361,20 @@ def prepare_geometry(

# Read waterMask, LayoverShadowMask, xybbox:
geo_ds = read_subset(metaFile, bbox, geometry=True)
dem_subset_array, slant_range, incidence_angle = read_and_interpolate_geometry(metaFile, demFile, geo_ds['xybbox'])
dem_subset_array, slant_range, incidence_angle, mask = read_and_interpolate_geometry(metaFile, demFile,
geo_ds['xybbox'], mask_file=maskFile)

length, width = dem_subset_array.shape

ds_name_dict = {
"height": [np.float32, (length, width), dem_subset_array],
"incidenceAngle": [np.float32, (length, width), incidence_angle],
"slantRangeDistance": [np.float32, (length, width), slant_range],
"shadowMask": [np.bool_, (length, width), geo_ds['layover_shadow_mask']],
"waterMask": [np.bool_, (length, width), geo_ds['water_mask']],
#"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"] = "geometry"
Expand Down
8 changes: 7 additions & 1 deletion src/mintpy/tropo_pyaps3.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -483,7 +483,13 @@ def dload_grib_files(grib_files, tropo_model='ERA5', snwe=None):

# Download grib file using PyAPS
if len(date_list2dload) > 0:
hour = re.findall(r'\d{8}[-_]\d{2}', os.path.basename(grib_files2dload[0]))[0].replace('-', '_').split('_')[1]
# This was the default for a file pattern of "ERA5_N30_N50_W130_W110_20200403_14.grb"
pattern = re.findall(r'\d{8}[-_]\d{2}', os.path.basename(grib_files2dload[0]))
if len(pattern) == 0:
# This is an exception for when the file has a naming format like: "ERA5_N30_N40_W130_W110_20080125T000000_06.grb"
pattern = re.findall(r'\d{8}T\d{6}[-_]\d{2}', os.path.basename(grib_files2dload[0]))
hour = pattern[0].replace('-', '_').split('_')[1]
#hour = re.findall(r'\d{8}[-_]\d{2}', os.path.basename(grib_files2dload[0]))[0].replace('-', '_').split('_')[1]
grib_dir = os.path.dirname(grib_files2dload[0])

# Check for non-empty account info in PyAPS config file
Expand Down

0 comments on commit ccc946c

Please sign in to comment.