Allow GPS class to work with other image projections

src/mintpy/objects/gps.py

@@ -17,8 +17,13 @@
 import numpy as np
 from pyproj import Geod
 
-from mintpy.objects.coord import coordinate
-from mintpy.utils import ptime, readfile, time_func, utils1 as ut
+from mintpy.utils import (
+    ptime,
+    readfile,
+    time_func,
+    utils0 as ut0,
+    utils1 as ut,
+)
 
 UNR_SITE_LIST_FILE = 'http://geodesy.unr.edu/NGLStationPages/DataHoldings.txt'
 
@@ -457,7 +462,7 @@ def displacement_enu2los(self, inc_angle:float, az_angle:float, gps_comp='enu2lo
         """Convert displacement in ENU to LOS direction.
 
         Parameters: inc_angle     - float, LOS incidence angle in degree
-                    az_angle      - float, LOS aziuth angle in degree
+                    az_angle      - float, LOS azimuth angle in degree
                                     from the north, defined as positive in clock-wise direction
                     gps_comp      - str, GPS components used to convert to LOS direction
                     horz_az_angle - float, fault azimuth angle used to convert horizontal to fault-parallel
@@ -493,12 +498,8 @@ def get_los_geometry(self, geom_obj, print_msg=False):
         if isinstance(geom_obj, str):
             # geometry file
             atr = readfile.read_attribute(geom_obj)
-            coord = coordinate(atr, lookup_file=geom_obj)
-            y, x = coord.geo2radar(lat, lon, print_msg=print_msg)[0:2]
-            # check against image boundary
-            y = max(0, y);  y = min(int(atr['LENGTH'])-1, y)
-            x = max(0, x);  x = min(int(atr['WIDTH'])-1, x)
-            box = (x, y, x+1, y+1)
+            row, col = ut0.get_image_rowcol(atr, lat=lat, lon=lon)
+            box = (col, row, col+1, row+1)
             inc_angle = readfile.read(geom_obj, datasetName='incidenceAngle', box=box, print_msg=print_msg)[0][0,0]
             az_angle  = readfile.read(geom_obj, datasetName='azimuthAngle',   box=box, print_msg=print_msg)[0][0,0]
 
@@ -508,10 +509,30 @@ def get_los_geometry(self, geom_obj, print_msg=False):
             az_angle  = ut.heading2azimuth_angle(float(geom_obj['HEADING']))
 
         else:
-            raise ValueError(f'input geom_obj is neight str nor dict: {geom_obj}')
+            raise ValueError(f'input geom_obj is neither str nor dict: {geom_obj}')
 
-        return inc_angle, az_angle
+        return inc_angle or np.nan, az_angle or np.nan
 
+    def get_image_values(self, filename, datasetName=None, pad: int = 0, print_msg=False):
+        """Read the value from `filename` at the pixel nearest to the GPS station.
+
+        Parameters: atr         - dict, mintpy attributes that includes "EPSG"
+                    datasetName - str, If `filename` is an HDF5 file, dataset to read from
+                    pad         - int, default = 0. Number of pixels of padding to read around
+                                  the nearest pixel.
+                                  `pad=0` only reads the closes pixel.
+                                  `pad=1` will return a 3x3 ndarray of image values.
+        Returns:    scalar, (or ndarray if `pad > 1`)
+                                  The values closes to the GPS station in `filename`.
+        """
+        lat, lon = self.get_stat_lat_lon(print_msg=print_msg)
+
+        atr = readfile.read_attribute(filename)
+        row, col = ut0.get_image_rowcol(atr, lat=lat, lon=lon)
+
+        box = (col - pad, row - pad, col + pad + 1, row + pad + 1)
+        values = readfile.read(filename, datasetName=datasetName, box=box, print_msg=print_msg)[0]
+        return np.squeeze(values)
 
     def read_gps_los_displacement(self, geom_obj, start_date=None, end_date=None, ref_site=None,
                                   gps_comp='enu2los', horz_az_angle=-90., print_msg=False):

src/mintpy/utils/utils0.py

@@ -16,6 +16,7 @@
 # Recommend import:
 #   from mintpy.utils import utils as ut
 
+from __future__ import annotations
 
 import math
 import os
@@ -275,34 +276,74 @@ def utm_zone2epsg_code(utm_zone):
     return epsg_code
 
 
+def reproject(x, y, *, from_epsg: int, to_epsg: int):
+    """Convert x, y in the projection EPSG:`from_epsg` to EPSG:`to_epsg`.
+
+    For lon/lat, the EPSG code is 4326.
+
+    Parameters: x/y       - scalar or 1/2D np.ndarray, coordinates in x and y direction
+                from_epsg - int, EPSG code of `x/y`
+                to_epsg   - int, EPSG code of `new_x/y`
+    Returns:    new_x/y   - scalar or 1/2D np.ndarray, coordinates in new projection
+    """
+    from pyproj import CRS, Transformer
+
+    transformer = Transformer.from_crs(
+        CRS.from_epsg(from_epsg), CRS.from_epsg(to_epsg), always_xy=True
+    )
+    new_x, new_y = transformer.transform(x, y)
+    return new_x, new_y
+
+
+def get_image_rowcol(atr: dict, lat: float, lon: float, print_msg: bool = False):
+    """Get the (row, column) of `lat`,`lon` for an image with attributes `atr`.
+
+    For images not using EPSG:4326 (as in lat/lon), will reproject the latitude/longitude
+    to the same projection as `atr['EPSG']`.
+
+    Parameters: atr       - dict, mintpy attributes that includes "EPSG"
+                lat/lon   - float, latitude/longitude of point of interest
+                print_msg - bool, enable verbose printing
+    Returns:    row, col  - integers for the closest row/column in `atr`.
+    """
+    from mintpy.objects.coord import coordinate
+    file_epsg = int(atr["EPSG"])
+    if file_epsg != 4326:
+        # Convert the GPS position to the same projection as `geom_obj`
+        x, y = reproject(lon, lat, from_epsg=4326, to_epsg=file_epsg)
+    else:
+        x, y = lon, lat
+
+    coord = coordinate(atr)
+    row, col = coord.geo2radar(y, x, print_msg=print_msg)[0:2]
+
+    if row < 0 or col < 0 or row > int(atr['LENGTH'])-1 or col > int(atr['WIDTH'])-1:
+        raise ValueError(f"{lat = }, {lon = } is outside the image bounds")
+    return row, col
+
+
 def to_latlon(infile, x, y):
     """Convert x, y in the projection coordinates of the file to lat/lon in degree.
 
     Similar functionality also exists in utm.to_latlon() at:
         https://github.com/Turbo87/utm#utm-to-latitudelongitude
 
-    Parameters: infile - str, GDAL supported file path
-                x/y    - scalar or 1/2D np.ndarray, coordinates in x and y direction
-    Returns:    y/x    - scalar or 1/2D np.ndarray, coordinates in latitutde and longitude
+    Parameters: infile  - str, GDAL supported file path
+                x/y     - scalar or 1/2D np.ndarray, coordinates in x and y direction
+    Returns:    lat/lon - scalar or 1/2D np.ndarray, coordinates in latitutde and longitude
     """
     from osgeo import gdal
-    from pyproj import Proj, Transformer
 
     # read projection info using gdal
     ds = gdal.Open(infile)
-    srs = ds.GetSpatialRef()
+    epsg = ds.GetSpatialRef().GetAuthorityCode(None)
 
     # if input file is already in lat/lon, do nothing and return
-    if (not srs.IsProjected()) and (srs.GetAttrValue('unit') == 'degree'):
+    if int(epsg) == 4326:
         return y, x
 
-    # convert coordinates using pyproj
-    # note that Transform.from_proj(x, y, always_xy=True) convert the x, y to lon, lat
-    p_in = Proj(ds.GetProjection())
-    p_out = Proj('epsg:4326')
-    transformer = Transformer.from_proj(p_in, p_out)
-    y, x = transformer.transform(x, y)
-    return y, x
+    lon, lat = reproject(x, y, from_epsg=epsg, to_epsg=4326)
+    return lat, lon
 
 
 def utm2latlon(meta, easting, northing):