resample: refactor default output step size calculation (#997)

+ rename utils.utils0.auto_lat_lon_step_size() to objects.resample.resample.get_lat_lon_step()
   - move the location because this function is only used in the resample code so far
   - add the user-specified lalo_step code body into this function, for a more simple `prepare_geometry_definition_radar()`
   - replace try/except with if/else and explicit required metadata checking and warning msg

+ objects.resample.resample.prepare_geometry_definition_radar.find_valid_lat_lon():
   - rename mark_lat_lon_anomaly() to find_valid_lat_lon()
   - add more comments on the input/output args

+ cli.geocode: set default fillValue to zero for .int/unw files, to be consistent with isce/roipac.

+ docs/api/module_hierarchy.md: update for objects.resample.py

+ objects.progress.FileProgressObject(): display in int MB to be simple

docs/api/module_hierarchy.md

@@ -42,7 +42,7 @@ Hierarchy of sub-modules within MintPy. Level _N_ modules depends on level _N-1_
         s1_utils      (objects/{stack}, utils/{ptime, time_func})
 ------------------ level 3 --------------------
     /objects
-        resample      (utils/{utils0, ptime, readfile})
+        resample      (utils/{utils0, ptime, readfile}, constants)
         coord         (constants, utils/{utils0, utils1, readfile})
     /simulation
         iono          (utils/{utils0, readfile})

src/mintpy/cli/geocode.py

@@ -81,7 +81,7 @@ def create_parser(subparsers=None):
     interp.add_argument('-i', '--interp', dest='interpMethod', default='nearest', choices={'nearest', 'linear'},
                         help='interpolation/resampling method (default: %(default)s).')
     interp.add_argument('--fill', dest='fillValue', type=float, default=math.nan,
-                        help='Fill value for extrapolation (default: %(default)s).')
+                        help='Fill value for extrapolation (default: %(default)s or 0 for *.int/unw files).')
     interp.add_argument('-n','--nprocs', dest='nprocs', type=int, default=1,
                         help='number of processors to be used for calculation (default: %(default)s).\n'
                              'Note: Do not use more processes than available processor cores.')
@@ -109,6 +109,11 @@ def cmd_line_parse(iargs=None):
     # import
     from mintpy.utils import readfile, utils as ut
 
+    # save argv (to check the manually specified arguments)
+    # use iargs        for python call
+    # use sys.argv[1:] for command line call
+    inps.argv = iargs if iargs else sys.argv[1:]
+
     # check
     if inps.templateFile:
         inps = read_template2inps(inps.templateFile, inps)
@@ -169,6 +174,11 @@ def cmd_line_parse(iargs=None):
             print('ERROR: "--geo2radar" is NOT supported for "--software scipy"!')
             sys.exit(0)
 
+    # default: --fill (set default to zero for .int/unw file)
+    fext = os.path.splitext(inps.file[0])[1]
+    if '--fill' not in inps.argv and fext in ['.int', '.unw']:
+        inps.fillValue = 0
+
     return inps
 
 

src/mintpy/objects/progress.py

@@ -36,7 +36,7 @@ def __init__(self, path, *args, **kwargs):
 
     def read(self, size):
         perc = self.tell() / self._total_size * 100
-        msg = f"extracting: {self.tell()/(1024**2):.1f} of {self._total_size/(1024**2):.1f} MB - {perc:.0f}%"
+        msg = f"extracting: {self.tell()/(1024**2):.0f} of {self._total_size/(1024**2):.0f} MB - {perc:.0f}%"
         sys.stdout.write("\r" + msg)
         sys.stdout.flush()
         return io.FileIO.read(self, size)

src/mintpy/objects/resample.py

@@ -19,11 +19,10 @@
     warnings.simplefilter('ignore', UserWarning)
     import pyresample as pr
 
+from mintpy.constants import EARTH_RADIUS
 from mintpy.objects.cluster import split_box2sub_boxes
 from mintpy.utils import ptime, readfile, utils0 as ut
 
-EARTH_RADIUS = 6378122.65   # m
-
 
 class resample:
     """
@@ -283,14 +282,101 @@ def get_radius_of_influence(lalo_step, src_meta, ratio=3):
         return radius
 
 
+    def get_lat_lon_step(self, src_lat0, src_lat1, src_lon0, src_lon1):
+        """Get/check the lat/lon step size for geocoding.
+
+        Approach 1: ensure the same pixel area before / after resampling
+            Treat the pixel in radar coordinates as an rotated rectangle. Use the bounding
+                box of the rotated rectangle for the ratio between lat and lon steps. Then
+                scale the lat and lon step size to ensure the same area between the pixels
+                in radar and geo coordinates.
+            This is recommended, but it requires metadata on the pixel size info.
+            Link: https://math.stackexchange.com/questions/4001034
+
+        Approach 2: ensure the same matrix shape before / after resampling
+            This is the backup approach if approach 1 does not work.
+
+        Parameters: src_lat0/1   - float, max/min latitude  in degree of the input data file
+                    src_lon0/1   - float, min/max longitude in degree of the input data file
+        Returns:    lat/lon_step - float, output step size in latitude/longitude in degree
+        """
+
+        if self.lalo_step is not None:
+            # use input lat/lon step
+            lat_step, lon_step = self.lalo_step
+
+        else:
+            # calculate default lat/lon step
+            # approach 1: ensure the same pixel area
+            print('calculate output pixel size using approach 1 '
+                  '(same pixel area before/after resampling)')
+
+            # check required metadata
+            meta = {**self.lut_meta, **self.src_meta}
+            key_list = [
+                'AZIMUTH_PIXEL_SIZE',
+                'HEADING',
+                'HEIGHT',
+                'RANGE_PIXEL_SIZE',
+                'STARTING_RANGE',
+                'WIDTH',
+            ]
+            key_not_found = [x for x in key_list if x not in meta.keys()]
+            if len(key_not_found) == 0:
+
+                # azimuth angle (rotation angle) in radian
+                az_angle = np.deg2rad(abs(ut.heading2azimuth_angle(float(meta['HEADING']))))
+
+                # radar pixel size in meter
+                az_step = ut.azimuth_ground_resolution(meta)
+                rg_step = ut.range_ground_resolution(meta)
+
+                # geo pixel size in meter
+                x_step = rg_step * abs(np.cos(az_angle)) + az_step * abs(np.sin(az_angle))
+                y_step = rg_step * abs(np.sin(az_angle)) + az_step * abs(np.cos(az_angle))
+                scale_factor = np.sqrt((rg_step * az_step) / (x_step * y_step))
+                x_step *= scale_factor
+                y_step *= scale_factor
+
+                # geo pixel size in degree
+                lat_c = (src_lat0 + src_lat1) / 2.
+                lon_step = np.rad2deg(x_step / (EARTH_RADIUS * np.cos(np.deg2rad(lat_c))))
+                lat_step = np.rad2deg(y_step / EARTH_RADIUS) * -1.
+
+            else:
+                # approach 2: ensure the same matrix shape
+                msg = 'WARNING: NOT all required metadata are found. '
+                msg += f'Missing metadata: {key_not_found}. Switch to approach 2.\n'
+                msg += 'calculate output pixel size using approach 2 '
+                msg += '(same matrix shape before/after resampling)'
+                print(msg)
+
+                # ensure the same matrix shape before / after geocoding
+                # if not enough metadata found for the above
+                lut_len = int(self.lut_meta['LENGTH'])
+                lut_wid = int(self.lut_meta['WIDTH'])
+                lat_step = (src_lat1 - src_lat0) / (lut_len - 1)
+                lon_step = (src_lon1 - src_lon0) / (lut_wid - 1)
+
+        # ensure lat/lon step sign
+        lat_step = abs(lat_step) * -1.
+        lon_step = abs(lon_step)
+
+        return lat_step, lon_step
+
+
     ##------------------------------ resample using pyresample -------------------------------------##
 
     def prepare_geometry_definition_radar(self):
         """Get src_def and dest_def for lookup table in radar-coord (from ISCE, DORIS)"""
 
-        def mark_lat_lon_anomoly(lat, lon):
+        def find_valid_lat_lon(lat, lon):
             """mask pixels with abnormal values (0, etc.)
             This is found on sentinelStack multiple swath lookup table file.
+
+            Parameters: lat/lon - 2D np.ndarray in float32, latitude/longitude in degrees.
+            Returns:    lat/lon - 2D np.ndarray in float32, latitude/longitude in degrees.
+                        mask    - 2D np.ndarray in bool, 1/0 for valid/invalid pixels.
             """
             # ignore pixels with zero value
             zero_mask = np.multiply(lat != 0., lon != 0.)
@@ -327,7 +413,7 @@ def mark_lat_lon_anomoly(lat, lon):
         lon_file = self.lon_file if self.lon_file else self.lut_file
         lut_lat = readfile.read(lat_file, datasetName='latitude')[0].astype(np.float32)
         lut_lon = readfile.read(lon_file, datasetName='longitude')[0].astype(np.float32)
-        lut_lat, lut_lon, mask = mark_lat_lon_anomoly(lut_lat, lut_lon)
+        lut_lat, lut_lon, mask = find_valid_lat_lon(lut_lat, lut_lon)
 
         # radar2geo (with block-by-block support)
         if 'Y_FIRST' not in self.src_meta.keys():
@@ -339,24 +425,7 @@ def mark_lat_lon_anomoly(lat, lon):
             src_lon1 = np.nanmax(lut_lon[mask])
 
             # parameter 1 - lalo_step (output grid)
-            if self.lalo_step is None:
-                try:
-                    # ensure the same pixel area before / after geocoding
-                    merged_meta = {**self.lut_meta, **self.src_meta}
-                    lat_c = (src_lat0 + src_lat1) / 2.
-                    lat_step, lon_step = ut.auto_lat_lon_step_size(merged_meta, lat_c)
-
-                except KeyError:
-                    # ensure the same matrix shape before / after geocoding
-                    # if not enough metadata found for the above
-                    lat_step = (src_lat1 - src_lat0) / (lut_lat.shape[0] - 1)
-                    lon_step = (src_lon1 - src_lon0) / (lut_lat.shape[1] - 1)
-                self.lalo_step = (abs(lat_step) * -1., abs(lon_step))
-
-            else:
-                # ensure lat/lon step sign
-                self.lalo_step = (abs(self.lalo_step[0]) * -1.,
-                                  abs(self.lalo_step[1]) * 1.)
+            self.lalo_step = self.get_lat_lon_step(src_lat0, src_lat1, src_lon0, src_lon1)
             print(f'output pixel size in (lat, lon) in degree: {self.lalo_step}')
 
             # parameter 2 / 3 - SNWE (at pixel outer boundary; output grid) / length & width

src/mintpy/utils/utils0.py

@@ -24,9 +24,9 @@
 import numpy as np
 
 # global variables
-SPEED_OF_LIGHT = 299792458 # m/s
-EARTH_RADIUS = 6371e3      # Earth radius in meters
-K = 40.31                  # m^3/s^2, constant
+SPEED_OF_LIGHT = 299792458  # m/s
+EARTH_RADIUS = 6371.0088e3  # Earth radius in meters
+K = 40.31                   # m^3/s^2, constant
 
 
 #################################### InSAR ##########################################
@@ -125,7 +125,7 @@ def incidence_angle(atr, dem=None, dimension=2, print_msg=True):
     # Read Attributes
     range_n = float(atr['STARTING_RANGE'])
     dR = float(atr['RANGE_PIXEL_SIZE'])
-    r = float(atr['EARTH_RADIUS'])
+    r = float(atr.get('EARTH_RADIUS', EARTH_RADIUS))
     H = float(atr['HEIGHT'])
     width = int(atr['WIDTH'])
 
@@ -228,47 +228,6 @@ def azimuth_ground_resolution(atr):
     return az_step
 
 
-def auto_lat_lon_step_size(atr, lat_c=None):
-    """Get the default lat/lon step size for geocoding.
-
-    Treat the pixel in radar coordinates as an rotated rectangle. Use the bounding box
-    of the rotated rectangle for the ratio between lat and lon steps. Then scale the
-    lat and lon step size to ensure the same area between the pixels in radar and geo
-    coordinates.
-
-    Link: https://math.stackexchange.com/questions/4001034
-
-    Parameters: atr      - dict, standard mintpy metadata
-                lat_c    - float, central latitude in degree
-    Returns:    lat_step - float, latitude  step size in degree
-                lon_step - float, longitude step size in degree
-    """
-    # azimuth angle (rotation angle) in radian
-    az_angle = np.deg2rad(abs(heading2azimuth_angle(float(atr['HEADING']))))
-
-    # radar pixel size in meter
-    az_step = azimuth_ground_resolution(atr)
-    rg_step = range_ground_resolution(atr)
-
-    # geo pixel size in meter
-    x_step = rg_step * abs(np.cos(az_angle)) + az_step * abs(np.sin(az_angle))
-    y_step = rg_step * abs(np.sin(az_angle)) + az_step * abs(np.cos(az_angle))
-    scale_factor = np.sqrt((rg_step * az_step) / (x_step * y_step))
-    x_step *= scale_factor
-    y_step *= scale_factor
-
-    # geo pixel size in degree
-    if lat_c is None:
-        if 'LAT_REF1' in atr.keys():
-            lat_c = (float(atr['LAT_REF1']) + float(atr['LAT_REF3'])) / 2.
-        else:
-            lat_c = 0
-    lon_step = np.rad2deg(x_step / (EARTH_RADIUS * np.cos(np.deg2rad(lat_c))))
-    lat_step = np.rad2deg(y_step / EARTH_RADIUS) * -1.
-
-    return lat_step, lon_step
-
-
 
 #################################### File Operation ##########################################
 def touch(fname_list, times=None):