UTM support for view --show-gps + dual-coord input for `lalo2yx/yx2…

…lalo()` (#1153)

+ utils.plot.py:
   - set `plt.rcParams["axes.formatter.limits"]` to disable scientific notation for UTM coordinates.
   - `plot_gps()`: convert between lat/lon and UTM coordinates to support files in UTM coordinates for `view.py --show-gps` option

+ objects.coord.py:
   - `yx2lalo()` and `lalo2yx()`: use dual-coord inputs, instead of single-coord inputs, for more convenient usage, and update all the instances throughout the repo
   - `lalo2yx()` and `geo2radar()`: attempt to convert lat/lon into UTM coordinates, if file/metadata is in UTM while the inputs are in WGS 84
   - add `_clean_coord()` for type checking / cleaning and easy re-use

---------

Co-authored-by: Zhang Yunjun <[email protected]>

src/mintpy/asc_desc2horz_vert.py

@@ -186,8 +186,7 @@ def run_asc_desc2horz_vert(inps):
     for i, (atr, fname) in enumerate(zip(atr_list, inps.file)):
         # overlap SNWE --> box to read for each specific file
         coord = ut.coordinate(atr)
-        x0 = coord.lalo2yx(W, coord_type='lon')
-        y0 = coord.lalo2yx(N, coord_type='lat')
+        y0, x0 = coord.lalo2yx(N, W)
         box = (x0, y0, x0 + width, y0 + length)
 
         # read data

src/mintpy/cli/plot_transection.py

@@ -126,8 +126,8 @@ def cmd_line_parse(iargs=None):
 
     # check: --start/end-lalo (start/end_lalo --> start/end_yx)
     if inps.start_lalo and inps.end_lalo:
-        [y0, y1] = inps.coord.lalo2yx([inps.start_lalo[0], inps.end_lalo[0]], coord_type='lat')
-        [x0, x1] = inps.coord.lalo2yx([inps.start_lalo[1], inps.end_lalo[1]], coord_type='lon')
+        [y0, y1], [x0, x1] = inps.coord.lalo2yx([inps.start_lalo[0], inps.end_lalo[0]],
+                                                [inps.start_lalo[1], inps.end_lalo[1]])
         inps.start_yx = [y0, x0]
         inps.end_yx = [y1, x1]
 

src/mintpy/legacy/transect_legacy2.py

@@ -233,8 +233,7 @@ def transect_yx(z, atr, start_yx, end_yx, interpolation='nearest'):
 def transect_lalo(z, atr, start_lalo, end_lalo, interpolation='nearest'):
     """Extract 2D matrix (z) value along the line [start_lalo, end_lalo]"""
     coord = ut.coordinate(atr)
-    [y0, y1] = coord.lalo2yx([start_lalo[0], end_lalo[0]], coord_type='lat')
-    [x0, x1] = coord.lalo2yx([start_lalo[1], end_lalo[1]], coord_type='lon')
+    [y0, y1], [x0, x1] = coord.lalo2yx([start_lalo[0], end_lalo[0]], [start_lalo[1], end_lalo[1]])
     transect = transect_yx(z, atr, [y0, x0], [y1, x1], interpolation)
     return transect
 
@@ -298,8 +297,8 @@ def plot_transect_location(ax, inps):
 
     coord = ut.coordinate(atr0)
     if inps.start_lalo and inps.end_lalo:
-        [y0, y1] = coord.lalo2yx([inps.start_lalo[0], inps.end_lalo[0]], coord_type='lat')
-        [x0, x1] = coord.lalo2yx([inps.start_lalo[1], inps.end_lalo[1]], coord_type='lon')
+        [y0, y1], [x0, x1] = coord.lalo2yx([inps.start_lalo[0], inps.end_lalo[0]],
+                                           [inps.start_lalo[1], inps.end_lalo[1]])
         inps.start_yx = [y0, x0]
         inps.end_yx = [y1, x1]
 
@@ -316,8 +315,7 @@ def format_coord(x, y):
         if 0 <= col < data0.shape[1] and 0 <= row < data0.shape[0]:
             z = data0[row, col]
             if 'X_FIRST' in atr0.keys():
-                lat = coord.yx2lalo(row, coord_type='row')
-                lon = coord.yx2lalo(col, coord_type='col')
+                lat, lon = coord.yx2lalo(row, col)
                 return f'lon={lon:.4f}, lat={lat:.4f}, x={x:.0f},  y={y:.0f},  value={z:.4f}'
             else:
                 return f'x={x:.0f},  y={y:.0f},  value={z:.4f}'

src/mintpy/objects/coord.py

@@ -54,6 +54,7 @@ def __init__(self, metadata, lookup_file=None):
         self.lut_y = None
         self.lut_x = None
 
+
     def open(self):
         self.earth_radius = float(self.src_metadata.get('EARTH_RADIUS', EARTH_RADIUS))
 
@@ -68,82 +69,95 @@ def open(self):
             if self.lookup_file:
                 self.lut_metadata = readfile.read_attribute(self.lookup_file[0])
 
-    def lalo2yx(self, coord_in, coord_type):
-        """convert geo coordinates into radar coordinates for Geocoded file only
-        Parameters: geoCoord   - list / tuple / 1D np.ndarray / float, coordinate(s) in latitude or longitude
-                    metadata   - dict, dictionary of file attributes
-                    coord_type - str, coordinate type: latitude or longitude
-        Example:    300 = coordinate.lalo2yx(32.104990,    metadata,'lat')
-                    [1000,1500] = coordinate.lalo2yx([130.5,131.4],metadata,'lon')
-        """
-        self.open()
-        if not self.geocoded:
-            raise ValueError('Input file is not geocoded.')
+        # remove empty attributes [to facilitate checking laterward]
+        key_list = ['UTM_ZONE']
+        for key in key_list:
+            if key in self.src_metadata and not self.src_metadata[key]:
+                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)
 
-        # input format
         if isinstance(coord_in, np.ndarray):
             coord_in = coord_in.tolist()
-        # note: np.float128 is not supported on Windows OS, use np.longdouble as a platform neutral syntax
-        if isinstance(coord_in, (float, np.float16, np.float32, np.float64, np.longdouble)):
+        if isinstance(coord_in, int_types + float_types):
             coord_in = [coord_in]
         coord_in = list(coord_in)
 
+        return coord_in
+
+
+    def lalo2yx(self, lat_in, lon_in):
+        """convert geo coordinates into radar coordinates for Geocoded file only.
+
+        Parameters: lat_in    - list / tuple / 1D np.ndarray / float, coordinate(s) in latitude / northing
+                    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, 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)
+
+        # attempts to convert lat/lon to utm coordinates if needed.
+        if ('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))
+
         # convert coordinates
-        coord_type = coord_type.lower()
-        coord_out = []
-        for coord_i in coord_in:
+        y_out = []
+        x_out = []
+        for lat_i, lon_i in zip(lat_in, lon_in):
             # plus 0.01 to be more robust in practice
-            if coord_type.startswith('lat'):
-                coord_o = int(np.floor((coord_i - self.lat0) / self.lat_step + 0.01))
-            elif coord_type.startswith('lon'):
-                coord_o = int(np.floor((coord_i - self.lon0) / self.lon_step + 0.01))
-            else:
-                raise ValueError('Unrecognized coordinate type: '+coord_type)
-            coord_out.append(coord_o)
+            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)))
 
         # output format
-        if len(coord_out) == 1:
-            coord_out = coord_out[0]
-        elif isinstance(coord_in, tuple):
-            coord_out = tuple(coord_out)
+        if len(y_out) == 1 and len(x_out) == 1:
+            coord_out = tuple([y_out[0], x_out[0]])
+        else:
+            coord_out = tuple([y_out, x_out])
+
         return coord_out
 
 
-    def yx2lalo(self, coord_in, coord_type):
+    def yx2lalo(self, y_in, x_in):
         """convert radar coordinates into geo coordinates (pixel center)
-            for Geocoded file only
-        Parameters: coord_in   _ list / tuple / 1D np.ndarray / int, coordinate(s) in row or col in int
-                    metadata   _ dict, dictionary of file attributes
-                    coord_type _ str, coordinate type: row / col / y / x
-        Example:    32.104990 = coord_yx2lalo(300, metadata, 'y')
-                    [130.5,131.4] = coord_yx2lalo([1000,1500], metadata, 'x')
+            for Geocoded file only.
+
+        Parameters: y_in      - list / tuple / 1D np.ndarray / int, coordinate(s) in row number
+                    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, 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.')
+            raise ValueError('Input file is NOT geocoded.')
 
         # Convert to List if input is String
-        if isinstance(coord_in, np.ndarray):
-            coord_in = coord_in.tolist()
-        if isinstance(coord_in, (int, np.int16, np.int32, np.int64)):
-            coord_in = [coord_in]
-        coord_in = list(coord_in)
+        y_in = self._clean_coord(y_in)
+        x_in = self._clean_coord(x_in)
 
-        coord_type = coord_type.lower()
-        coord_out = []
-        for i in range(len(coord_in)):
-            if coord_type.startswith(('row', 'y', 'az', 'azimuth')):
-                coord = (coord_in[i] + 0.5) * self.lat_step + self.lat0
-            elif coord_type.startswith(('col', 'x', 'rg', 'range')):
-                coord = (coord_in[i] + 0.5) * self.lon_step + self.lon0
-            else:
-                raise ValueError('Unrecognized coordinate type: '+coord_type)
-            coord_out.append(coord)
+        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)
+
+        if len(lat_out) == 1 and len(lon_out) == 1:
+            coord_out = tuple([lat_out[0], lon_out[0]])
+        else:
+            coord_out = tuple([lat_out, lon_out])
 
-        if len(coord_out) == 1:
-            coord_out = coord_out[0]
-        elif isinstance(coord_in, tuple):
-            coord_out = tuple(coord_out)
         return coord_out
 
 
@@ -197,6 +211,7 @@ def read_lookup_table(self, print_msg=True):
         self.lut_x = readfile.read(self.lookup_file[1], datasetName=ds_name_x, print_msg=print_msg)[0]
         return self.lut_y, self.lut_x
 
+
     def _read_geo_lut_metadata(self):
         """Read lat/lon0, lat/lon_step_deg, lat/lon_step into a Namespace - lut"""
         lut = Namespace()
@@ -212,14 +227,16 @@ def _read_geo_lut_metadata(self):
         lut.lon_step = lut.lon_step_deg * np.pi/180.0 * self.earth_radius * np.cos(lat_c * np.pi/180)
         return lut
 
+
     def geo2radar(self, lat, lon, print_msg=True, debug_mode=False):
         """Convert geo coordinates into radar coordinates.
+
         Parameters: lat/lon   - np.array / float, latitude/longitude
         Returns:    az/rg     - np.array / int, range/azimuth pixel number
                     az/rg_res - float, residul/uncertainty of coordinate conversion
         """
-        # ensure longitude convention to be consistent with src_metadata
-        if 'X_FIRST' in self.src_metadata.keys():
+        # check 1: ensure longitude convention to be consistent with src_metadata [for WGS84 coordinates]
+        if 'X_FIRST' in self.src_metadata.keys() and 'UTM_ZONE' not in self.src_metadata.keys():
             width = int(self.src_metadata['WIDTH'])
             lon_step = float(self.src_metadata['X_STEP'])
             min_lon = float(self.src_metadata['X_FIRST'])
@@ -244,10 +261,15 @@ def geo2radar(self, lat, lon, print_msg=True, debug_mode=False):
                 else:
                     lon[lon > 180.] -= 360
 
+        # check 2: attempts to convert lat/lon to utm coordinates if needed
+        if ('UTM_ZONE' in self.src_metadata
+                and np.max(np.abs(lat)) <= 90
+                and np.max(np.abs(lon)) <= 360):
+            lat, lon = ut0.latlon2utm(np.array(lat), np.array(lon))
+
         self.open()
         if self.geocoded:
-            az = self.lalo2yx(lat, coord_type='lat')
-            rg = self.lalo2yx(lon, coord_type='lon')
+            az, rg = self.lalo2yx(lat, lon)
             return az, rg, 0, 0
 
         if not isinstance(lat, np.ndarray):
@@ -336,8 +358,7 @@ def radar2geo(self, az, rg, print_msg=True, debug_mode=False):
         """
         self.open()
         if self.geocoded:
-            lat = self.yx2lalo(az, coord_type='az')
-            lon = self.yx2lalo(rg, coord_type='rg')
+            lat, lon = self.yx2lalo(az, rg)
             return lat, lon, 0, 0
 
         if not isinstance(az, np.ndarray):
@@ -398,14 +419,14 @@ def radar2geo(self, az, rg, print_msg=True, debug_mode=False):
             lon_resid = 0.
         return lat, lon, lat_resid, lon_resid
 
+
     def box_pixel2geo(self, pixel_box):
         """Convert pixel_box to geo_box in UL corner
         Parameters: pixel_box - list/tuple of 4 int   in (x0, y0, x1, y1)
         Returns:    geo_box   - tuple      of 4 float in (W, N, E, S)
         """
         try:
-            lat = self.yx2lalo([pixel_box[1], pixel_box[3]], coord_type='y')
-            lon = self.yx2lalo([pixel_box[0], pixel_box[2]], coord_type='x')
+            lat, lon = self.yx2lalo([pixel_box[1], pixel_box[3]], [pixel_box[0], pixel_box[2]])
             # shift lat from pixel center to the UL corner
             lat = [i - self.lat_step / 2.0 for i in lat]
             lon = [i - self.lon_step / 2.0 for i in lon]
@@ -414,19 +435,20 @@ def box_pixel2geo(self, pixel_box):
             geo_box = None
         return geo_box
 
+
     def box_geo2pixel(self, geo_box):
         """Convert geo_box to pixel_box
         Parameters: geo_box   - tuple      of 4 float in (W, N, E, S)
         Returns:    pixel_box - list/tuple of 4 int   in (x0, y0, x1, y1)
         """
         try:
-            y = self.lalo2yx([geo_box[1], geo_box[3]], coord_type='latitude')
-            x = self.lalo2yx([geo_box[0], geo_box[2]], coord_type='longitude')
+            y, x = self.lalo2yx([geo_box[1], geo_box[3]], [geo_box[0], geo_box[2]])
             pixel_box = (x[0], y[0], x[1], y[1])
         except:
             pixel_box = None
         return pixel_box
 
+
     def bbox_radar2geo(self, pix_box, print_msg=False):
         """Calculate bounding box in lat/lon for file in geo coord, based on input radar/pixel box
         Parameters: pix_box - tuple of 4 int, in (x0, y0, x1, y1)
@@ -440,6 +462,7 @@ def bbox_radar2geo(self, pix_box, print_msg=False):
                    np.max(lon) + buf, np.min(lat) - buf)
         return geo_box
 
+
     def bbox_geo2radar(self, geo_box, 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
@@ -454,6 +477,7 @@ def bbox_geo2radar(self, geo_box, print_msg=False):
                    np.max(x) + buf, np.max(y) + buf)
         return pix_box
 
+
     def check_box_within_data_coverage(self, pixel_box, print_msg=True):
         """Check the subset box's conflict with data coverage
         Parameters:  pixel_box - 4-tuple of int, indicating y/x coordinates of subset

src/mintpy/plot_coherence_matrix.py

@@ -211,8 +211,7 @@ def format_coord(x, y):
     def update_coherence_matrix(self, event):
         if event.inaxes == self.ax_img:
             if self.fig_coord == 'geo':
-                yx = [self.coord.lalo2yx(event.ydata, coord_type='lat'),
-                      self.coord.lalo2yx(event.xdata, coord_type='lon')]
+                yx = self.coord.lalo2yx(event.ydata, event.xdata)
             else:
                 yx = [int(event.ydata+0.5),
                       int(event.xdata+0.5)]

src/mintpy/plot_transection.py

@@ -168,8 +168,7 @@ def draw_line(self, start_yx, end_yx):
 
         # convert coordinates accordingly
         if 'Y_FIRST' in self.atr.keys():
-            ys = self.coord.yx2lalo([start_yx[0], end_yx[0]], coord_type='y')
-            xs = self.coord.yx2lalo([start_yx[1], end_yx[1]], coord_type='x')
+            ys, xs = self.coord.yx2lalo([start_yx[0], end_yx[0]], [start_yx[1], end_yx[1]])
         else:
             ys = [start_yx[0], end_yx[0]]
             xs = [start_yx[1], end_yx[1]]

src/mintpy/reference_point.py

@@ -171,8 +171,7 @@ def reference_point_attribute(atr, y, x):
     atrNew['REF_X'] = str(x)
     coord = ut.coordinate(atr)
     if 'X_FIRST' in atr.keys():
-        atrNew['REF_LAT'] = str(coord.yx2lalo(y, coord_type='y'))
-        atrNew['REF_LON'] = str(coord.yx2lalo(x, coord_type='x'))
+        atrNew['REF_LAT'], atrNew['REF_LON'] = (str(val) for val in coord.yx2lalo(y, x))
     return atrNew
 
 

src/mintpy/subset.py

@@ -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'], coord_type='latitude')
+        sub_y = coord.lalo2yx(subset_dict['subset_lat'], 0)[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(subset_dict['subset_lon'], coord_type='longitude')
+        sub_x = coord.lalo2yx(0, subset_dict['subset_lon'])[1]
     elif subset_dict['subset_x']:
         sub_x = subset_dict['subset_x']
     else: