Skip to content

Commit

Permalink
convert coord btw. UTM and lat/lon (#1052)
Browse files Browse the repository at this point in the history
+ utils.utils0.py: add the following functions to convert coordinates between UTM and lat/lon:
   - utm2latlon()
   - latlon2utm()
   - utm_zone2epsg_code(): bugfix if the latitude letter is not N or S
   - get_lat_lon(): convert UTM to lat/lon to ensure lat/lon output for hyp3 products

+ add utm to the dependency files, since 1) it depends on numpy only; and 2) it's available on PyPI and conda-forge
   - requirements.txt
   - setup.py
  • Loading branch information
yunjunz authored Aug 1, 2023
1 parent fec802b commit 308024a
Show file tree
Hide file tree
Showing 3 changed files with 51 additions and 10 deletions.
1 change: 1 addition & 0 deletions requirements.txt
Original file line number Diff line number Diff line change
Expand Up @@ -23,3 +23,4 @@ pysolid
rich
scikit-image
scipy
utm
1 change: 1 addition & 0 deletions setup.py
Original file line number Diff line number Diff line change
Expand Up @@ -71,6 +71,7 @@
"setuptools",
"scikit-image",
"scipy",
"utm",
],
extras_require={
"cli": ["argcomplete"],
Expand Down
59 changes: 49 additions & 10 deletions src/mintpy/utils/utils0.py
Original file line number Diff line number Diff line change
Expand Up @@ -259,21 +259,25 @@ def touch(fname_list, times=None):
################################## Coordinate ##########################################
def utm_zone2epsg_code(utm_zone):
"""Convert UTM Zone string to EPSG code.
Parameters: utm_zone - str, atr['UTM_ZONE']
Returns: epsg - str, EPSG code
Examples: epsg = utm_zone2epsg_code('11N')
Parameters: utm_zone - str, atr['UTM_ZONE']
Returns: epsg_code - str, EPSG code
Examples: epsg_code = utm_zone2epsg_code('11N')
"""
from pyproj import CRS
crs = CRS.from_dict({'proj': 'utm',
'zone': int(utm_zone[:-1]),
'south': utm_zone[-1] == 'S',
})
epsg = crs.to_authority()[1]
return epsg

# N is the first letter in the northern hemisphere
crs = CRS.from_dict({
'proj': 'utm',
'zone': int(utm_zone[:-1]),
'south': utm_zone[-1].upper() < 'N',
})
epsg_code = crs.to_authority()[1]
return epsg_code


def to_latlon(infile, x, y):
"""Convert x, y in the projection coordinates of the file to lon/lat in degree.
"""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
Expand Down Expand Up @@ -302,6 +306,37 @@ def to_latlon(infile, x, y):
return y, x


def utm2latlon(meta, easting, northing):
"""Convert UTM easting/northing in meters to lat/lon in degrees.
Parameters: meta - dict, mintpy attributes that includes:
UTM_ZONE
easting - scalar or 1/2D np.ndarray, UTM coordiantes in x direction
northing - scalar or 1/2D np.ndarray, UTM coordiantes in y direction
Returns: lat - scalar or 1/2D np.ndarray, WGS 84 coordiantes in y direction
lon - scalar or 1/2D np.ndarray, WGS 84 coordiantes in x direction
"""
import utm
zone_num = int(meta['UTM_ZONE'][:-1])
lat_band = meta['UTM_ZONE'][-1].upper()
lat, lon = utm.to_latlon(easting, northing, zone_num, lat_band)
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 coordiantes in y direction
lon - scalar or 1/2D np.ndarray, WGS 84 coordiantes in x direction
Returns: easting - scalar or 1/2D np.ndarray, UTM coordiantes in x direction
northing - scalar or 1/2D np.ndarray, UTM coordiantes in y direction
zone_num - int, UTM zone number, from 1 to 60
lat_band - str, MGRS latitude band, C-X omitting I and O
"""
import utm
return utm.from_latlon(lat, lon)


def snwe_to_wkt_polygon(snwe):
"""Convert the input bounding box in SNWE into WKT format POLYGON.
Expand Down Expand Up @@ -369,6 +404,10 @@ def get_lat_lon(meta, geom_file=None, box=None, dimension=2, ystep=1, xstep=1):
else:
raise ValueError(f'un-supported dimension = {dimension}')

# UTM to lat/lon
if not meta['Y_UNIT'].startswith('deg') and 'UTM_ZONE' in meta.keys():
lats, lons = utm2latlon(meta, easting=lons, northing=lats)

else:
msg = 'Can not get pixel-wise lat/lon!'
msg += '\nmeta dict is not geocoded and/or geometry file does not contains latitude/longitude dataset.'
Expand Down

0 comments on commit 308024a

Please sign in to comment.