From 39a2d6619ef4340e6ebe3cc4e15e766856a629da Mon Sep 17 00:00:00 2001 From: mirzaees Date: Mon, 11 Mar 2024 10:32:25 -0700 Subject: [PATCH 01/11] add mask to inputs --- src/mintpy/cli/prep_nisar.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/mintpy/cli/prep_nisar.py b/src/mintpy/cli/prep_nisar.py index 67143e7e7..07004de28 100755 --- a/src/mintpy/cli/prep_nisar.py +++ b/src/mintpy/cli/prep_nisar.py @@ -41,12 +41,12 @@ def create_parser(subparsers=None): ) parser.add_argument( - "-o", - "--out-dir", - dest="out_dir", + "-m", + "--mask", + dest="mask_file", type=str, - default="./mintpy", - help="output directory (default: %(default)s).", + default=None, + help="path to the mask (default: %(default)s).", ) parser.add_argument( From d22e198029f798ebf3899ebc07f159123e452214 Mon Sep 17 00:00:00 2001 From: mirzaees Date: Mon, 11 Mar 2024 10:33:04 -0700 Subject: [PATCH 02/11] remove output dir from prep_nisar --- src/mintpy/load_data.py | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/src/mintpy/load_data.py b/src/mintpy/load_data.py index 01ec011a2..fd10cfcc2 100644 --- a/src/mintpy/load_data.py +++ b/src/mintpy/load_data.py @@ -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' From 98ccb11ca6d79ed702f22cbff956ea36c395eeeb Mon Sep 17 00:00:00 2001 From: mirzaees Date: Mon, 11 Mar 2024 10:33:43 -0700 Subject: [PATCH 03/11] update for mask and correct lat/lon subset --- src/mintpy/prep_nisar.py | 93 +++++++++++++++++++++++++++++----------- 1 file changed, 68 insertions(+), 25 deletions(-) diff --git a/src/mintpy/prep_nisar.py b/src/mintpy/prep_nisar.py index 12f266037..b9b896806 100644 --- a/src/mintpy/prep_nisar.py +++ b/src/mintpy/prep_nisar.py @@ -27,8 +27,8 @@ '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", + # 'waterMask' : f"{DATASET_ROOT_UNW}/mask", 'epsg' : f"{DATASET_ROOT_UNW}/projection", 'xSpacing' : f"{DATASET_ROOT_UNW}/xCoordinateSpacing", 'ySpacing' : f"{DATASET_ROOT_UNW}/yCoordinateSpacing", @@ -69,8 +69,8 @@ def load_nisar(inps): metadata, bounds = extract_metadata(input_files, bbox) # output filename - stack_file = os.path.join(inps.out_dir, "inputs/ifgramStack.h5") - geometry_file = os.path.join(inps.out_dir, "inputs/geometryGeo.h5") + stack_file = os.path.join("./inputs/ifgramStack.h5") + geometry_file = os.path.join("./inputs/geometryGeo.h5") # get date info: date12_list date12_list = _get_date_pairs(input_files) @@ -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( @@ -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 @@ -195,8 +204,8 @@ def common_raster_bound(input_files, utm_bbox=None): x_bounds.append([west, east]) y_bounds.append([south, north]) if not utm_bbox is None: - x_bounds.append([utm_bbox[0], utm_bbox[2]]) - y_bounds.append([utm_bbox[1], utm_bbox[3]]) + x_bounds.append([utm_bbox[1], utm_bbox[3]]) + y_bounds.append([utm_bbox[0], utm_bbox[2]]) bounds = max(x_bounds)[0], max(y_bounds)[0], min(x_bounds)[1], min(y_bounds)[1] return bounds @@ -225,8 +234,8 @@ def read_subset(inp_file, bbox, geometry=False): 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['water_mask'] = ds[DATASETS['waterMask']][row1:row2, col1:col2] * 0 + 1 + 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] @@ -237,13 +246,16 @@ 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) + geotransform = dem_dataset.GetGeoTransform() + proj = gdal.osr.SpatialReference(wkt=dem_dataset.GetProjection()) + dem_src_epsg = int(proj.GetAttrValue('AUTHORITY', 1)) + dem_raster_array = dem_dataset.ReadAsArray() + del dem_dataset + rdr_coords = {} with h5py.File(gunw_file, 'r') as ds: @@ -263,8 +275,8 @@ def read_and_interpolate_geometry(gunw_file, dem_file, xybbox): # dem_subset_array = np.zeros((subset_rows, subset_cols), dtype=raster_array.dtype) Y_2d, X_2d = np.meshgrid(ycoord, xcoord, indexing='ij') - if not src_epsg == dst_epsg: - coord_transform = Transformer.from_crs(dst_epsg, src_epsg, always_xy=True) + if not dem_src_epsg == dst_epsg: + coord_transform = Transformer.from_crs(dst_epsg, dem_src_epsg, always_xy=True) x_dem, y_dem = coord_transform.transform(X_2d.flatten(), Y_2d.flatten()) else: x_dem, y_dem = X_2d.flatten(), Y_2d.flatten() @@ -272,11 +284,38 @@ def read_and_interpolate_geometry(gunw_file, dem_file, xybbox): cols = ((y_dem - geotransform[3]) / geotransform[5]).astype(int) rows = ((x_dem - geotransform[0]) / geotransform[1]).astype(int) - dem_subset_array = raster_array[cols.reshape(subset_rows, subset_cols), rows.reshape(subset_rows, subset_cols)] + dem_subset_array = dem_raster_array[cols.reshape(subset_rows, subset_cols), rows.reshape(subset_rows, subset_cols)] slant_range, incidence_angle = interpolate_geometry(X_2d, Y_2d, dem_subset_array, rdr_coords) - return dem_subset_array, slant_range, incidence_angle + + 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) + geotransform = mask_dataset.GetGeoTransform() + proj = gdal.osr.SpatialReference(wkt=mask_dataset.GetProjection()) + mask_src_epsg = int(proj.GetAttrValue('AUTHORITY', 1)) + mask_raster_array = mask_dataset.ReadAsArray() + del mask_dataset + + if not mask_src_epsg == dst_epsg: + coord_transform = Transformer.from_crs(dst_epsg, mask_src_epsg, always_xy=True) + x_mask, y_mask = coord_transform.transform(X_2d.flatten(), Y_2d.flatten()) + else: + x_mask, y_mask = X_2d.flatten(), Y_2d.flatten() + + cols = ((y_mask - geotransform[3]) / geotransform[5]).astype(int) + rows = ((x_mask - geotransform[0]) / geotransform[1]).astype(int) + mask_subset_array = mask_raster_array[cols.reshape(subset_rows, subset_cols), + rows.reshape(subset_rows, subset_cols)] + except: + raise IOError('*** Mask is not gdal readable ***') + + + return dem_subset_array, slant_range, incidence_angle, mask_subset_array def interpolate_geometry(X_2d, Y_2d, dem, rdr_coords): @@ -315,7 +354,8 @@ def prepare_geometry( metaFile, metadata, bbox, - demFile + demFile, + maskFile ): """Prepare the geometry file.""" print("-" * 50) @@ -326,7 +366,8 @@ 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 @@ -334,9 +375,11 @@ def prepare_geometry( "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" From ba10fab4ec4e1bf4b886727c4572db643b7eacad Mon Sep 17 00:00:00 2001 From: mirzaees Date: Mon, 11 Mar 2024 10:35:18 -0700 Subject: [PATCH 04/11] add condition for other grib file pattern --- src/mintpy/tropo_pyaps3.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/mintpy/tropo_pyaps3.py b/src/mintpy/tropo_pyaps3.py index 49867cd0d..de5571b70 100644 --- a/src/mintpy/tropo_pyaps3.py +++ b/src/mintpy/tropo_pyaps3.py @@ -419,7 +419,11 @@ 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] + pattern = re.findall(r'\d{8}[-_]\d{2}', os.path.basename(grib_files2dload[0])) + if len(pattern) == 0: + 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 From 1b2bb4759c9af462dd2c9d097bc3f2e486e2be70 Mon Sep 17 00:00:00 2001 From: mirzaees Date: Mon, 11 Mar 2024 10:41:55 -0700 Subject: [PATCH 05/11] pre-commit changes --- src/mintpy/load_data.py | 2 +- src/mintpy/prep_nisar.py | 14 +++++++------- 2 files changed, 8 insertions(+), 8 deletions(-) diff --git a/src/mintpy/load_data.py b/src/mintpy/load_data.py index fd10cfcc2..5ef0bb842 100644 --- a/src/mintpy/load_data.py +++ b/src/mintpy/load_data.py @@ -635,7 +635,7 @@ def prepare_metadata(iDict): 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] diff --git a/src/mintpy/prep_nisar.py b/src/mintpy/prep_nisar.py index b9b896806..c0cfe5575 100644 --- a/src/mintpy/prep_nisar.py +++ b/src/mintpy/prep_nisar.py @@ -247,9 +247,9 @@ def read_subset(inp_file, bbox, geometry=False): 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, + """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) + dem_dataset = gdal.Open(dem_file, gdal.GA_ReadOnly) geotransform = dem_dataset.GetGeoTransform() proj = gdal.osr.SpatialReference(wkt=dem_dataset.GetProjection()) dem_src_epsg = int(proj.GetAttrValue('AUTHORITY', 1)) @@ -306,14 +306,14 @@ def read_and_interpolate_geometry(gunw_file, dem_file, xybbox, mask_file=None): x_mask, y_mask = coord_transform.transform(X_2d.flatten(), Y_2d.flatten()) else: x_mask, y_mask = X_2d.flatten(), Y_2d.flatten() - + cols = ((y_mask - geotransform[3]) / geotransform[5]).astype(int) rows = ((x_mask - geotransform[0]) / geotransform[1]).astype(int) - mask_subset_array = mask_raster_array[cols.reshape(subset_rows, subset_cols), + mask_subset_array = mask_raster_array[cols.reshape(subset_rows, subset_cols), rows.reshape(subset_rows, subset_cols)] except: - raise IOError('*** Mask is not gdal readable ***') - + raise OSError('*** Mask is not gdal readable ***') + return dem_subset_array, slant_range, incidence_angle, mask_subset_array @@ -366,7 +366,7 @@ def prepare_geometry( # Read waterMask, LayoverShadowMask, xybbox: geo_ds = read_subset(metaFile, bbox, geometry=True) - dem_subset_array, slant_range, incidence_angle, mask = read_and_interpolate_geometry(metaFile, demFile, + 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 From c1290230629f2fa985edc2e25167dda367798d20 Mon Sep 17 00:00:00 2001 From: mirzaees Date: Tue, 12 Mar 2024 11:17:17 -0700 Subject: [PATCH 06/11] correct order of bounds --- src/mintpy/prep_nisar.py | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/src/mintpy/prep_nisar.py b/src/mintpy/prep_nisar.py index c0cfe5575..a18332eaa 100644 --- a/src/mintpy/prep_nisar.py +++ b/src/mintpy/prep_nisar.py @@ -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) @@ -164,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 @@ -204,8 +203,8 @@ def common_raster_bound(input_files, utm_bbox=None): x_bounds.append([west, east]) y_bounds.append([south, north]) if not utm_bbox is None: - x_bounds.append([utm_bbox[1], utm_bbox[3]]) - y_bounds.append([utm_bbox[0], utm_bbox[2]]) + x_bounds.append([utm_bbox[0], utm_bbox[2]]) + y_bounds.append([utm_bbox[1], utm_bbox[3]]) bounds = max(x_bounds)[0], max(y_bounds)[0], min(x_bounds)[1], min(y_bounds)[1] return bounds From ee9696906afd6b947ea2d4989e7a2f842eacdb1a Mon Sep 17 00:00:00 2001 From: mirzaees Date: Tue, 12 Mar 2024 11:17:38 -0700 Subject: [PATCH 07/11] add comment --- src/mintpy/tropo_pyaps3.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/mintpy/tropo_pyaps3.py b/src/mintpy/tropo_pyaps3.py index de5571b70..9b50b7545 100644 --- a/src/mintpy/tropo_pyaps3.py +++ b/src/mintpy/tropo_pyaps3.py @@ -419,8 +419,10 @@ def dload_grib_files(grib_files, tropo_model='ERA5', snwe=None): # Download grib file using PyAPS if len(date_list2dload) > 0: + # 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] From aa8231f705691521cb409ad4881362eab1abb5d0 Mon Sep 17 00:00:00 2001 From: mirzaees Date: Wed, 13 Mar 2024 14:47:20 -0700 Subject: [PATCH 08/11] use gdalwarp to transform projection --- src/mintpy/cli/prep_nisar.py | 9 ++++++ src/mintpy/prep_nisar.py | 63 +++++++++++++++++------------------- 2 files changed, 38 insertions(+), 34 deletions(-) diff --git a/src/mintpy/cli/prep_nisar.py b/src/mintpy/cli/prep_nisar.py index 07004de28..8068f8f3b 100755 --- a/src/mintpy/cli/prep_nisar.py +++ b/src/mintpy/cli/prep_nisar.py @@ -49,6 +49,15 @@ def create_parser(subparsers=None): help="path to the mask (default: %(default)s).", ) + parser.add_argument( + "-o", + "--out-dir", + dest="out_dir", + type=str, + default=".", + help="output directory (default: %(default)s).", + ) + parser.add_argument( '--force', dest='update_mode', diff --git a/src/mintpy/prep_nisar.py b/src/mintpy/prep_nisar.py index a18332eaa..ded8f4ed7 100644 --- a/src/mintpy/prep_nisar.py +++ b/src/mintpy/prep_nisar.py @@ -28,7 +28,6 @@ 'cor' : f"{DATASET_ROOT_UNW}/POL/coherenceMagnitude", 'connComp' : f"{DATASET_ROOT_UNW}/POL/connectedComponents", 'mask' : f"{DATASET_ROOT_UNW}/mask", - # 'waterMask' : f"{DATASET_ROOT_UNW}/mask", 'epsg' : f"{DATASET_ROOT_UNW}/projection", 'xSpacing' : f"{DATASET_ROOT_UNW}/xCoordinateSpacing", 'ySpacing' : f"{DATASET_ROOT_UNW}/yCoordinateSpacing", @@ -69,8 +68,8 @@ def load_nisar(inps): metadata, bounds = extract_metadata(input_files, bbox) # output filename - stack_file = os.path.join("./inputs/ifgramStack.h5") - geometry_file = os.path.join("./inputs/geometryGeo.h5") + stack_file = os.path.join(inps.out_dir, "inputs/ifgramStack.h5") + geometry_file = os.path.join(inps.out_dir, "inputs/geometryGeo.h5") # get date info: date12_list date12_list = _get_date_pairs(input_files) @@ -232,8 +231,6 @@ 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['mask'] = ds[DATASETS['mask']][row1:row2, col1:col2] dataset['xybbox'] = (col1, row1, col2, row2) else: @@ -249,10 +246,8 @@ 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) - geotransform = dem_dataset.GetGeoTransform() proj = gdal.osr.SpatialReference(wkt=dem_dataset.GetProjection()) dem_src_epsg = int(proj.GetAttrValue('AUTHORITY', 1)) - dem_raster_array = dem_dataset.ReadAsArray() del dem_dataset rdr_coords = {} @@ -271,45 +266,45 @@ def read_and_interpolate_geometry(gunw_file, dem_file, xybbox, mask_file=None): 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') - - if not dem_src_epsg == dst_epsg: - coord_transform = Transformer.from_crs(dst_epsg, dem_src_epsg, always_xy=True) - x_dem, y_dem = coord_transform.transform(X_2d.flatten(), Y_2d.flatten()) - else: - x_dem, y_dem = X_2d.flatten(), Y_2d.flatten() - - cols = ((y_dem - geotransform[3]) / geotransform[5]).astype(int) - rows = ((x_dem - geotransform[0]) / geotransform[1]).astype(int) - - dem_subset_array = dem_raster_array[cols.reshape(subset_rows, subset_cols), rows.reshape(subset_rows, subset_cols)] - + 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' ) + 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: try: mask_dataset = gdal.Open(mask_file, gdal.GA_ReadOnly) - geotransform = mask_dataset.GetGeoTransform() proj = gdal.osr.SpatialReference(wkt=mask_dataset.GetProjection()) mask_src_epsg = int(proj.GetAttrValue('AUTHORITY', 1)) - mask_raster_array = mask_dataset.ReadAsArray() del mask_dataset - if not mask_src_epsg == dst_epsg: - coord_transform = Transformer.from_crs(dst_epsg, mask_src_epsg, always_xy=True) - x_mask, y_mask = coord_transform.transform(X_2d.flatten(), Y_2d.flatten()) - else: - x_mask, y_mask = X_2d.flatten(), Y_2d.flatten() + # 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() - cols = ((y_mask - geotransform[3]) / geotransform[5]).astype(int) - rows = ((x_mask - geotransform[0]) / geotransform[1]).astype(int) - mask_subset_array = mask_raster_array[cols.reshape(subset_rows, subset_cols), - rows.reshape(subset_rows, subset_cols)] except: raise OSError('*** Mask is not gdal readable ***') @@ -447,9 +442,9 @@ def prepare_stack( # read/write perpendicular baseline file f['bperp'][i] = dataset['pbase'] - + prog_bar.update(i + 1, suffix=date12_list[i]) prog_bar.close() - + print(f"finished writing to HDF5 file: {outfile}") return outfile From ec1dbb5df46cda5e0620ba2928dec26178136ee4 Mon Sep 17 00:00:00 2001 From: mirzaees Date: Wed, 13 Mar 2024 14:47:20 -0700 Subject: [PATCH 09/11] use gdalwarp to transform projection --- src/mintpy/prep_nisar.py | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/src/mintpy/prep_nisar.py b/src/mintpy/prep_nisar.py index ded8f4ed7..cdf64ee2b 100644 --- a/src/mintpy/prep_nisar.py +++ b/src/mintpy/prep_nisar.py @@ -266,16 +266,16 @@ def read_and_interpolate_geometry(gunw_file, dem_file, xybbox, mask_file=None): 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 + + # 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, + srcSRS=input_projection, dstSRS=output_projection, resampleAlg=gdal.GRA_Bilinear, width=subset_cols, height=subset_rows, options=['COMPRESS=DEFLATE']) @@ -295,11 +295,11 @@ def read_and_interpolate_geometry(gunw_file, dem_file, xybbox, mask_file=None): mask_src_epsg = int(proj.GetAttrValue('AUTHORITY', 1)) del mask_dataset - # Warp mask to the interferograms grid + # 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, + srcSRS=input_projection, dstSRS=output_projection, resampleAlg=gdal.GRA_Byte, width=subset_cols, height=subset_rows, options=['COMPRESS=DEFLATE']) @@ -442,9 +442,9 @@ def prepare_stack( # read/write perpendicular baseline file f['bperp'][i] = dataset['pbase'] - + prog_bar.update(i + 1, suffix=date12_list[i]) prog_bar.close() - + print(f"finished writing to HDF5 file: {outfile}") return outfile From b1976d2c5f64d61c5243867253c4a10f60446735 Mon Sep 17 00:00:00 2001 From: mirzaees Date: Fri, 22 Mar 2024 11:14:42 -0700 Subject: [PATCH 10/11] change to new GUNW structure --- src/mintpy/prep_nisar.py | 19 ++++++++++--------- 1 file changed, 10 insertions(+), 9 deletions(-) diff --git a/src/mintpy/prep_nisar.py b/src/mintpy/prep_nisar.py index cdf64ee2b..925bdc022 100644 --- a/src/mintpy/prep_nisar.py +++ b/src/mintpy/prep_nisar.py @@ -18,22 +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", - 'mask' : f"{DATASET_ROOT_UNW}/mask", + #'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' : f"/science/LSAR/GUNW/grids/frequencyA/listOfPolarizations", + 'range_look' : f"{PARAMETERS}/numberOfRangeLooks", + 'azimuth_look' : f"{PARAMETERS}/numberOfAzimuthLooks", } PROCESSINFO = { 'centerFrequency': "/science/LSAR/GUNW/grids/frequencyA/centerFrequency", @@ -231,7 +232,7 @@ def read_subset(inp_file, bbox, geometry=False): col1, row1, col2, row2 = get_rows_cols(xcoord, ycoord, bbox) if geometry: - dataset['mask'] = ds[DATASETS['mask']][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] @@ -369,7 +370,7 @@ def prepare_geometry( "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['mask']], + #"shadowMask": [np.bool_, (length, width), geo_ds['mask']], #"waterMask": [np.bool_, (length, width), geo_ds['water_mask']], } if maskFile: From 0ec622d1d49d52f75651ea42bc4763b90eeb40e1 Mon Sep 17 00:00:00 2001 From: mirzaees Date: Mon, 25 Mar 2024 07:42:49 -0700 Subject: [PATCH 11/11] remove f string --- src/mintpy/prep_nisar.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/mintpy/prep_nisar.py b/src/mintpy/prep_nisar.py index 925bdc022..f9cd6b6f7 100644 --- a/src/mintpy/prep_nisar.py +++ b/src/mintpy/prep_nisar.py @@ -32,7 +32,7 @@ 'epsg' : f"{DATASET_ROOT_UNW}/projection", 'xSpacing' : f"{DATASET_ROOT_UNW}/xCoordinateSpacing", 'ySpacing' : f"{DATASET_ROOT_UNW}/yCoordinateSpacing", - 'polarization' : f"/science/LSAR/GUNW/grids/frequencyA/listOfPolarizations", + 'polarization' : "/science/LSAR/GUNW/grids/frequencyA/listOfPolarizations", 'range_look' : f"{PARAMETERS}/numberOfRangeLooks", 'azimuth_look' : f"{PARAMETERS}/numberOfAzimuthLooks", }