diff --git a/src/mintpy/objects/gps.py b/src/mintpy/objects/gps.py index 76e97d413..4658a95ba 100644 --- a/src/mintpy/objects/gps.py +++ b/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): diff --git a/src/mintpy/utils/utils0.py b/src/mintpy/utils/utils0.py index d00b334fe..db640663b 100644 --- a/src/mintpy/utils/utils0.py +++ b/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):