read/write/resample: improved support for GDAL complex file (#1000)

+ `readfile`:
    - `read_gdal()`: add `.cos` to the know list (for TerraSAR-X SLC data)
    - `read_binary_file()`: support returning complex data using `datasetName=complex`.
    - `read_binary_file()`: move `cpx_band` checking to expand its applicability from isce2 to generic files, including GDAL complex files.

+ `resample`: adjust size calculation in block-by-block IO for complex data type

+ `save_gdal`:
   - refactor `array2raster()` to `write_gdal()`, to be consistent with `utils.writefile.write()` in terms of input args.
   - use `readfile.DATA_TYPE_NUMPY2GDAL` to support non-float32 types, such as complex64

+ docs:
    - update Yunjun's contact email from Gmail to Outlook
    - README: simplify styles for Yunjun et al. (2019)
    - CONTRIBUTING: add links to the to-do list and roadmap

Author: yunjunz

docs/CODE_OF_CONDUCT.md

@@ -55,7 +55,7 @@ further defined and clarified by project maintainers.
 ## Enforcement
 
 Instances of abusive, harassing, or otherwise unacceptable behavior may be
-reported by contacting the project team at `yunjunzgeo at gmail dot com` or `hersh.fattahi at gmail dot com`. All
+reported by contacting the project team at `yunjunz at outlook dot com` or `hersh.fattahi at gmail dot com`. All
 complaints will be reviewed and investigated and will result in a response that
 is deemed necessary and appropriate to the circumstances. The project team is
 obligated to maintain confidentiality with regard to the reporter of an incident.

docs/CONTRIBUTING.md

@@ -108,3 +108,7 @@ It takes about 15 mins to finish.
 ## Things you should NOT do ##
 
 (For anyone with push rights to github.com/insarlab/MintPy) Never modify a commit or the history of anything that has been committed to the `main` branch.
+
+## Looking for ideas to contribute? ##
+
+Feel ready and look for something to try? Check our [to-do list](https://github.com/insarlab/MintPy/wiki/To-do-list) and [roadmap](https://github.com/insarlab/MintPy/wiki/version-2-roadmap) for the next major version 2.0!

docs/README.md

@@ -6,7 +6,7 @@
 [![Version](https://img.shields.io/github/v/release/insarlab/MintPy?color=yellow&label=version&style=flat-square)](https://github.com/insarlab/MintPy/releases)
 [![Forum](https://img.shields.io/badge/forum-Google%20Groups-orange.svg?style=flat-square)](https://groups.google.com/g/mintpy)
 [![License](https://img.shields.io/badge/license-GPLv3+-blue.svg?style=flat-square)](https://github.com/insarlab/MintPy/blob/main/LICENSE)
-[![Citation](https://img.shields.io/badge/DOI-10.1016%2Fj.cageo.2019.104331-blue?style=flat-square)](https://doi.org/10.1016/j.cageo.2019.104331)
+[![Citation](https://img.shields.io/badge/doi-10.1016%2Fj.cageo.2019.104331-blue?style=flat-square)](https://doi.org/10.1016/j.cageo.2019.104331)
 
 ## MintPy ##
 
@@ -110,6 +110,6 @@ _This disclaimer was adapted from the [MetPy project](https://github.com/Unidata
 
 ### 6. Citing this work ###
 
-Yunjun, Z., Fattahi, H., and Amelung, F. (2019), Small baseline InSAR time series analysis: Unwrapping error correction and noise reduction, _Computers & Geosciences_, _133_, 104331, doi:[10.1016/j.cageo.2019.104331](https://doi.org/10.1016/j.cageo.2019.104331), [arXiv](https://eartharxiv.org/9sz6m/), [data](https://doi.org/10.5281/zenodo.3464190), [notebooks](https://github.com/geodesymiami/Yunjun_et_al-2019-MintPy).
+Yunjun, Z., Fattahi, H., and Amelung, F. (2019), Small baseline InSAR time series analysis: Unwrapping error correction and noise reduction, _Computers & Geosciences_, _133_, 104331. [ [doi](https://doi.org/10.1016/j.cageo.2019.104331) \| [arxiv](https://doi.org/10.31223/osf.io/9sz6m) \| [data](https://doi.org/10.5281/zenodo.3464190) \| [notebook](https://github.com/geodesymiami/Yunjun_et_al-2019-MintPy) ]
 
 In addition to the above, we recommend that you cite the original publications that describe the algorithms used in your specific analysis. They are noted briefly in the [default template file](../mintpy/defaults/smallbaselineApp.cfg) and listed in the [references.md file](./references.md).

docs/_config.yml

@@ -5,5 +5,5 @@ google_analytics: UA-104225904-1
 theme: jekyll-theme-cayman
 author:
   name: Zhang Yunjun, Heresh Fattahi
-  url: https://yunjunzhang.wordpress.com
-author_email: yunjunzgeo@gmail.com
+  url: https://yunjunz.github.io/
+author_email: yunjunz@outlook.com

docs/dask.md

@@ -5,7 +5,7 @@ Most computations in MintPy are operated in either a pixel-by-pixel or a epoch-b
 + **local cluster:** on a single machine (laptop or computing node) with multiple CPU cores, suitable for laptops, local cluster/stations and distributed High Performance Cluster (HPC). No job scheduler is required.
 + **non-local cluster:** on a distributed HPC with job scheduler installed, including PBS, LSF and SLURM.
 
-Below is a brief description of the required options and recommended best practices for each cluster/scheduler.
+[Here](https://github.com/2gotgrossman/dask-rsmas-presentation) is an entry-level presentation on parallel computing using Dask by David Grossman. Below we brief describe for each cluster/scheduler the required options and recommended best practices.
 
 ## 1. local cluster ##
 

setup.py

@@ -29,7 +29,7 @@
     url="https://github.com/insarlab/MintPy",
     download_url=(f"https://github.com/insarlab/MintPy/archive/v{version}.tar.gz"),
     author="Zhang Yunjun, Heresh Fattahi",
-    author_email="yunjunzgeo@gmail.com",
+    author_email="yunjunz@outlook.com",
     license="GPL-3.0-or-later",
     license_files=("LICENSE",),
 

src/mintpy/cli/save_gdal.py

@@ -15,6 +15,7 @@
   save_gdal.py geo/geo_timeseries_ERA5_demErr.h5  -d 20200505_20200517 --of ENVI
   save_gdal.py geo/geo_ifgramStack.h5 -d unwrapPhase-20101120_20110220 --of ISCE
   save_gdal.py geo/geo_ifgramStack.h5 -d   coherence-20101120_20110220 --of ISCE
+  save_gdal.py geo_20230225.slc
 """
 
 

src/mintpy/objects/resample.py

@@ -252,18 +252,22 @@ def get_num_box(src_file=None, max_memory=0, scale_fac=3.0):
         # auto split into list of boxes if max_memory > 0
         if src_file and os.path.isfile(src_file) and max_memory > 0:
             # get max dataset shape
+            atr = readfile.read_attribute(src_file)
             fext = os.path.splitext(src_file)[1]
             if fext in ['.h5', '.he5']:
                 with h5py.File(src_file, 'r') as f:
                     ds_shapes = [f[i].shape for i in f.keys()
                                  if isinstance(f[i], h5py.Dataset)]
                     max_ds_size = max(np.prod(i) for i in ds_shapes)
             else:
-                atr = readfile.read_attribute(src_file)
                 max_ds_size = int(atr['LENGTH']) * int(atr['WIDTH'])
 
+            # scale the size for complex arrays
+            if atr.get('DATA_TYPE', 'float32').startswith('complex'):
+                max_ds_size *= 6 if atr['DATA_TYPE'].endswith('128') else 3
+
             # calc num_box
-            num_box = int(np.ceil((max_ds_size * 4 * 4) / (max_memory * 1024**3)))
+            num_box = int(np.ceil((max_ds_size * 4 * 6) / (max_memory * 1024**3)))
 
         return num_box
 

src/mintpy/save_gdal.py

@@ -24,37 +24,66 @@
 
 
 ##############################################################################
-def array2raster(array, out_file, transform, epsg, out_fmt='GTiff'):
-    """Write 2D matrix into gdal raster.
-
-    Parameters: array     - 2D np.ndarray
-                out_file  - str, output file path
-                transform - tuple(float), geotransform to connect image- to the geo-coordinates.
-                            https://gdal.org/tutorials/geotransforms_tut.html
-                epsg      - int, EPSG code for the Spatial Reference System (SRS)
-                out_fmt   - str, gdal driver name
-    Returns:    out_file  - str, output file path
+def write_gdal(data, meta, out_file, out_fmt='GTiff'):
+    """Write 2D matrix into GDAL raster.
+
+    Parameters: data     - 2D np.ndarray
+                meta     - dict, metadata used to calculate the transform / epsg info:
+                           X/Y_FIRST/STEP
+                           EPSG / UTM_ZONE
+                out_file - str, output file path
+                out_fmt  - str, output GDAL file format (driver name)
+                           https://gdal.org/drivers/raster/index.html
+    Returns:    out_file - str, output file path
     """
 
+    # prepare geotransform to connect image- to the geo-coordinates.
+    # link: https://gdal.org/tutorials/geotransforms_tut.html
+    transform = (
+        float(meta['X_FIRST']), float(meta['X_STEP']), 0,
+        float(meta['Y_FIRST']), 0, float(meta['Y_STEP']),
+    )
+
+    # prepare EPSG code for the Spatial Reference System (SRS)
+    if 'EPSG' in meta.keys():
+        epsg = int(meta['EPSG'])
+    elif 'UTM_ZONE' in meta.keys():
+        epsg = int(ut.utm_zone2epsg_code(meta['UTM_ZONE']))
+    else:
+        epsg = 4326
+        msg = 'No EPSG or UTM_ZONE metadata found! '
+        msg += 'Assume EPSG = 4326 (WGS84) and continue.'
+        warnings.warn(msg)
+
+    # write file
     driver = gdal.GetDriverByName(out_fmt)
     print(f'initiate GDAL driver: {driver.LongName}')
 
-    rows, cols = array.shape
-    print('create raster band')
-    print(f'raster row / column number: {rows}, {cols}')
-    print(f'raster transform info: {transform}')
-    raster = driver.Create(out_file, cols, rows, 1, gdal.GDT_Float32)
-    raster.SetGeoTransform(transform)
+    rows, cols = data.shape
+    dtype = readfile.DATA_TYPE_NUMPY2GDAL[str(data.dtype)]
+    print('create raster band:')
+    print(f'  raster row / column number: {rows}, {cols}')
+    print(f'  raster data type: {dtype} ({data.dtype})')
+    raster = driver.Create(
+        out_file,
+        xsize=cols,
+        ysize=rows,
+        bands=1,
+        eType=dtype,
+    )
 
-    print('write data to raster band')
-    band = raster.GetRasterBand(1)
-    band.WriteArray(array)
+    print(f'set transform info: {transform}')
+    raster.SetGeoTransform(transform)
 
     print(f'set projection as: EPSG {epsg}')
     srs = osr.SpatialReference()
     srs.ImportFromEPSG(epsg)
     raster.SetProjection(srs.ExportToWkt())
 
+    print('write data to raster band')
+    band = raster.GetRasterBand(1)
+    band.WriteArray(data)
+
     band.FlushCache()
     print(f'finished writing to {out_file}')
 
@@ -74,44 +103,20 @@ def save_gdal(inps):
 
     ds_name = inps.dset if inps.dset else 'data'
     print(f'read {ds_name} from file: {inps.file}')
-    data, atr = readfile.read(inps.file, datasetName=inps.dset)
+    data, meta = readfile.read(inps.file, datasetName=inps.dset)
 
     if ftype == 'timeseries' and inps.ref_date:
         print(f'read {inps.ref_date} from file: {inps.file}')
         data -= readfile.read(inps.file, datasetName=inps.ref_date)[0]
 
-    ## prepare the output
+    ## write file
     # output file name
     if not inps.outfile:
         fbase = pp.auto_figure_title(inps.file, inps.dset, vars(inps))
-        inps.outfile = fbase + GDAL_DRIVER2EXT.get(inps.out_format, '')
-    else:
-        inps.outfile = os.path.abspath(inps.outfile)
-
-    # geotransform
-    transform = (
-        float(atr['X_FIRST']), float(atr['X_STEP']), 0,
-        float(atr['Y_FIRST']), 0, float(atr['Y_STEP']),
-    )
+        fext = GDAL_DRIVER2EXT.get(inps.out_format, '')
+        inps.outfile = fbase + fext
+    inps.outfile = os.path.abspath(inps.outfile)
 
-    # epsg
-    if 'EPSG' in atr.keys():
-        epsg = int(atr['EPSG'])
-    elif 'UTM_ZONE' in atr.keys():
-        epsg = int(ut.utm_zone2epsg_code(atr['UTM_ZONE']))
-    else:
-        epsg = 4326
-        msg = 'No EPSG or UTM_ZONE metadata found! '
-        msg += 'Assume EPSG = 4326 (WGS84) and continue.'
-        warnings.warn(msg)
-
-    ## write gdal raster
-    array2raster(
-        data,
-        out_file=inps.outfile,
-        transform=transform,
-        epsg=epsg,
-        out_fmt=inps.out_format,
-    )
+    write_gdal(data, meta, out_file=inps.outfile, out_fmt=inps.out_format)
 
     return

src/mintpy/utils/readfile.py

@@ -119,6 +119,7 @@
 }
 
 # 2 - GDAL
+# link: https://gdal.org/api/raster_c_api.html#_CPPv412GDALDataType
 DATA_TYPE_GDAL2NUMPY = {
     1 : 'uint8',
     2 : 'uint16',
@@ -131,21 +132,25 @@
     9 : 'cint32',       # for translation purpose only, as numpy does not support complex int
     10: 'complex64',
     11: 'complex128',
+    12: 'uint64',
+    13: 'int64',
 }
 
 DATA_TYPE_NUMPY2GDAL = {
-    "uint8"     : 1,
-    "int8"      : 1,
-    "uint16"    : 2,
-    "int16"     : 3,
-    "uint32"    : 4,
-    "int32"     : 5,
-    "float32"   : 6,
-    "float64"   : 7,
-    "cint16"    : 8,    # for translation purpose only, as numpy does not support complex int
-    "cint32"    : 9,    # for translation purpose only, as numpy does not support complex int
-    "complex64" : 10,
-    "complex128": 11,
+    "uint8"     : 1,    # GDT_Byte
+    "int8"      : 1,    # GDT_Int8   (GDAL >= 3.7)
+    "uint16"    : 2,    # GDT_UInt16
+    "int16"     : 3,    # GDT_Int16
+    "uint32"    : 4,    # GDT_UInt32
+    "int32"     : 5,    # GDT_Int32
+    "float32"   : 6,    # GDT_Float32
+    "float64"   : 7,    # GDT_Float64
+    "cint16"    : 8,    # GDT_CInt16, for translation purpose only, as numpy does not support complex int
+    "cint32"    : 9,    # GDT_CInt32, for translation purpose only, as numpy does not support complex int
+    "complex64" : 10,   # GDT_CFloat32
+    "complex128": 11,   # GDT_CFloat64
+    "uint64"    : 12,   # GDT_UInt64 (GDAL >= 3.5)
+    "int64"     : 13,   # GDT_Int64  (GDAL >= 3.5)
 }
 
 # 3 - ISCE
@@ -173,7 +178,8 @@
 }
 
 # single file (data + attributes) supported by GDAL
-GDAL_FILE_EXTS = ['.tiff', '.tif', '.grd']
+# .cos file - TerraSAR-X complex SAR data (https://gdal.org/drivers/raster/cosar.html)
+GDAL_FILE_EXTS = ['.tiff', '.tif', '.grd', '.cos']
 
 ENVI_BAND_INTERLEAVE = {
     'BAND' : 'BSQ',
@@ -515,7 +521,20 @@ def read_binary_file(fname, datasetName=None, box=None, xstep=1, ystep=1):
 
     # default data to read
     band = 1
-    cpx_band = 'phase'
+    if datasetName:
+        if datasetName.startswith(('mag', 'amp')):
+            cpx_band = 'magnitude'
+        elif datasetName in ['phase', 'angle']:
+            cpx_band = 'phase'
+        elif datasetName.lower() == 'real':
+            cpx_band = 'real'
+        elif datasetName.lower().startswith('imag'):
+            cpx_band = 'imag'
+        else:
+            cpx_band = 'complex'
+    else:
+        # use phase as default value, since it's the most common one.
+        cpx_band = 'phase'
 
     # ISCE
     if processor in ['isce']:
@@ -553,18 +572,6 @@ def read_binary_file(fname, datasetName=None, box=None, xstep=1, ystep=1):
                 band = 2
             elif datasetName.lower() == 'band3':
                 band = 3
-
-            elif datasetName.startswith(('mag', 'amp')):
-                cpx_band = 'magnitude'
-            elif datasetName in ['phase', 'angle']:
-                cpx_band = 'phase'
-            elif datasetName.lower() == 'real':
-                cpx_band = 'real'
-            elif datasetName.lower().startswith('imag'):
-                cpx_band = 'imag'
-            elif datasetName.startswith(('cpx', 'complex')):
-                cpx_band = 'complex'
-
             else:
                 # flexible band list
                 ds_list = get_slice_list(fname)
@@ -766,6 +773,8 @@ def get_hdf5_2d_dataset(name, obj):
     # Binary Files
     else:
         num_band = int(atr.get('BANDS', '1'))
+        dtype = atr.get('DATA_TYPE', 'float32')
+
         if fext in ['.trans', '.utm_to_rdc']:
             # roipac / gamma lookup table
             slice_list = ['rangeCoord', 'azimuthCoord']
@@ -777,7 +786,7 @@ def get_hdf5_2d_dataset(name, obj):
         elif fext in ['.unw', '.ion']:
             slice_list = ['magnitude', 'phase']
 
-        elif fext in ['.int', '.slc']:
+        elif fext in ['.int', '.slc'] or (dtype.startswith('c') and num_band == 1):
             if no_complex:
                 slice_list = ['magnitude', 'phase']
             else:
@@ -1979,8 +1988,8 @@ def read_gdal(fname, box=None, band=1, cpx_band='phase', xstep=1, ystep=1):
             data = data.imag
         elif cpx_band.startswith('pha'):
             data = np.angle(data)
-        elif cpx_band.startswith('mag'):
-            data = np.absolute(data)
+        elif cpx_band.startswith(('mag', 'amp')):
+            data = np.abs(data)
         elif cpx_band.startswith(('cpx', 'complex')):
             pass
         else: