From 082e8ea1f0c87d6b9c76fc98c91337315a117c4b Mon Sep 17 00:00:00 2001 From: Scott Staniewicz Date: Fri, 18 Aug 2023 19:44:41 -0700 Subject: [PATCH 1/9] utils0: add a general epsg-to-epsg reprojection --- src/mintpy/utils/utils0.py | 37 ++++++++++++++++++++++++------------- 1 file changed, 24 insertions(+), 13 deletions(-) diff --git a/src/mintpy/utils/utils0.py b/src/mintpy/utils/utils0.py index b5523358a..06a41a3e7 100644 --- a/src/mintpy/utils/utils0.py +++ b/src/mintpy/utils/utils0.py @@ -275,6 +275,24 @@ 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, use 4326. + + Parameters: x/y - scalar or 1/2D np.ndarray, coordinates in x and y direction + from_epsg - int, EPSG code of `x/y` + Returns: 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 to_latlon(infile, x, y): """Convert x, y in the projection coordinates of the file to lat/lon in degree. @@ -282,27 +300,20 @@ def to_latlon(infile, x, y): https://github.com/Turbo87/utm#utm-to-latitudelongitude Parameters: infile - str, GDAL supported file path - x/y - scalar or 1/2D np.ndarray, coordiantes in x and y direction - Returns: y/x - scalar or 1/2D np.ndarray, coordinates in latitutde and longitude + 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 latitude 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 coordiantes 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, epsg, 4326) + return lat, lon def utm2latlon(meta, easting, northing): From 70b515a90e2c66c369fd8ad71ec606251a0ee05f Mon Sep 17 00:00:00 2001 From: Scott Staniewicz Date: Fri, 18 Aug 2023 19:45:51 -0700 Subject: [PATCH 2/9] gps.py: convert the GPS lat/lon to same EPSG as geometry file --- src/mintpy/objects/gps.py | 23 ++++++++++++++++++----- 1 file changed, 18 insertions(+), 5 deletions(-) diff --git a/src/mintpy/objects/gps.py b/src/mintpy/objects/gps.py index 9f6d9f3b5..6c1460241 100644 --- a/src/mintpy/objects/gps.py +++ b/src/mintpy/objects/gps.py @@ -18,7 +18,13 @@ 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' @@ -493,12 +499,19 @@ def get_los_geometry(self, geom_obj, print_msg=False): if isinstance(geom_obj, str): # geometry file atr = readfile.read_attribute(geom_obj) + file_epsg = int(atr["EPSG"]) + if file_epsg != 4326: + # Convert the GPS position to the same projection as `geom_obj` + x, y = ut0.reproject(lon, lat, from_epsg=4326, to_epsg=file_epsg) + else: + x, y = lon, lat + coord = coordinate(atr, lookup_file=geom_obj) - y, x = coord.geo2radar(lat, lon, print_msg=print_msg)[0:2] + row, col = coord.geo2radar(y, x, 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 = max(0, row); row = min(int(atr['LENGTH'])-1, row) + col = max(0, col); col = min(int(atr['WIDTH'])-1, col) + 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] From 55631317955452831255a03b7a536b91e8eee161 Mon Sep 17 00:00:00 2001 From: Scott Staniewicz Date: Sat, 26 Aug 2023 17:58:26 -0700 Subject: [PATCH 3/9] fix positional args for reproject --- src/mintpy/utils/utils0.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/mintpy/utils/utils0.py b/src/mintpy/utils/utils0.py index 06a41a3e7..07beb46f5 100644 --- a/src/mintpy/utils/utils0.py +++ b/src/mintpy/utils/utils0.py @@ -312,7 +312,7 @@ def to_latlon(infile, x, y): if int(epsg) == 4326: return y, x - lon, lat = reproject(x, y, epsg, 4326) + lon, lat = reproject(x, y, from_epsg=epsg, to_epsg=4326) return lat, lon From 2b18f0fe086f9bb8d992f2bb42100164a07637a6 Mon Sep 17 00:00:00 2001 From: Scott Staniewicz Date: Mon, 28 Aug 2023 17:55:23 -0700 Subject: [PATCH 4/9] convert `Path`s to `str` in `read_attribute` --- src/mintpy/utils/readfile.py | 1 + 1 file changed, 1 insertion(+) diff --git a/src/mintpy/utils/readfile.py b/src/mintpy/utils/readfile.py index b739c1a8e..fde11a3df 100644 --- a/src/mintpy/utils/readfile.py +++ b/src/mintpy/utils/readfile.py @@ -963,6 +963,7 @@ def read_attribute(fname, datasetName=None, metafile_ext=None): ... Returns: atr : dict, attributes dictionary """ + fname = str(fname) # convert from possible Path fdir = os.path.dirname(fname) fbase, fext = os.path.splitext(os.path.basename(fname)) fext = fext.lower() From 2d64373fde80716389234ae8f748a9b2fe847ac9 Mon Sep 17 00:00:00 2001 From: Scott Staniewicz Date: Mon, 28 Aug 2023 18:04:52 -0700 Subject: [PATCH 5/9] add `get_image_rowcol` and `read_image_values to clean up reused functionality --- src/mintpy/objects/gps.py | 35 ++++++++++++++++++++++------------- src/mintpy/utils/utils0.py | 29 ++++++++++++++++++++++++++++- 2 files changed, 50 insertions(+), 14 deletions(-) diff --git a/src/mintpy/objects/gps.py b/src/mintpy/objects/gps.py index 6c1460241..0831ca407 100644 --- a/src/mintpy/objects/gps.py +++ b/src/mintpy/objects/gps.py @@ -499,18 +499,7 @@ def get_los_geometry(self, geom_obj, print_msg=False): if isinstance(geom_obj, str): # geometry file atr = readfile.read_attribute(geom_obj) - file_epsg = int(atr["EPSG"]) - if file_epsg != 4326: - # Convert the GPS position to the same projection as `geom_obj` - x, y = ut0.reproject(lon, lat, from_epsg=4326, to_epsg=file_epsg) - else: - x, y = lon, lat - - coord = coordinate(atr, lookup_file=geom_obj) - row, col = coord.geo2radar(y, x, print_msg=print_msg)[0:2] - # check against image boundary - row = max(0, row); row = min(int(atr['LENGTH'])-1, row) - col = max(0, col); col = min(int(atr['WIDTH'])-1, col) + 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] @@ -523,8 +512,28 @@ def get_los_geometry(self, geom_obj, print_msg=False): else: raise ValueError(f'input geom_obj is neight 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 07beb46f5..b646ad6e1 100644 --- a/src/mintpy/utils/utils0.py +++ b/src/mintpy/utils/utils0.py @@ -15,7 +15,7 @@ # Math / Statistics # Recommend import: # from mintpy.utils import utils as ut - +from __future__ import annotations import math import os @@ -293,6 +293,33 @@ def reproject(x, y, *, from_epsg: int, to_epsg: int): 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, will reproject the latitude/longitude + to the same projection as `atr`. + + 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. From 4ca6363207f34fa776014a2c1c8ab33a9ae1f2b2 Mon Sep 17 00:00:00 2001 From: Scott Staniewicz Date: Mon, 28 Aug 2023 18:05:55 -0700 Subject: [PATCH 6/9] fix error message typo --- src/mintpy/objects/gps.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/mintpy/objects/gps.py b/src/mintpy/objects/gps.py index 0831ca407..9f9c2b06a 100644 --- a/src/mintpy/objects/gps.py +++ b/src/mintpy/objects/gps.py @@ -510,7 +510,7 @@ 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 or np.nan, az_angle or np.nan From c5103184b7061af6e1cada4a4184ebf2bdaf8768 Mon Sep 17 00:00:00 2001 From: Scott Staniewicz Date: Mon, 28 Aug 2023 18:07:18 -0700 Subject: [PATCH 7/9] other typo, unused import --- src/mintpy/objects/gps.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/mintpy/objects/gps.py b/src/mintpy/objects/gps.py index 9f9c2b06a..1c0d4b8d1 100644 --- a/src/mintpy/objects/gps.py +++ b/src/mintpy/objects/gps.py @@ -17,7 +17,6 @@ import numpy as np from pyproj import Geod -from mintpy.objects.coord import coordinate from mintpy.utils import ( ptime, readfile, @@ -463,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 From 4874bda9285b6507c9132221595c992dc3a2b284 Mon Sep 17 00:00:00 2001 From: Zhang Yunjun Date: Wed, 30 Aug 2023 17:42:02 +0800 Subject: [PATCH 8/9] utils0.to_latlon(): update comment & its indentation --- src/mintpy/utils/utils0.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/src/mintpy/utils/utils0.py b/src/mintpy/utils/utils0.py index 4ce8ff977..d18a8f8f7 100644 --- a/src/mintpy/utils/utils0.py +++ b/src/mintpy/utils/utils0.py @@ -15,6 +15,7 @@ # Math / Statistics # Recommend import: # from mintpy.utils import utils as ut + from __future__ import annotations import math @@ -326,15 +327,16 @@ def to_latlon(infile, x, y): 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 # read projection info using gdal ds = gdal.Open(infile) epsg = ds.GetSpatialRef().GetAuthorityCode(None) + # if input file is already in lat/lon, do nothing and return if int(epsg) == 4326: return y, x From 045ea8786ca09c06fe79edcb89ccf97cc680b5ff Mon Sep 17 00:00:00 2001 From: Zhang Yunjun Date: Wed, 30 Aug 2023 18:02:05 +0800 Subject: [PATCH 9/9] utils0.reproject/get_image_rowcol(): update comment --- src/mintpy/utils/utils0.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/src/mintpy/utils/utils0.py b/src/mintpy/utils/utils0.py index d18a8f8f7..db640663b 100644 --- a/src/mintpy/utils/utils0.py +++ b/src/mintpy/utils/utils0.py @@ -279,11 +279,12 @@ def utm_zone2epsg_code(utm_zone): 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, use 4326. + 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` - Returns: x/y - scalar or 1/2D np.ndarray, coordinates in new projection + 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 @@ -297,8 +298,8 @@ def reproject(x, y, *, from_epsg: int, to_epsg: int): 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, will reproject the latitude/longitude - to the same projection as `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