Sb/ro prepifg

.travis.yml

@@ -79,7 +79,7 @@ deploy:
   verbose: true
   on:
     branch: master
-    condition: $GDALVERSION="3.0.2" && $TRAVIS_PYTHON_VERSION=3.8.*
+    python: 3.8
   github_token: $GITHUB_TOKEN2
   local_dir: docs/_build/html
   project_name: PyRate

pyrate/conv2tif.py

@@ -22,6 +22,7 @@
 from typing import Tuple, List
 from joblib import Parallel, delayed
 import numpy as np
+from pathlib import Path
 
 from pyrate.core.prepifg_helper import PreprocessError
 from pyrate.core import shared, mpiops, config as cf, gamma, roipac
@@ -104,6 +105,7 @@ def _geotiff_multiprocessing(unw_path: MultiplePaths, params: dict) -> Tuple[str
         else:
             raise PreprocessError('Processor must be ROI_PAC (0) or GAMMA (1)')
         shared.write_fullres_geotiff(header, unw_path.unwrapped_path, dest, nodata=params[cf.NO_DATA_VALUE])
+        Path(dest).chmod(0o444)  # readonly output
         return dest, True
     else:
         log.warning(f"Full-res geotiff already exists in {dest}! Returning existing geotiff!")

pyrate/core/gdal_python.py

@@ -276,9 +276,9 @@ def _gdalwarp_width_and_height(max_x, max_y, min_x, min_y, geo_trans):
 
 
 def crop_resample_average(
-        input_tif, extents: Union[List, Tuple], new_res, output_file, thresh,
-        out_driver_type='GTiff',
-        match_pyrate=False, hdr=None, coherence_path=None, coherence_thresh=None):
+        input_tif, extents: Union[List, Tuple], new_res, output_file, thresh, hdr, out_driver_type='GTiff',
+        match_pyrate=False, coherence_path=None, coherence_thresh=None
+        ):
     """
     Crop, resample, and average a geotiff image.
 
@@ -325,11 +325,8 @@ def crop_resample_average(
     gt = dst_ds.GetGeoTransform()
     wkt = dst_ds.GetProjection()
 
-    # TEST HERE IF EXISTING FILE HAS PYRATE METADATA. IF NOT ADD HERE
-    if ifc.DATA_TYPE not in dst_ds.GetMetadata() and hdr is not None:
-        md = shared.collate_metadata(hdr)
-    else:
-        md = dst_ds.GetMetadata()
+    # insert metadata from the header
+    md = shared.collate_metadata(hdr)
 
     # update metadata for output
 
@@ -346,11 +343,6 @@ def crop_resample_average(
                 md.update({ifc.DATA_TYPE: ifc.MLOOKED_DEM})
             elif v == ifc.INCIDENCE:
                 md.update({ifc.DATA_TYPE: ifc.MLOOKED_INC})
-            elif (v == ifc.COHERENCE) and (coherence_path is None):
-                # during orbital fit multilooking
-                pass
-            elif v == ifc.MULTILOOKED:
-                pass
             else:
                 raise TypeError(f'Data Type metadata {v} not recognised')
 

pyrate/core/orbital.py

@@ -68,7 +68,7 @@
 PART_CUBIC = cf.PART_CUBIC
 
 
-def remove_orbital_error(ifgs: Iterable, params: dict, preread_ifgs=None) -> None:
+def remove_orbital_error(ifgs: Iterable, params: dict, headers: List[dict], preread_ifgs=None) -> None:
     """
     Wrapper function for PyRate orbital error removal functionality.
 
@@ -96,7 +96,9 @@ def remove_orbital_error(ifgs: Iterable, params: dict, preread_ifgs=None) -> Non
             xlooks=params[cf.ORBITAL_FIT_LOOKS_X],
             ylooks=params[cf.ORBITAL_FIT_LOOKS_Y],
             thresh=params[cf.NO_DATA_AVERAGING_THRESHOLD],
-            write_to_disc=False)
+            write_to_disc=False,
+            headers=headers
+        )
         mlooked = [Ifg(m[1]) for m in mlooked_dataset]
 
         for m in mlooked:

pyrate/core/prepifg_helper.py

@@ -144,7 +144,7 @@ def _get_extents(ifgs, crop_opt, user_exts=None):
     return extents
 
 
-def prepare_ifg(raster_path, xlooks, ylooks, exts, thresh, crop_opt, write_to_disk=True, out_path=None, header=None,
+def prepare_ifg(raster_path, xlooks, ylooks, exts, thresh, crop_opt, header, write_to_disk=True, out_path=None,
                 coherence_path=None, coherence_thresh=None):
     """
     Open, resample, crop and optionally save to disk an interferogram or DEM.
@@ -202,7 +202,7 @@ def prepare_ifg(raster_path, xlooks, ylooks, exts, thresh, crop_opt, write_to_di
 
 
 # TODO: crop options 0 = no cropping? get rid of same size
-def prepare_ifgs(raster_data_paths, crop_opt, xlooks, ylooks, thresh=0.5, user_exts=None, write_to_disc=True,
+def prepare_ifgs(raster_data_paths, crop_opt, xlooks, ylooks, headers, thresh=0.5, user_exts=None, write_to_disc=True,
                  out_path=None):
     """
     Wrapper function to prepare a sequence of interferogram files for
@@ -228,8 +228,8 @@ def prepare_ifgs(raster_data_paths, crop_opt, xlooks, ylooks, thresh=0.5, user_e
     # use metadata check to check whether it's a dem or ifg
     rasters = [dem_or_ifg(r) for r in raster_data_paths]
     exts = get_analysis_extent(crop_opt, rasters, xlooks, ylooks, user_exts)
-
-    return [prepare_ifg(d, xlooks, ylooks, exts, thresh, crop_opt, write_to_disc, out_path) for d in raster_data_paths]
+    return [prepare_ifg(d, xlooks, ylooks, exts, thresh, crop_opt, h, write_to_disc, out_path) for d, h
+            in zip(raster_data_paths, headers)]
 
 
 # TODO: Not being used. Remove in future?

pyrate/core/shared.py

@@ -807,9 +807,9 @@ def collate_metadata(header):
             for k in [ifc.MASTER_TIME, ifc.SLAVE_TIME, ifc.PYRATE_INCIDENCE_DEGREES]:
                 md.update({k: str(header[k])})
     elif _is_incidence(header):
-        md.update({ifc.DATA_TYPE:ifc.INCIDENCE})
-    else: # must be dem
-        md.update({ifc.DATA_TYPE:ifc.DEM})
+        md.update({ifc.DATA_TYPE: ifc.INCIDENCE})
+    else:  # must be dem
+        md.update({ifc.DATA_TYPE: ifc.DEM})
 
     return md
 
@@ -1285,13 +1285,13 @@ def extract_epochs_from_filename(filename_with_epochs: str) -> List[str]:
 
 def mpi_vs_multiprocess_logging(step, params):
     if mpiops.size > 1:  # Over-ride input options if this is an MPI job
-        log.info(f"Running {step} step using mpi processing. Disabling parallel processing.")
+        log.info(f"Running {step} step using MPI processing. Disabling parallel processing.")
         params[cf.PARALLEL] = 0
     else:
         if params[cf.PARALLEL] == 1:
-            log.info(f"Running {step} using {params[cf.PROCESSES]} processes")
+            log.info(f"Running {step} step in parallel using {params[cf.PROCESSES]} processes")
         else:
-            log.info(f"Running {step} serially")
+            log.info(f"Running {step} step in serial")
 
 
 def dem_or_ifg(data_path):
@@ -1308,4 +1308,4 @@ def dem_or_ifg(data_path):
     if ifc.MASTER_DATE in md:  # ifg
         return Ifg(data_path)
     else:
-        return DEM(data_path)
+        return DEM(data_path)

pyrate/prepifg.py

@@ -206,14 +206,7 @@ def _prepifg_multiprocessing(path, xlooks, ylooks, exts, thresh, crop, params):
     """
     Multiprocessing wrapper for prepifg
     """
-    processor = params[cf.PROCESSOR]  # roipac, gamma or geotif
-    if (processor == GAMMA) or (processor == GEOTIF):
-        header = gamma.gamma_header(path, params)
-    elif processor == ROIPAC:
-        log.info("Warning: ROI_PAC support will be deprecated in a future PyRate release")
-        header = roipac.roipac_header(path, params)
-    else:
-        raise PreprocessError('Processor must be ROI_PAC (0) or GAMMA (1)')
+    header = find_header(path, params)
 
     # If we're performing coherence masking, find the coherence file for this IFG.
     if params[cf.COH_MASK] and shared._is_interferogram(header):
@@ -230,3 +223,15 @@ def _prepifg_multiprocessing(path, xlooks, ylooks, exts, thresh, crop, params):
     else:
         prepifg_helper.prepare_ifg(path, xlooks, ylooks, exts, thresh, crop, out_path=params[cf.OUT_DIR],
                                    header=header, coherence_path=coherence_path, coherence_thresh=coherence_thresh)
+
+
+def find_header(path, params):
+    processor = params[cf.PROCESSOR]  # roipac, gamma or geotif
+    if (processor == GAMMA) or (processor == GEOTIF):
+        header = gamma.gamma_header(path, params)
+    elif processor == ROIPAC:
+        log.info("Warning: ROI_PAC support will be deprecated in a future PyRate release")
+        header = roipac.roipac_header(path, params)
+    else:
+        raise PreprocessError('Processor must be ROI_PAC (0) or GAMMA (1)')
+    return header

pyrate/process.py

@@ -17,7 +17,6 @@
 """
 This Python module runs the main PyRate processing workflow
 """
-import logging
 import os
 from os.path import join
 import pickle as cp
@@ -31,6 +30,7 @@
 from pyrate.core.aps import wrap_spatio_temporal_filter
 from pyrate.core.shared import Ifg, PrereadIfg, get_tiles, mpi_vs_multiprocess_logging
 from pyrate.core.logger import pyratelogger as log
+from pyrate.prepifg import find_header
 
 MASTER_PROCESS = 0
 
@@ -186,7 +186,8 @@ def _orb_fit_calc(ifg_paths, params, preread_ifgs=None):
         # A performance comparison should be made for saving multilooked
         # files on disc vs in memory single process multilooking
         if mpiops.rank == MASTER_PROCESS:
-            orbital.remove_orbital_error(ifg_paths, params, preread_ifgs)
+            headers = [find_header(p, params) for p in ifg_paths]
+            orbital.remove_orbital_error(ifg_paths, params, headers, preread_ifgs=preread_ifgs)
     mpiops.comm.barrier()
     log.debug('Finished Orbital error correction')
 

tests/common.py

@@ -274,9 +274,12 @@ def reconstruct_mst(shape, tiles, output_dir):
     return mst
 
 
-def move_files(source_dir, dest_dir, file_type='*.tif'):
+def move_files(source_dir, dest_dir, file_type='*.tif', copy=False):
     for filename in glob.glob(os.path.join(source_dir, file_type)):
-        shutil.move(filename, dest_dir)
+        if copy:
+            shutil.copy(filename, dest_dir)
+        else:
+            shutil.move(filename, dest_dir)
 
 
 def assert_ifg_phase_equal(ifg_path1, ifg_path2):

tests/conftest.py

@@ -46,6 +46,11 @@ def fin():
     return mpiops.comm
 
 
+@pytest.fixture(params=[0, 1])
+def coh_mask(request):
+    return request.param
+
+
 @pytest.fixture(params=[1, 2])
 def ref_est_method(request):
     return request.param
@@ -107,7 +112,7 @@ def roipac_or_gamma_conf(request):
 def gamma_conf(request):
     params = Configuration(TEST_CONF_GAMMA).__dict__
     yield request.param
-    shutil.rmtree(params[cf.OUT_DIR])
+    # shutil.rmtree(params[cf.OUT_DIR])
 
 
 @pytest.fixture

tests/test_coherence.py

@@ -1,4 +1,5 @@
 import os
+import stat
 import tempfile
 import numpy as np
 from osgeo import osr
@@ -24,14 +25,14 @@ def test_small_data_coherence(gamma_params):
     conv2tif.main(gamma_params)
 
     for i in ifg_multilist:
-        p = i.converted_path
-        ifg = pyrate.core.shared.dem_or_ifg(data_path=p)
+        p = Path(i.converted_path)
+        p.chmod(0o664)  # assign write permission as conv2tif output is readonly
+        ifg = pyrate.core.shared.dem_or_ifg(data_path=p.as_posix())
         if not isinstance(ifg, Ifg):
             continue
         ifg.open()
-
         # now do coherence masking and compare
-        ifg = pyrate.core.shared.dem_or_ifg(data_path=p)
+        ifg = pyrate.core.shared.dem_or_ifg(data_path=p.as_posix())
         ifg.open()
         converted_coh_file_path = cf.coherence_paths_for(p, gamma_params, tif=True)
         gdal_python.coherence_masking(ifg.dataset,

tests/test_conv2tif.py

@@ -2,6 +2,7 @@
 import pytest
 import glob
 import copy
+from pathlib import Path
 
 import pyrate.core.config as cf
 from pyrate import conv2tif, prepifg
@@ -33,6 +34,8 @@ def test_num_gamma_tifs_equals_num_unws(gamma_params):
     # 17 unws + dem
     assert len(gtifs) == 18
 
+    for g, _ in gtifs:   # assert all output from conv2tfi are readonly
+        assert Path(g).stat().st_mode == 33060
 
 def test_num_roipac_tifs_equals_num_unws(roipac_params):
     gtifs = conv2tif.main(roipac_params)

tests/test_covariance.py

@@ -30,6 +30,8 @@
 from pyrate.core.covariance import cvd, get_vcmt, RDist
 from pyrate.configuration import Configuration
 import pyrate.core.orbital
+from pyrate.core import roipac
+from pyrate.core.config import parse_namelist
 from tests import common
 from tests.common import (small5_mock_ifgs, small5_ifgs, TEST_CONF_ROIPAC,
     small_data_setup, prepare_ifgs_without_phase)
@@ -179,9 +181,9 @@ class LegacyEqualityTest(unittest.TestCase):
 
     @classmethod
     def setUpClass(cls):
-
         params = Configuration(TEST_CONF_ROIPAC).__dict__
         cls.temp_out_dir = tempfile.mkdtemp()
+        common.move_files(common.SML_TEST_OBS, cls.temp_out_dir, file_type='*.rsc', copy=True)
         sys.argv = ['prepifg.py', TEST_CONF_ROIPAC]
         params[cf.OUT_DIR] = cls.temp_out_dir
         params[cf.TMPDIR] = os.path.join(cls.temp_out_dir, cf.TMPDIR)
@@ -191,12 +193,12 @@ def setUpClass(cls):
         prepifg.main(params)
         cls.params = params
         xlks, ylks, crop = cf.transform_params(params)
-        base_ifg_paths = cf.original_ifg_paths(params[cf.IFG_FILE_LIST],
-                                               params[cf.OBS_DIR])
+        base_ifg_paths = parse_namelist(params[cf.IFG_FILE_LIST])
         dest_paths = cf.get_dest_paths(base_ifg_paths, crop, params, xlks)
         ifgs = common.pre_prepare_ifgs(dest_paths, params)
         refx, refy = process._ref_pixel_calc(dest_paths, params)
-        pyrate.core.orbital.remove_orbital_error(ifgs, params)
+        headers = [roipac.roipac_header(i, cls.params) for i in base_ifg_paths]
+        pyrate.core.orbital.remove_orbital_error(ifgs, params, headers)
         ifgs = prepare_ifgs_without_phase(dest_paths, params)
         for ifg in ifgs:
             ifg.close()
@@ -205,8 +207,7 @@ def setUpClass(cls):
         r_dist = RDist(ifgs[0])()
         ifgs[0].close()
         # Calculate interferogram noise
-        cls.maxvar = [cvd(i, params, r_dist, calc_alpha=True,
-                          save_acg=True, write_vals=True)[0] for i in dest_paths]
+        cls.maxvar = [cvd(i, params, r_dist, calc_alpha=True, save_acg=True, write_vals=True)[0] for i in dest_paths]
         cls.vcmt = get_vcmt(ifgs, cls.maxvar)
         for ifg in ifgs:
             ifg.close()
@@ -226,8 +227,7 @@ def test_legacy_maxvar_equality_small_test_files(self):
     def test_legacy_vcmt_equality_small_test_files(self):
         from tests.common import SML_TEST_DIR
         LEGACY_VCM_DIR = os.path.join(SML_TEST_DIR, 'vcm')
-        legacy_vcm = np.genfromtxt(os.path.join(LEGACY_VCM_DIR,
-                                   'vcmt.csv'), delimiter=',')
+        legacy_vcm = np.genfromtxt(os.path.join(LEGACY_VCM_DIR, 'vcmt.csv'), delimiter=',')
         np.testing.assert_array_almost_equal(legacy_vcm, self.vcmt, decimal=3)
 
     def test_metadata(self):
@@ -243,8 +243,7 @@ def test_save_cvd_data(self):
             if not ifg.is_open:
                 ifg.open()
             data_file = join(self.params[cf.TMPDIR],
-                             'cvd_data_{b}.npy'.format(
-                                 b=basename(ifg.data_path).split('.')[0]))
+                             'cvd_data_{b}.npy'.format(b=basename(ifg.data_path).split('.')[0]))
             assert isfile(data_file)
 
 if __name__ == "__main__":

tests/test_data/merge/stack_rate.png.aux.xml

@@ -1,5 +1,5 @@
 <PAMDataset>
-  <SRS dataAxisToSRSAxisMapping="2,1">GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]</SRS>
+  <SRS>GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433],AUTHORITY["EPSG","4326"]]</SRS>
   <GeoTransform>  1.5091000000000000e+02,  8.3333300000000001e-04,  0.0000000000000000e+00, -3.4170000000000002e+01,  0.0000000000000000e+00, -8.3333300000000001e-04</GeoTransform>
   <Metadata domain="IMAGE_STRUCTURE">
     <MDI key="INTERLEAVE">PIXEL</MDI>