From febf7c5ff3071da56ed3d1a81ea4a14f19339b76 Mon Sep 17 00:00:00 2001 From: Falk Amelung Date: Thu, 25 Jul 2024 18:39:38 +0200 Subject: [PATCH] subset.py: support of --lat --lon in radar coordinates --- src/mintpy/cli/subset.py | 2 ++ src/mintpy/objects/coord.py | 8 +++-- src/mintpy/subset.py | 60 ++++++++++++++++++++++++++----------- 3 files changed, 51 insertions(+), 19 deletions(-) diff --git a/src/mintpy/cli/subset.py b/src/mintpy/cli/subset.py index d53546008..18aa974a8 100755 --- a/src/mintpy/cli/subset.py +++ b/src/mintpy/cli/subset.py @@ -20,6 +20,7 @@ EXAMPLE = """example: subset.py inputs/ifgramStack.h5 -y 400 1500 -x 200 600 subset.py geo_velocity.h5 -l 30.5 30.8 -L 130.3 130.9 + subset.py velocity.h5 --lat 25.894 25.896 --lon -80.124 -80.122 subset.py 030405_090801.unw -t SinabungT495F50AlosA.template subset.py demLat*.dem.wgs84 --lat 32.5 33.0 --lon 130.2 130.6 -o srtm1.h5 @@ -28,6 +29,7 @@ # multiple files input subset.py *velocity*.h5 timeseries*.h5 -y 400 1500 -x 200 600 + subset.py geometryRadar.h5 slcStack.h5 --lat 25.894 25.896 --lon -80.124 -80.122 # crop to larger area with custom fill value subset.py geo_velocity.h5 -l 32.2 33.5 --outfill-nan diff --git a/src/mintpy/objects/coord.py b/src/mintpy/objects/coord.py index 2b8a23195..46eb334bc 100644 --- a/src/mintpy/objects/coord.py +++ b/src/mintpy/objects/coord.py @@ -490,16 +490,20 @@ def bbox_radar2geo(self, pix_box, print_msg=False): return geo_box - def bbox_geo2radar(self, geo_box, print_msg=False): + def bbox_geo2radar(self, geo_box, buf=None, print_msg=False): """Calculate bounding box in x/y for file in radar coord, based on input geo box. Parameters: geo_box - tuple of 4 float, indicating the UL/LR lon/lat + buf - buffer [number of pixel] Returns: pix_box - tuple of 4 int, indicating the UL/LR x/y of the bounding box in radar coord for the corresponding lat/lon coverage. """ lat = np.array([geo_box[3], geo_box[3], geo_box[1], geo_box[1]]) lon = np.array([geo_box[0], geo_box[2], geo_box[0], geo_box[2]]) y, x, y_res, x_res = self.geo2radar(lat, lon, print_msg=print_msg) - buf = 2 * np.max(np.abs([x_res, y_res])) + + if buf is None: + buf = 2 * np.max(np.abs([x_res, y_res])) + pix_box = (np.min(x) - buf, np.min(y) - buf, np.max(x) + buf, np.max(y) + buf) return pix_box diff --git a/src/mintpy/subset.py b/src/mintpy/subset.py index fe6689fee..d7556f718 100644 --- a/src/mintpy/subset.py +++ b/src/mintpy/subset.py @@ -183,27 +183,53 @@ def subset_input_dict2box(subset_dict, meta_dict): # Use subset_lat/lon input if existed, priority: lat/lon > y/x > len/wid coord = ut.coordinate(meta_dict) - if subset_dict.get('subset_lat', None): - sub_y = coord.lalo2yx(subset_dict['subset_lat'], None)[0] - elif subset_dict['subset_y']: - sub_y = subset_dict['subset_y'] - else: - sub_y = [0, length] - if subset_dict.get('subset_lon', None): - sub_x = coord.lalo2yx(None, subset_dict['subset_lon'])[1] - elif subset_dict['subset_x']: - sub_x = subset_dict['subset_x'] + if 'Y_FIRST' in meta_dict.keys(): + geocoded_flag =True else: - sub_x = [0, width] + geocoded_flag =False - # Get subset box in y/x - sub_x = sorted(sub_x) - sub_y = sorted(sub_y) - pix_box = (sub_x[0], sub_y[0], sub_x[1], sub_y[1]) + if geocoded_flag: + if subset_dict.get('subset_lat', None): + sub_y = coord.lalo2yx(subset_dict['subset_lat'], None)[0] + else: + sub_y = [0, length] - # Get subset box in lat/lon from subset box in y/x - geo_box = coord.box_pixel2geo(pix_box) + if subset_dict.get('subset_lon', None): + sub_x = coord.lalo2yx(None, subset_dict['subset_lon'])[1] + else: + sub_x = [0, width] + + # Get subset box in y/x + sub_x = sorted(sub_x) + sub_y = sorted(sub_y) + pix_box = (sub_x[0], sub_y[0], sub_x[1], sub_y[1]) + + # Get subset box in lat/lon from subset box in y/x + geo_box = coord.box_pixel2geo(pix_box) + + # for data in radar coordinates + else: + if subset_dict.get('subset_lat', None): + geo_box = (subset_dict['subset_lon'][0], subset_dict['subset_lat'][1], subset_dict['subset_lon'][1], subset_dict['subset_lat'][0]) + pix_box = coord.bbox_geo2radar(geo_box, buf=0) + else: + if subset_dict['subset_y']: + sub_y = subset_dict['subset_y'] + else: + sub_y = [0, length] + + if subset_dict['subset_x']: + sub_x = subset_dict['subset_x'] + else: + sub_x = [0, width] + # Get subset box in y/x + sub_x = sorted(sub_x) + sub_y = sorted(sub_y) + pix_box = (sub_x[0], sub_y[0], sub_x[1], sub_y[1]) + + # Get subset box in lat/lon from subset box in y/x + geo_box = coord.box_pixel2geo(pix_box) return pix_box, geo_box