Skip to content

Commit

Permalink
support LAT/LON_REF1/2/3/4 in UTM coordinates
Browse files Browse the repository at this point in the history
+ docs/api/attributes.md: ensure all LAT/LON attributes represents both latitude/northing and longitude/easting in degrees/meters

+ utils.utils0.utm2latlon/latlon2utm(): support list/tuple input types

+ prep_hyp3:
   - do not convert UTM to lat/lon for LAT/LON_REF1/2/3/4
   - add HDF-EOS5 metadata (`beam_mode` and `unwrap_method`) from the HyP3 metadata file

+ save_hdfeos5.metadata_mintpy2unavco(): convert UTM to lat/lon while calculating the scene_footprint

+ tests/configs/RidgecrestSenDT71: turn ON save.hdfEos5 option

+ tropo_pyaps3.get_bounding_box(): convert UTM to lat/lon
  • Loading branch information
yunjunz committed Mar 18, 2024
1 parent 0c45019 commit 51fe70e
Show file tree
Hide file tree
Showing 6 changed files with 58 additions and 33 deletions.
7 changes: 4 additions & 3 deletions docs/api/attributes.md
Original file line number Diff line number Diff line change
Expand Up @@ -9,7 +9,7 @@ If using ROI_PAC as the InSAR processor, both **baseline parameter RSC** file (i
+ X/Y_FIRST = (for geocoded product) Longitude/easting/X and latitude/northing/Y coordinate in degrees/meters of the upper left corner of the first pixel.
+ X/Y_STEP = (for geocoded product) Ground resolution in degrees/meters in X/Y direction.
+ X/Y_UNIT = (for geocoded product) Coordinate unit in X/Y direction: degrees or meters.
+ LAT/LON_REF1/2/3/4 = Latitude/longitude at corner 1/2/3/4 (degree), used in save_unavco, PyAPS (DEM file in radar coord), not accurate; number named in order of first line near/far range, last line near/far range.
+ LAT/LON_REF1/2/3/4 = Latitude/northing and longitude/easting at corner 1/2/3/4 (in degrees or meters), used in save_unavco, PyAPS (DEM file in radar coord), not accurate; number named in order of first line near/far range, last line near/far range.
+ WAVELENGTH = Radar wavelength in meters.
+ RANGE_PIXEL_SIZE = Slant range pixel size (search for pixel_ratio to convert to ground size, in m), used in dem_error, incidence_angle, multilook, transect.
+ EARTH_RADIUS = Best fitting spheroid radius in meters, used in dem_error, incidence_angle, convert2mat.
Expand All @@ -30,7 +30,7 @@ The following attributes vary for each interferogram:

+ ANTENNA_SIDE = -1 for right looking radar, used in save_unavco
+ AZIMUTH_PIXEL_SIZE = Azimuth pixel size at orbital altitude (multiply by Re/(Re+h) for ground size (m), where Re is the local earth radius), used in baseline_error/trop and multilook.
+ HEADING = Spacecraft heading at peg point (degree), measured from the north with clock-wise as positive, used in asc_desc
+ HEADING = Spacecraft heading at peg point (degrees), measured from the north with clock-wise as positive, used in asc_desc
+ PRF = Pulse repetition frequency (Hz), used in save_unavco

### Self-generated attributes ###
Expand All @@ -47,7 +47,8 @@ The following attributes vary for each interferogram:
+ NO_DATA_VALUE = No data value, value that should be ignored.
+ UNIT = data unit, i.e. m, m/yr, radian, and 1 for file without unit, such as coherence [[source]](https://github.com/insarlab/MintPy/blob/main/src/mintpy/objects/stack.py#L75)
+ REF_DATE = reference date
+ REF_X/Y/LAT/LON = column/row/latitude/longitude of reference point
+ REF_X/Y = column/row of the reference point
+ REF_LAT/LON = latitude/northing and longitude/easting of the reference point (in degrees or meters)
+ SUBSET_XMIN/XMAX/YMIN/YMAX = start/end column/row number of subset in the original coverage
+ 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)
Expand Down
9 changes: 5 additions & 4 deletions src/mintpy/prep_hyp3.py
Original file line number Diff line number Diff line change
Expand Up @@ -67,10 +67,6 @@ def add_hyp3_metadata(fname, meta, is_ifg=True):
S = N + float(meta['Y_STEP']) * int(meta['LENGTH'])
E = W + float(meta['X_STEP']) * int(meta['WIDTH'])

# convert UTM to lat/lon
N, W = ut.utm2latlon(meta, W, N)
S, E = ut.utm2latlon(meta, E, S)

if meta['ORBIT_DIRECTION'] == 'ASCENDING':
meta['LAT_REF1'] = str(S)
meta['LAT_REF2'] = str(S)
Expand Down Expand Up @@ -110,6 +106,11 @@ def add_hyp3_metadata(fname, meta, is_ifg=True):
meta['P_BASELINE_TOP_HDR'] = hyp3_meta['Baseline']
meta['P_BASELINE_BOTTOM_HDR'] = hyp3_meta['Baseline']

# HDF-EOS5 metadata
if hyp3_meta['ReferenceGranule'].startswith('S1'):
meta['beam_mode'] = 'IW'
meta['unwrap_method'] = hyp3_meta['Unwrappingtype']

return(meta)


Expand Down
28 changes: 16 additions & 12 deletions src/mintpy/save_hdfeos5.py
Original file line number Diff line number Diff line change
Expand Up @@ -115,20 +115,24 @@ def metadata_mintpy2unavco(meta_in, dateList, geom_file):
unavco_meta['last_date'] = dt.datetime.strptime(dateList[-1], '%Y%m%d').isoformat()[0:10]

# footprint
lons = [meta['LON_REF1'],
meta['LON_REF3'],
meta['LON_REF4'],
meta['LON_REF2'],
meta['LON_REF1']]

lats = [meta['LAT_REF1'],
meta['LAT_REF3'],
meta['LAT_REF4'],
meta['LAT_REF2'],
meta['LAT_REF1']]
lons = [float(meta['LON_REF1']),
float(meta['LON_REF3']),
float(meta['LON_REF4']),
float(meta['LON_REF2']),
float(meta['LON_REF1'])]

lats = [float(meta['LAT_REF1']),
float(meta['LAT_REF3']),
float(meta['LAT_REF4']),
float(meta['LAT_REF2']),
float(meta['LAT_REF1'])]

# convert UTM to lat/lon
if not meta_in.get('Y_UNIT', 'degrees').lower().startswith('deg'):
lats, lons = ut.utm2latlon(meta_in, easting=lons, northing=lats)

unavco_meta['scene_footprint'] = "POLYGON((" + ",".join(
[lon+' '+lat for lon, lat in zip(lons, lats)]) + "))"
[f'{lon} {lat}' for lon, lat in zip(lons, lats)]) + "))"

unavco_meta['history'] = dt.datetime.utcnow().isoformat()[0:10]

Expand Down
10 changes: 6 additions & 4 deletions src/mintpy/tropo_pyaps3.py
Original file line number Diff line number Diff line change
Expand Up @@ -275,10 +275,8 @@ def get_bounding_box(meta, geom_file=None):
lat1 = lat0 + lat_step * (length - 1)
lon1 = lon0 + lon_step * (width - 1)

# 'Y_FIRST' not in 'degree'
# e.g. meters for UTM projection from ASF HyP3
y_unit = meta.get('Y_UNIT', 'degrees').lower()
if not y_unit.startswith('deg'):
# for UTM projection, e.g. ASF HyP3
if not meta.get('Y_UNIT', 'degrees').lower().startswith('deg'):
lat0, lon0 = ut.utm2latlon(meta, easting=lon0, northing=lat0)
lat1, lon1 = ut.utm2latlon(meta, easting=lon1, northing=lat1)

Expand All @@ -300,8 +298,12 @@ def get_bounding_box(meta, geom_file=None):
lon1 = np.nanmax(lons)

else:
# use the rough (not accurate) lat/lon info of the four corners
lats = [float(meta[f'LAT_REF{i}']) for i in [1,2,3,4]]
lons = [float(meta[f'LON_REF{i}']) for i in [1,2,3,4]]
# for UTM projection, e.g. ASF HyP3
if not meta.get('Y_UNIT', 'degrees').lower().startswith('deg'):
lats, lons = ut.utm2latlon(meta, easting=lons, northing=lats)
lat0 = np.mean(lats[0:2])
lat1 = np.mean(lats[2:4])
lon0 = np.mean(lons[0:3:2])
Expand Down
33 changes: 23 additions & 10 deletions src/mintpy/utils/utils0.py
Original file line number Diff line number Diff line change
Expand Up @@ -345,10 +345,10 @@ def utm2latlon(meta, easting, northing):
Parameters: meta - dict, mintpy attributes that includes:
UTM_ZONE
easting - scalar or 1/2D np.ndarray, UTM coordinates in x direction
northing - scalar or 1/2D np.ndarray, UTM coordinates in y direction
Returns: lat - scalar or 1/2D np.ndarray, WGS 84 coordinates in y direction
lon - scalar or 1/2D np.ndarray, WGS 84 coordinates in x direction
easting - scalar/list/tuple/1-2D np.ndarray, UTM coordinates in x direction
northing - scalar/list/tuple/1-2D np.ndarray, UTM coordinates in y direction
Returns: lat - scalar/list/tuple/1-2D np.ndarray, WGS 84 coordinates in y direction
lon - scalar/list/tuple/1-2D np.ndarray, WGS 84 coordinates in x direction
"""
import utm
zone_num = int(meta['UTM_ZONE'][:-1])
Expand All @@ -357,20 +357,33 @@ def utm2latlon(meta, easting, northing):
# which can be common for large area analysis, e.g. the Norwegian mapping authority
# publishes a height data in UTM zone 33 coordinates for the whole country, even though
# most of it is technically outside zone 33.
lat, lon = utm.to_latlon(easting, northing, zone_num, northern=northern, strict=False)
lat, lon = utm.to_latlon(np.array(easting), np.array(northing), zone_num,
northern=northern, strict=False)

# output format
if any(isinstance(x, (list, tuple)) for x in [easting, northing]):
lat = lat.tolist()
lon = lon.tolist()

return lat, lon


def latlon2utm(lat, lon):
"""Convert latitude/longitude in degrees to UTM easting/northing in meters.
Parameters: lat - scalar or 1/2D np.ndarray, WGS 84 coordinates in y direction
lon - scalar or 1/2D np.ndarray, WGS 84 coordinates in x direction
Returns: easting - scalar or 1/2D np.ndarray, UTM coordinates in x direction
northing - scalar or 1/2D np.ndarray, UTM coordinates in y direction
Parameters: lat - scalar/list/tuple/1-2D np.ndarray, WGS 84 coordinates in y direction
lon - scalar/list/tuple/1-2D np.ndarray, WGS 84 coordinates in x direction
Returns: easting - scalar/list/tuple/1-2D np.ndarray, UTM coordinates in x direction
northing - scalar/list/tuple/1-2D np.ndarray, UTM coordinates in y direction
"""
import utm
easting, northing = utm.from_latlon(lat, lon)[:2]
easting, northing = utm.from_latlon(np.array(lat), np.array(lon))[:2]

# output format
if any(isinstance(x, (list, tuple)) for x in [lat, lon]):
easting = easting.tolist()
northing = northing.tolist()

return northing, easting


Expand Down
4 changes: 4 additions & 0 deletions tests/configs/RidgecrestSenDT71.txt
Original file line number Diff line number Diff line change
Expand Up @@ -19,4 +19,8 @@ mintpy.troposphericDelay.method = pyaps
mintpy.topographicResidual = yes
mintpy.topographicResidual.stepFuncDate = 20190706T0320
mintpy.topographicResidual.pixelwiseGeometry = no
mintpy.save.hdfEos5 = yes
mintpy.plot.maxMemory = 2

##————————————————————————————— HDF-EOS5 Attributes -—————————————————————————##
relative_orbit = 71

0 comments on commit 51fe70e

Please sign in to comment.