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

Jump to bottom

Create Mintpy geometryGeo.h5 file in prep_mintpy.py #24

Merged
merged 3 commits into from
Aug 28, 2023
Merged

Create Mintpy geometryGeo.h5 file in prep_mintpy.py #24

Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
440014f
create a `geometryGeo.h5` file in mintpy script
scottstanie Aug 28, 2023
0d8518a
add a try/except for pending mintpy int8 support
scottstanie Aug 28, 2023
15ac037
add azimuth/incidence datasets as redundancy
scottstanie Aug 28, 2023
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
105 changes: 40 additions & 65 deletions scripts/prep_mintpy.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -20,37 +20,21 @@
from dolphin import io
from dolphin.utils import full_suffix, get_dates
from dolphin.workflows.config import OPERA_DATASET_ROOT
from mintpy.utils import arg_utils, isce_utils, ptime, readfile, writefile
from mintpy.utils import arg_utils, ptime, readfile, writefile
from mintpy.utils.utils0 import calc_azimuth_from_east_north_obs

####################################################################################
EXAMPLE = """example:

python3 ./prep_sweets.py -u 'interferograms/stitched/*.unw' -m gslcs/t064_135516_iw3/20190704/static_layers_t064_135516_iw3.h5
python ./prep_mintpy.py -u 'interferograms/unwrapped/*.unw.tif' -c 'interferograms/stitched/*.cor' -m gslcs/

## example commands after prep_sweets.py
reference_point.py timeseries.h5 -y 500 -x 1150
generate_mask.py temporalCoherence.h5 -m 0.7 -o maskTempCoh.h5
tropo_pyaps3.py -f timeseries.h5 -g inputs/geometryRadar.h5
remove_ramp.py timeseries_ERA5.h5 -m maskTempCoh.h5 -s linear
dem_error.py timeseries_ERA5_ramp.h5 -g inputs/geometryRadar.h5
timeseries2velocity.py timeseries_ERA5_ramp_demErr.h5
geocode.py velocity.h5 -l inputs/geometryRadar.h5
""" # noqa: E501

# """
# Scott TODO:
# - dont need the ifg sample file, the unwrapped files should have everything
# - UTM_ZONE, EPSG from the stitched IFG (it won't work to get a single GSLC burst)
# - pixel size is wrong since we're taking range/azimuth size, instead of geocoded size
# - HEIGHT: do we wanna try to get that from the saved orbit info?
# - cor_files = [x.split(".unw")[0] + ".cor" for x in unw_files]
# - This might be wrong for dolphin? to verify

# Scott did:
# - remove bbox
# - "set all data less than 0" thats a fringe specific thing to do -1, -2 for tcoh

# """


def _create_parser():
Expand Down Expand Up @@ -83,7 +67,6 @@ def _create_parser():
parser.add_argument(
"-m",
"--meta-file",
dest="metaFile",
type=str,
help="GSLC metadata file or directory",
)
Expand All @@ -98,7 +81,6 @@ def _create_parser():
parser.add_argument(
"-o",
"--out-dir",
dest="outDir",
type=str,
default="./mintpy",
help="output directory (default: %(default)s).",
Expand Down Expand Up @@ -134,7 +116,7 @@ def _create_parser():
),
)

parser = arg_utils.add_subset_argument(parser, geo=False)
parser = arg_utils.add_subset_argument(parser, geo=True)

return parser

Expand All @@ -145,14 +127,11 @@ def cmd_line_parse(iargs=None):
inps = parser.parse_args(args=iargs)

# in case meta_file is input as wildcard
inps.metaFile = sorted(glob.glob(inps.metaFile))[0]
inps.meta_file = sorted(glob.glob(inps.meta_file))[0]

return inps


####################################################################################


def prepare_metadata(meta_file, int_file, nlks_x=1, nlks_y=1):
"""Get the metadata from the GSLC metadata file and the unwrapped interferogram."""
print("-" * 50)
Expand Down Expand Up @@ -327,7 +306,6 @@ def prepare_timeseries(
outfile,
unw_files,
metadata,
processor,
baseline_dir=None,
):
"""Prepare the timeseries file."""
Expand All @@ -353,16 +331,7 @@ def prepare_timeseries(
# baseline info
pbase = np.zeros(num_date, dtype=np.float32)
if baseline_dir is not None:
# read baseline data
# TODO: this isn't implemented yet
raise NotImplementedError
# baseline_dict = isce_utils.read_baseline_timeseries(
# baseline_dir, processor=processor, ref_date=ref_date
# )
# # dict to array
# for i in range(num_date):
# pbase_top, pbase_bottom = baseline_dict[date_list[i]]
# pbase[i] = (pbase_top + pbase_bottom) / 2.0

# size info
cols, rows = io.get_raster_xysize(unw_files[0])
Expand Down Expand Up @@ -400,7 +369,7 @@ def prepare_timeseries(
return outfile


def prepare_geometry(outfile, geom_dir, metadata, box, water_mask_file=None):
def prepare_geometry(outfile, geom_dir, metadata, water_mask_file=None):
"""Prepare the geometry file."""
print("-" * 50)
print(f"preparing geometry file: {outfile}")
Expand All @@ -411,9 +380,9 @@ def prepare_geometry(outfile, geom_dir, metadata, box, water_mask_file=None):
meta["FILE_TYPE"] = "geometry"

file_to_path = {
"los_east": geom_path / "los_east.tif",
"los_north": geom_path / "los_north.tif",
"height": geom_path / "height.tif",
"incidenceAngle": geom_path / "incidence_angle.tif",
"azimuthAngle": geom_path / "heading_angle.tif",
"shadowMask": geom_path / "layover_shadow_mask.tif",
}

Expand All @@ -422,11 +391,27 @@ def prepare_geometry(outfile, geom_dir, metadata, box, water_mask_file=None):

dsDict = {}
for dsName, fname in file_to_path.items():
dsDict[dsName] = readfile.read(fname, datasetName=dsName, box=box)[0]
try:
data = readfile.read(fname, datasetName=dsName)[0]
# TODO: add general functionality to handle nodata into Mintpy
data[data == 0] = np.nan
dsDict[dsName] = data

# write data to HDF5 file
writefile.write(dsDict, outfile, metadata=meta)
# write data to HDF5 file
except KeyError as e: # https://github.com/insarlab/MintPy/issues/1081
print(f"Skipping {fname}: {e}")

# Compute the azimuth and incidence angles from east/north coefficients
east = dsDict["los_east"]
north = dsDict["los_north"]
azimuth_angle = calc_azimuth_from_east_north_obs(east, north)
dsDict["azimuthAngle"] = azimuth_angle

up = np.sqrt(1 - east**2 - north**2)
incidence_angle = np.rad2deg(np.arccos(up))
dsDict["incidenceAngle"] = incidence_angle

writefile.write(dsDict, outfile, metadata=meta)
return outfile


Expand Down Expand Up @@ -471,8 +456,7 @@ def prepare_stack(
unw_files,
cor_files,
metadata,
processor,
baseline_dir=None,
# baseline_dir=None,
):
"""Prepare the input unw stack."""
print("-" * 50)
Expand Down Expand Up @@ -508,20 +492,8 @@ def prepare_stack(
# get date info: date12_list
date12_list = _get_date_pairs(unw_files)

# TODO: compute the spatial baseline using COMPASS metadata
pbase = np.zeros(num_pair, dtype=np.float32)
# prepare baseline info
if baseline_dir is not None:
# read baseline timeseries
baseline_dict = isce_utils.read_baseline_timeseries(
baseline_dir, processor=processor
)

# calc baseline for each pair
print("calc perp baseline pairs from time-series")
# pbase = np.zeros(num_pair, dtype=np.float32)
for i, date12 in enumerate(date12_list):
[date1, date2] = date12.split("_")
pbase[i] = np.subtract(baseline_dict[date2], baseline_dict[date1]).mean()

# size info
cols, rows = io.get_raster_xysize(unw_files[0])
Expand Down Expand Up @@ -577,10 +549,11 @@ def main(iargs=None):
print(f"Found {len(unw_files)} unwrapped files")
cor_files = sorted(glob.glob(inps.cor_file_glob))
print(f"Found {len(cor_files)} correlation files")

# translate input options
processor = "sweets" # isce_utils.get_processor(inps.metaFile)
processor = "sweets" # isce_utils.get_processor(inps.meta_file)
# metadata
meta_file = Path(inps.metaFile)
meta_file = Path(inps.meta_file)
if meta_file.is_dir():
# Search for the line of sight static_layers file
try:
Expand All @@ -594,11 +567,12 @@ def main(iargs=None):
)

# output directory
for dname in [inps.outDir, os.path.join(inps.outDir, "inputs")]:
for dname in [inps.out_dir, os.path.join(inps.out_dir, "inputs")]:
os.makedirs(dname, exist_ok=True)

stack_file = os.path.join(inps.outDir, "inputs/ifgramStack.h5")
ts_file = os.path.join(inps.outDir, "timeseries.h5")
stack_file = os.path.join(inps.out_dir, "inputs/ifgramStack.h5")
ts_file = os.path.join(inps.out_dir, "timeseries.h5")
geom_file = os.path.join(inps.out_dir, "geometryGeo.h5")

if inps.single_reference:
# time-series (if inputs are all single-reference)
Expand All @@ -607,17 +581,18 @@ def main(iargs=None):
unw_files=unw_files,
metadata=meta,
processor=processor,
baseline_dir=inps.baselineDir,
# baseline_dir=inps.baseline_dir,
)

prepare_geometry(geom_file, geom_dir=inps.geom_dir, metadata=meta)

# prepare ifgstack with connected components
prepare_stack(
outfile=stack_file,
unw_files=unw_files,
cor_files=cor_files,
metadata=meta,
processor=processor,
baseline_dir=inps.baselineDir,
# baseline_dir=inps.baseline_dir,
)

print("Done.")
Expand Down
Loading