Skip to content

Commit

Permalink
bugfix for coord.lalo2yx() when lat/lon have diff sizes (insarlab#1160
Browse files Browse the repository at this point in the history
)

+ `objects.coord.coordinate`:
   - `_clean_coord()`: 1) take two arg instead of one, and ensure the two coords have the same size; 2) support `None` as input
   - `lalo2yx/yx2lalo()`: support `None` as inputs

+ `subset.subset_input_dict2box()`: use `None` as input to convert one type of coord only.

+ `utils.arg_utils.add_subset_argument()`: remove the unused option name alias (`--suby/x/lat/lon`), as they are not very readable and not frequently used.
  • Loading branch information
yunjunz authored and ehavazli committed Mar 19, 2024
1 parent fe6f84b commit 6d27181
Show file tree
Hide file tree
Showing 5 changed files with 61 additions and 33 deletions.
2 changes: 1 addition & 1 deletion docs/api/attributes.md
Original file line number Diff line number Diff line change
Expand Up @@ -52,7 +52,7 @@ The following attributes vary for each interferogram:
+ MODIFICATION_TIME = dataset modification time, exists in ifgramStack.h5 file for 3D dataset, used for "--update" option of unwrap error corrections.
+ NCORRLOOKS = number of independent looks, as explained in [SNAPHU](https://web.stanford.edu/group/radar/softwareandlinks/sw/snaphu/snaphu.conf.full)
+ UTM_ZONE = [UTM zone](https://docs.up42.com/data/reference/utm#utm-wgs84), comprises a zone number and a hemisphere, e.g. 11N, 60S, for geocoded file with UTM projection only.
+ EPSG = EPSG code for coordinate systems, for geocoded files only. Check [here] for its relationship with UTM zone.
+ EPSG = EPSG code for coordinate systems, for geocoded files only. Check [here](https://docs.up42.com/data/reference/utm#utm-wgs84) for its relationship with UTM zone.
+ CENTER_INCIDENCE_ANGLE = incidence angle in degrees at the scene center, read from the 2D incidence angle matrix, for isce2 files only.

### Reference ###
Expand Down
69 changes: 48 additions & 21 deletions src/mintpy/objects/coord.py
Original file line number Diff line number Diff line change
Expand Up @@ -76,18 +76,39 @@ def open(self):
self.src_metadata.pop(key)


def _clean_coord(self, coord_in):
# note: np.float128 is not supported on Windows OS, use np.longdouble as a platform neutral syntax
float_types = (float, np.float16, np.float32, np.float64, np.longdouble)
int_types = (int, np.int16, np.int32, np.int64)

if isinstance(coord_in, np.ndarray):
coord_in = coord_in.tolist()
if isinstance(coord_in, int_types + float_types):
coord_in = [coord_in]
coord_in = list(coord_in)
def _clean_coord(self, coord1, coord2):
"""Ensure input coordinate as list type and same size."""

def _to_list(x):
# convert numpy array to list or scalar
if isinstance(x, np.ndarray):
x = x.tolist()
# convert scalar / None to list
# Note: np.float128 is not supported on Windows OS,
# use np.longdouble as a platform neutral syntax
float_types = (float, np.float16, np.float32, np.float64, np.longdouble)
int_types = (int, np.int16, np.int32, np.int64)
if isinstance(x, int_types + float_types):
x = [x]
elif x is None:
x = [None]
x = list(x)
return x

# convert to list type
coord1 = _to_list(coord1)
coord2 = _to_list(coord2)

# ensure the same size
if len(coord1) != len(coord2):
if len(coord1) == 1 and len(coord2) > 1:
coord1 *= len(coord2)
elif len(coord1) > 1 and len(coord2) == 1:
coord2 *= len(coord1)
else:
raise ValueError('Input two coordinates do NOT have the same size!')

return coord_in
return coord1, coord2


def lalo2yx(self, lat_in, lon_in):
Expand All @@ -97,17 +118,18 @@ def lalo2yx(self, lat_in, lon_in):
lon_in - list / tuple / 1D np.ndarray / float, coordinate(s) in longitude / easting
Returns: coord_out - tuple(list / float / 1D np.ndarray), coordinates(s) in (row, col) numbers
Example: 300, 1000 = coordinate.lalo2yx(32.1, 130.5)
300 = coordinate.lalo2yx(32.1, None)[0]
([300, 301], [1000, 1001]) = coordinate.lalo2yx((32.1, 32.101), (130.5, 130.501))
"""
self.open()
if not self.geocoded:
raise ValueError('Input file is NOT geocoded.')

lat_in = self._clean_coord(lat_in)
lon_in = self._clean_coord(lon_in)
lat_in, lon_in = self._clean_coord(lat_in, lon_in)

# attempts to convert lat/lon to utm coordinates if needed.
if ('UTM_ZONE' in self.src_metadata
if (lat_in is not None and lon_in is not None
and 'UTM_ZONE' in self.src_metadata
and np.max(np.abs(lat_in)) <= 90
and np.max(np.abs(lon_in)) <= 360):
lat_in, lon_in = ut0.latlon2utm(np.array(lat_in), np.array(lon_in))
Expand All @@ -117,8 +139,10 @@ def lalo2yx(self, lat_in, lon_in):
x_out = []
for lat_i, lon_i in zip(lat_in, lon_in):
# plus 0.01 to be more robust in practice
y_out.append(int(np.floor((lat_i - self.lat0) / self.lat_step + 0.01)))
x_out.append(int(np.floor((lon_i - self.lon0) / self.lon_step + 0.01)))
y_i = None if lat_i is None else int(np.floor((lat_i - self.lat0) / self.lat_step + 0.01))
x_i = None if lon_i is None else int(np.floor((lon_i - self.lon0) / self.lon_step + 0.01))
y_out.append(y_i)
x_out.append(x_i)

# output format
if len(y_out) == 1 and len(x_out) == 1:
Expand All @@ -137,22 +161,25 @@ def yx2lalo(self, y_in, x_in):
x_in - list / tuple / 1D np.ndarray / int, coordinate(s) in col number
Returns: coord_out - tuple(list / 1D np.ndarray / int), coordinate(s) in (lat/northing, lon/easting)
Example: 32.1, 130.5 = coordinate.yx2lalo(300, 1000)
32.1 = coordinate.yx2lalo(300, None)[0]
([32.1, 32.101], [130.5, 130.501]) = coordinate.lalo2yx([(300, 301), (1000, 1001)])
"""
self.open()
if not self.geocoded:
raise ValueError('Input file is NOT geocoded.')

# Convert to List if input is String
y_in = self._clean_coord(y_in)
x_in = self._clean_coord(x_in)
y_in, x_in = self._clean_coord(y_in, x_in)

# convert coordinates
lat_out = []
lon_out = []
for y_i, x_i in zip(y_in, x_in):
lat_out.append((y_i + 0.5) * self.lat_step + self.lat0)
lon_out.append((x_i + 0.5) * self.lon_step + self.lon0)
lat_i = None if y_i is None else (y_i + 0.5) * self.lat_step + self.lat0
lon_i = None if x_i is None else (x_i + 0.5) * self.lon_step + self.lon0
lat_out.append(lat_i)
lon_out.append(lon_i)

# output format
if len(lat_out) == 1 and len(lon_out) == 1:
coord_out = tuple([lat_out[0], lon_out[0]])
else:
Expand Down
4 changes: 2 additions & 2 deletions src/mintpy/subset.py
Original file line number Diff line number Diff line change
Expand Up @@ -184,14 +184,14 @@ 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'], 0)[0]
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(0, subset_dict['subset_lon'])[1]
sub_x = coord.lalo2yx(None, subset_dict['subset_lon'])[1]
elif subset_dict['subset_x']:
sub_x = subset_dict['subset_x']
else:
Expand Down
8 changes: 4 additions & 4 deletions src/mintpy/utils/arg_utils.py
Original file line number Diff line number Diff line change
Expand Up @@ -447,14 +447,14 @@ def add_save_argument(parser):
def add_subset_argument(parser, geo=True):
"""Argument group parser for subset options"""
sub = parser.add_argument_group('Subset', 'Display dataset in subset range')
sub.add_argument('--sub-x','--subx','--subset-x', dest='subset_x', type=int, nargs=2,
sub.add_argument('--sub-x','--subset-x', dest='subset_x', type=int, nargs=2,
metavar=('XMIN', 'XMAX'), help='subset display in x/cross-track/range direction')
sub.add_argument('--sub-y','--suby','--subset-y', dest='subset_y', type=int, nargs=2,
sub.add_argument('--sub-y','--subset-y', dest='subset_y', type=int, nargs=2,
metavar=('YMIN', 'YMAX'), help='subset display in y/along-track/azimuth direction')
if geo:
sub.add_argument('--sub-lat','--sublat','--subset-lat', dest='subset_lat', type=float, nargs=2,
sub.add_argument('--sub-lat','--subset-lat', dest='subset_lat', type=float, nargs=2,
metavar=('LATMIN', 'LATMAX'), help='subset display in latitude')
sub.add_argument('--sub-lon','--sublon','--subset-lon', dest='subset_lon', type=float, nargs=2,
sub.add_argument('--sub-lon','--subset-lon', dest='subset_lon', type=float, nargs=2,
metavar=('LONMIN', 'LONMAX'), help='subset display in longitude')
return parser

Expand Down
11 changes: 6 additions & 5 deletions src/mintpy/utils/utils0.py
Original file line number Diff line number Diff line change
Expand Up @@ -270,7 +270,9 @@ def touch(fname_list, times=None):
def utm_zone2epsg_code(utm_zone):
"""Convert UTM Zone string to EPSG code.
Reference: https://docs.up42.com/data/reference/utm#utm-wgs84
Reference:
https://docs.up42.com/data/reference/utm#utm-wgs84
https://pyproj4.github.io/pyproj/stable/examples.html#initializing-crs
Parameters: utm_zone - str, atr['UTM_ZONE'], comprises a zone number
and a hemisphere, e.g. 11N, 36S, etc.
Expand All @@ -290,7 +292,9 @@ def utm_zone2epsg_code(utm_zone):
def epsg_code2utm_zone(epsg_code):
"""Convert EPSG code to UTM Zone string.
Reference: https://docs.up42.com/data/reference/utm#utm-wgs84
Reference:
https://docs.up42.com/data/reference/utm#utm-wgs84
https://pyproj4.github.io/pyproj/stable/examples.html#initializing-crs
Parameters: epsg_code - str / int, EPSG code
Returns: utm_zone - str, atr['UTM_ZONE'], comprises a zone number
Expand All @@ -299,9 +303,6 @@ def epsg_code2utm_zone(epsg_code):
Examples: utm_zone = epsg_code2utm_zone('32736')
"""
from pyproj import CRS

# link: https://pyproj4.github.io/pyproj/stable/examples.html#initializing-crs
# alternatively, use CRS.from_user_input(epsg_code) or CRS.from_string(f'EPSG:{epsg_code}')
crs = CRS.from_epsg(epsg_code)
utm_zone = crs.utm_zone
if not utm_zone:
Expand Down

0 comments on commit 6d27181

Please sign in to comment.