Merge pull request #41 from ASFHyP3/develop

Release color browse images and S2 patches

CHANGELOG.md

@@ -6,7 +6,15 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
 and this project adheres to [PEP 440](https://www.python.org/dev/peps/pep-0440/)
 and uses [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
-## [0.3.0](https://github.com/ASFHyP3/hyp3-autorift/compare/v0.1.0...v0.2.0)
+## [0.3.1](https://github.com/ASFHyP3/hyp3-autorift/compare/v0.3.0...v0.3.1)
+
+### Changed
+* Browse images and thumbnails will be generated using the [ITS_LIVE](https://its-live.jpl.nasa.gov/) colormap
+
+### Fixed
+* Failures due to stable surface misclassification
+
+## [0.3.0](https://github.com/ASFHyP3/hyp3-autorift/compare/v0.2.0...v0.3.0)
 
 ### Added
 * Installed autoRIFT v1.0.8 for processing optical scenes (in addition to v1.0.7 already installed as part of ISCE)

conda-env.yml

@@ -8,6 +8,7 @@ dependencies:
   - python=3.8
   - pip
   # For packaging, and testing
+  - pillow
   - pytest
   - pytest-console-scripts
   - pytest-cov
@@ -21,7 +22,7 @@ dependencies:
   - isce2=2.4.2
   - autorift=1.0.8
   - boto3
-  - importlib_metadata
+  - matplotlib-base
   - netCDF4
   - numpy
   - psycopg2  # missing hyp3proclib dep

hyp3_autorift/image.py

@@ -0,0 +1,49 @@
+from pathlib import Path
+from typing import Optional
+
+import matplotlib.pyplot as plt
+import numpy as np
+from matplotlib.colors import LinearSegmentedColormap
+from scipy.interpolate import PchipInterpolator
+
+COLOR_MAP = np.array([
+    # data value, R, G, B, A
+    [0, 255, 255, 255, 0],
+    [2, 166, 238, 255, 255],
+    [4, 97, 195, 219, 255],
+    [9, 84, 169, 254, 255],
+    [16, 84, 130, 254, 255],
+    [25, 84, 85, 254, 255],
+    [36, 50, 119, 220, 255],
+    [49, 16, 153, 186, 255],
+    [64, 16, 186, 153, 255],
+    [81, 50, 220, 119, 255],
+    [100, 84, 254, 85, 255],
+    [121, 118, 221, 51, 255],
+    [144, 153, 186, 16, 255],
+    [169, 187, 152, 17, 255],
+    [196, 221, 118, 51, 255],
+    [225, 255, 85, 85, 255],
+    [289, 255, 25, 85, 255],
+    [324, 213, 1, 72, 255],
+    [361, 158, 1, 66, 255],
+    [400, 140, 0, 51, 255],
+    [441, 122, 0, 166, 255],
+    [484, 140, 0, 191, 255],
+    [529, 159, 0, 217, 255],
+    [576, 213, 0, 255, 255],
+    [625, 255, 0, 138, 255],
+])
+
+
+def make_browse(out_file: Path, data: np.ndarray,
+                min_value: Optional[float] = None, max_value: Optional[float] = 625.) -> Path:
+    data_values = COLOR_MAP[:, 0]
+    pchip = PchipInterpolator(data_values, np.linspace(0, 1, len(data_values)))
+    image = pchip(np.clip(data, min_value, max_value))
+
+    rgb_values = COLOR_MAP[:, 1:] / 255
+    cmap = LinearSegmentedColormap.from_list('its-live', rgb_values)
+    plt.imsave(out_file, image, cmap=cmap)
+
+    return out_file

hyp3_autorift/netcdf_output.py

@@ -10,6 +10,13 @@
 import hyp3_autorift
 
 
+def v_error_cal(vx_error, vy_error):
+    vx = np.random.normal(0, vx_error, 1000000)
+    vy = np.random.normal(0, vy_error, 1000000)
+    v = np.sqrt(vx**2 + vy**2)
+    return np.std(v)
+
+
 def netCDF_packaging(VX, VY, DX, DY, INTERPMASK, CHIPSIZEX, CHIPSIZEY, SSM, SX, SY,
                      offset2vx_1, offset2vx_2, offset2vy_1, offset2vy_2, MM, VXref, VYref,
                      rangePixelSize, azimuthPixelSize, dt, epsg, srs, tran, out_nc_filename, pair_type,
@@ -19,11 +26,11 @@ def netCDF_packaging(VX, VY, DX, DY, INTERPMASK, CHIPSIZEX, CHIPSIZEY, SSM, SX,
         vx_mean_shift = offset2vx_1 * dx_mean_shift + offset2vx_2 * dy_mean_shift
         temp = vx_mean_shift
         temp[np.logical_not(SSM)] = np.nan
-        vx_mean_shift = np.median(temp[(temp > -500) & (temp < 500)])
+        vx_mean_shift = np.median(temp[np.logical_not(np.isnan(temp))])
         vy_mean_shift = offset2vy_1 * dx_mean_shift + offset2vy_2 * dy_mean_shift
         temp = vy_mean_shift
         temp[np.logical_not(SSM)] = np.nan
-        vy_mean_shift = np.median(temp[(temp > -500) & (temp < 500)])
+        vy_mean_shift = np.median(temp[np.logical_not(np.isnan(temp))])
     else:
         vx_mean_shift = 0.0
         vy_mean_shift = 0.0
@@ -125,7 +132,7 @@ def netCDF_packaging(VX, VY, DX, DY, INTERPMASK, CHIPSIZEX, CHIPSIZEY, SSM, SX,
     institution = 'NASA Jet Propulsion Laboratory (JPL), California Institute of Technology'
 
     isce_version = subprocess.check_output('conda list | grep isce | awk \'{print $2}\'', shell=True, text=True)
-    autorift_version = '1.0.7'
+    autorift_version = '1.0.8'
     source = f'ASF DAAC HyP3 {datetime.datetime.now().year} using the {hyp3_autorift.__name__} plugin version' \
              f' {hyp3_autorift.__version__} running autoRIFT version {autorift_version} as distributed with ISCE ' \
              f'version {isce_version.strip()}. Contains modified Copernicus Sentinel data ' \
@@ -332,6 +339,7 @@ def netCDF_packaging(VX, VY, DX, DY, INTERPMASK, CHIPSIZEX, CHIPSIZEY, SSM, SX,
     var[:] = np.round(np.clip(V, -32768, 32767)).astype(np.int16)
     var.setncattr('missing_value', np.int16(NoDataValue))
 
+    v_error = v_error_cal(vx_error, vy_error)
     varname = 'v_error'
     datatype = np.dtype('int16')
     dimensions = ('y', 'x')
@@ -346,6 +354,7 @@ def netCDF_packaging(VX, VY, DX, DY, INTERPMASK, CHIPSIZEX, CHIPSIZEY, SSM, SX,
         var.setncattr('description', 'velocity magnitude error')
     var.setncattr('units', 'm/y')
     V_error = np.sqrt((vx_error * VX / V) ** 2 + (vy_error * VY / V) ** 2)
+    V_error[V==0] = v_error
     V_error[noDataMask] = NoDataValue
     var[:] = np.round(np.clip(V_error, -32768, 32767)).astype(np.int16)
     var.setncattr('missing_value', np.int16(NoDataValue))

hyp3_autorift/process.py

@@ -17,11 +17,11 @@
 from hyp3lib.fetch import download_file
 from hyp3lib.file_subroutines import mkdir_p
 from hyp3lib.get_orb import downloadSentinelOrbitFile
-from hyp3lib.makeAsfBrowse import makeAsfBrowse
 from hyp3lib.scene import get_download_url
-from osgeo import gdal
+from netCDF4 import Dataset
 
 from hyp3_autorift import geometry
+from hyp3_autorift import image
 from hyp3_autorift import io
 
 log = logging.getLogger(__name__)
@@ -194,35 +194,22 @@ def process(reference: str, secondary: str, polarization: str = 'hh', band: str
             cmd = f'testautoRIFT.py -r {reference_url} -s {secondary_url} {autorift_parameters} -nc S2 -fo 1 -urlflag 1'
             execute(cmd, logfile=f, uselogging=True)
 
-    velocity_tif = gdal.Open('velocity.tif')
-    x_velocity = np.ma.masked_invalid(velocity_tif.GetRasterBand(1).ReadAsArray())
-    y_velocity = np.ma.masked_invalid(velocity_tif.GetRasterBand(2).ReadAsArray())
-    velocity = np.sqrt(x_velocity**2 + y_velocity**2)
-
-    browse_file = Path('velocity_browse.tif')
-    driver = gdal.GetDriverByName('GTiff')
-    browse_tif = driver.Create(
-        str(browse_file), velocity.shape[1], velocity.shape[0], 1, gdal.GDT_Byte, ['COMPRESS=LZW']
-    )
-    browse_tif.SetProjection(velocity_tif.GetProjection())
-    browse_tif.SetGeoTransform(velocity_tif.GetGeoTransform())
-    velocity_band = browse_tif.GetRasterBand(1)
-    velocity_band.WriteArray(velocity)
-
-    del velocity_band, browse_tif, velocity_tif
-
-    product_name = get_product_name(reference, secondary, orbit_files=(reference_state_vec, secondary_state_vec),
-                                    band=band)
-    makeAsfBrowse(str(browse_file), product_name)
-
     netcdf_files = glob.glob('*.nc')
     if not netcdf_files:
         raise Exception('Processing failed! Output netCDF file not found')
     if len(netcdf_files) > 1:
         log.warning(f'Too many netCDF files found; using first:\n    {netcdf_files}')
 
+    product_name = get_product_name(
+        reference, secondary, orbit_files=(reference_state_vec, secondary_state_vec), band=band
+    )
     product_file = Path(f'{product_name}.nc').resolve()
-    shutil.move(netcdf_files[0], f'{product_name}.nc')
+    shutil.move(netcdf_files[0], str(product_file))
+
+    with Dataset(product_file) as nc:
+        velocity = nc.variables['v']
+        data = np.ma.masked_values(velocity, -32767.).filled(0)
+    image.make_browse(product_file.with_suffix('.png'), data)
 
     return product_file