Sb/refpixel lat lon

input_parameters.conf

@@ -16,7 +16,7 @@ apsest:     0
 tscal:      1
 
 # Optional save of numpy array files for output products (MERGE)
-savenpy:    1
+savenpy:    0
 
 #%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 # Multi-threading parameters used by stacking/timeseries/prepifg
@@ -88,8 +88,8 @@ nan_conversion: 1
 # refnx/y: number of search grid points in x/y image dimensions
 # refchipsize: size of the data window at each search grid point
 # refminfrac: minimum fraction of valid (non-NaN) pixels in the data window
-refx:          
-refy:         
+refx:          150.941666654
+refy:          -34.218333314
 refnx:         5
 refny:         5
 refchipsize:   5

pyrate/core/ifgconstants.py

@@ -64,6 +64,14 @@
 DATA_UNITS = 'DATA_UNITS'
 INPUT_TYPE = 'INPUT_TYPE'
 
+PYRATE_REFPIX_X = 'REF_PIX_X'
+PYRATE_REFPIX_Y = 'REF_PIX_Y'
+PYRATE_REFPIX_LAT = 'REF_PIX_LAT'
+PYRATE_REFPIX_LON = 'REF_PIX_LON'
+PYRATE_MEAN_REF_AREA = 'REF_AREA_MEAN'
+PYRATE_STDDEV_REF_AREA = 'REF_AREA_STDDEV'
+SEQUENCE_POSITION = 'SEQUENCE_POSITION'
+
 DAYS_PER_YEAR = 365.25  # span of year, not a calendar year
 YEARS_PER_DAY = 1 / DAYS_PER_YEAR
 SPEED_OF_LIGHT_METRES_PER_SECOND = 3e8

pyrate/core/refpixel.py

@@ -25,12 +25,104 @@
 from numpy import isnan, std, mean, sum as nsum
 from joblib import Parallel, delayed
 
+from osgeo import gdal
+
 import pyrate.core.config as cf
+from pyrate.core import ifgconstants as ifc
+from pyrate.core import mpiops
 from pyrate.core.shared import Ifg
 from pyrate.core.shared import joblib_log_level
 from pyrate.core.logger import pyratelogger as log
 
 
+def update_refpix_metadata(ifg_paths, refx, refy, transform, params):
+    """
+    Function that adds metadata about the chosen reference pixel to each interferogram.
+    """
+
+    pyrate_refpix_lon, pyrate_refpix_lat = mpiops.run_once(convert_pixel_value_to_geographic_coordinate,
+                                                           refx, refy, transform)
+
+    process_ifgs_paths = mpiops.array_split(ifg_paths)
+
+    for ifg_file in process_ifgs_paths:
+        log.debug("Updating metadata for: "+ifg_file)
+        ifg = Ifg(ifg_file)
+        log.debug("Open dataset")
+        ifg.open(readonly=True)
+        log.debug("Set no data value")
+        ifg.nodata_value = params["noDataValue"]
+        log.debug("Update no data values in dataset")
+        ifg.convert_to_nans()
+        log.debug("Convert mm")
+        ifg.convert_to_mm()
+        half_patch_size = params["refchipsize"] // 2
+        x, y = refx, refy
+        log.debug("Extract reference pixel windows")
+        data = ifg.phase_data[y - half_patch_size: y + half_patch_size + 1,
+                              x - half_patch_size: x + half_patch_size + 1]
+        log.debug("Calculate standard deviation for reference window")
+        stddev_ref_area = np.nanstd(data)
+        log.debug("Calculate mean for reference window")
+        mean_ref_area = np.nanmean(data)
+        ifg.add_metadata(**{
+            ifc.PYRATE_REFPIX_X: str(refx),
+            ifc.PYRATE_REFPIX_Y: str(refy),
+            ifc.PYRATE_REFPIX_LAT: str(pyrate_refpix_lat),
+            ifc.PYRATE_REFPIX_LON: str(pyrate_refpix_lon),
+            ifc.PYRATE_MEAN_REF_AREA: str(mean_ref_area),
+            ifc.PYRATE_STDDEV_REF_AREA: str(stddev_ref_area)
+        })
+
+        ifg.close()
+
+
+def convert_pixel_value_to_geographic_coordinate(refx, refy, transform):
+    """
+    Converts a pixel coordinate to a latitude/longitude coordinate given the
+    geotransform of the image.
+    Args:
+        refx: The pixel x coordinate.
+        refy: The pixel ye coordinate.
+        transform: The geotransform array of the image.
+    Returns:
+        Tuple of lon, lat geographic coordinate.
+    """
+
+    xOrigin = transform[0]
+    yOrigin = transform[3]
+    pixelWidth = transform[1]
+    pixelHeight = -transform[5]
+
+    lon = refx*pixelWidth + xOrigin
+    lat = yOrigin - refy*pixelHeight
+
+    return lon, lat
+
+
+def convert_geographic_coordinate_to_pixel_value(lon, lat, transform):
+    """
+    Converts a latitude/longitude coordinate to a pixel coordinate given the
+    geotransform of the image.
+    Args:
+        lon: Pixel longitude.
+        lat: Pixel latitude.
+        transform: The geotransform array of the image.
+    Returns:
+        Tuple of refx, refy pixel coordinates.
+    """
+
+    xOrigin = transform[0]
+    yOrigin = transform[3]
+    pixelWidth = transform[1]
+    pixelHeight = -transform[5]
+
+    refx = int((lon - xOrigin) / pixelWidth)
+    refy = int((yOrigin - lat) / pixelHeight)
+
+    return refx + 1, refy + 1
+
+
 # TODO: move error checking to config step (for fail fast)
 # TODO: this function is not used. Plan removal
 def ref_pixel(ifgs, params):

pyrate/core/shared.py

@@ -479,6 +479,14 @@ def write_modified_phase(self, data=None):
             self.dataset.SetMetadataItem(k, v)
         self.dataset.FlushCache()
 
+    def add_metadata(self, **kwargs):
+        if (not self.is_open) or self.is_read_only:
+            raise IOError("Ifg not open or readonly. Cannot write!")
+
+        for k, v in kwargs.items():
+            self.dataset.SetMetadataItem(k, v)
+        self.dataset.FlushCache()  # write to disc
+
 
 class IfgPart(object):
     """
@@ -942,21 +950,13 @@ def write_output_geotiff(md, gt, wkt, data, dest, nodata):
 
     # set other metadata
     ds.SetMetadataItem('DATA_TYPE', str(md['DATA_TYPE']))
+
     # sequence position for time series products
-    if "SEQUENCE_POSITION" in md:
-        ds.SetMetadataItem("SEQUENCE_POSITION", str(md["SEQUENCE_POSITION"]))
-    if "PYRATE_REFPIX_LAT" in md:
-        ds.SetMetadataItem("PYRATE_REFPIX_LAT", str(md["PYRATE_REFPIX_LAT"]))
-    if "PYRATE_REFPIX_LON" in md:
-        ds.SetMetadataItem("PYRATE_REFPIX_LON", str(md["PYRATE_REFPIX_LON"]))
-    if "PYRATE_REFPIX_X" in md:
-        ds.SetMetadataItem("PYRATE_REFPIX_X", str(md["PYRATE_REFPIX_X"]))
-    if "PYRATE_REFPIX_Y" in md:
-        ds.SetMetadataItem("PYRATE_REFPIX_Y", str(md["PYRATE_REFPIX_Y"]))
-    if "PYRATE_MEAN_REF_AREA" in md:
-        ds.SetMetadataItem("PYRATE_MEAN_REF_AREA", str(md["PYRATE_MEAN_REF_AREA"]))
-    if "STANDARD_DEVIATION_REF_AREA" in md:
-        ds.SetMetadataItem("PYRATE_STANDARD_DEVIATION_REF_AREA", str(md["PYRATE_STANDARD_DEVIATION_REF_AREA"]))
+
+    for k in [ifc.SEQUENCE_POSITION, ifc.PYRATE_REFPIX_X, ifc.PYRATE_REFPIX_Y, ifc.PYRATE_REFPIX_LAT,
+              ifc.PYRATE_REFPIX_LON, ifc.PYRATE_MEAN_REF_AREA, ifc.PYRATE_STDDEV_REF_AREA]:
+        if k in md:
+            ds.SetMetadataItem(k, str(md[k]))
 
     # write data to geotiff
     band = ds.GetRasterBand(1)

pyrate/process.py

@@ -124,13 +124,15 @@ def _ref_pixel_calc(ifg_paths, params):
     """
     Wrapper for reference pixel calculation
     """
-    refx = params[cf.REFX]
-    refy = params[cf.REFY]
+    lon = params[cf.REFX]
+    lat = params[cf.REFY]
 
     ifg = Ifg(ifg_paths[0])
     ifg.open(readonly=True)
+    # assume all interferograms have same projection and will share the same transform
+    transform = ifg.dataset.GetGeoTransform()
 
-    if refx == -1 or refy == -1:
+    if lon == -1 or lat == -1:
 
         log.info('Searching for best reference pixel location')
 
@@ -150,12 +152,19 @@ def _ref_pixel_calc(ifg_paths, params):
                 "Reference pixel calculation returned an all nan slice!\n"
                 "Cannot continue downstream computation. Please change reference pixel algorithm used before "
                 "continuing.")
+        refy, refx = refpixel_returned   # row first means first value is latitude
+        log.info('Selected reference pixel coordinate (x, y): ({}, {})'.format(refx, refy))
+        lon, lat = refpixel.convert_pixel_value_to_geographic_coordinate(refx, refy, transform)
+        log.info('Selected reference pixel coordinate (lon, lat): ({}, {})'.format(lon, lat))
 
-        refy, refx = refpixel_returned
-
-        log.info('Selected reference pixel coordinate: ({}, {})'.format(refx, refy))
     else:
-        log.info('Reusing reference pixel from config file: ({}, {})'.format(refx, refy))
+        log.info('Using reference pixel from config file (lon, lat): ({}, {})'.format(lon, lat))
+        log.warning("Ensure user supplied reference pixel values are in lon/lat")
+        refx, refy = refpixel.convert_geographic_coordinate_to_pixel_value(lon, lat, transform)
+        log.info('Converted reference pixel coordinate (x, y): ({}, {})'.format(refx, refy))
+
+    refpixel.update_refpix_metadata(ifg_paths, refx, refy, transform, params)
+
     log.debug("refpx, refpy: "+str(refx) + " " + str(refy))
     ifg.close()
     return int(refx), int(refy)

tests/common.py

@@ -160,6 +160,26 @@ def assert_tifs_equal(tif1, tif2):
     sds = None
 
 
+def copy_small_ifg_file_list():
+    temp_dir = tempfile.mkdtemp()
+    move_files(SML_TEST_TIF, temp_dir, file_type='*.tif', copy=True)
+    datafiles = glob.glob(join(temp_dir, "*.tif"))
+    for d in datafiles:
+        Path(d).chmod(0o664)  # assign write permission as conv2tif output is readonly
+    return temp_dir, datafiles
+
+
+def copy_and_setup_small_data():
+    temp_dir, datafiles = copy_small_ifg_file_list()
+    datafiles.sort()
+    ifgs = [dem_or_ifg(i) for i in datafiles]
+
+    for i in ifgs:
+        i.open()
+        i.nodata_value = 0
+    return temp_dir, ifgs
+
+
 def small_ifg_file_list(datafiles=None):
     """Returns the file list of all the .tif files after prepifg conversion
     input phase data is in radians; these ifgs are in radians - not converted to mm"""

tests/test_data/small_test/conf/pyrate_gamma_test.conf

@@ -62,8 +62,8 @@ ifgylast:     -34.22
 # refnx/y: number of search grid points in x/y direction
 # refchipsize: chip size of the data window at each search grid point
 # refminfrac: minimum fraction of valid (non-NaN) pixels in the data window
-refx:
-refy:
+refx:          150.941666654
+refy:          -34.218333314
 refnx:         5
 refny:         5
 refchipsize:   5

tests/test_mpi_vs_multiprocess_vs_single_process.py

@@ -87,17 +87,11 @@ def test_pipeline_parallel_vs_mpi(modified_config, gamma_conf):
 
     mr_conf, params_m = modified_config(gamma_conf, 1, 'multiprocess_conf.conf')
 
-    check_call(f"pyrate conv2tif -f {mr_conf}", shell=True)
-    check_call(f"pyrate prepifg -f {mr_conf}", shell=True)
-    check_call(f"pyrate process -f {mr_conf}", shell=True)
-    check_call(f"pyrate merge -f {mr_conf}", shell=True)
+    check_call(f"pyrate workflow -f {mr_conf}", shell=True)
 
     sr_conf, params_s = modified_config(gamma_conf, 0, 'singleprocess_conf.conf')
 
-    check_call(f"pyrate conv2tif -f {sr_conf}", shell=True)
-    check_call(f"pyrate prepifg -f {sr_conf}", shell=True)
-    check_call(f"pyrate process -f {sr_conf}", shell=True)
-    check_call(f"pyrate merge -f {sr_conf}", shell=True)
+    check_call(f"pyrate workflow -f {sr_conf}", shell=True)
 
     # convert2tif tests, 17 interferograms
     assert_same_files_produced(params[cf.OUT_DIR], params_m[cf.OUT_DIR], params_s[cf.OUT_DIR], "*_unw_ifg.tif", 17)
@@ -231,11 +225,7 @@ def test_stack_and_ts_mpi_vs_parallel_vs_serial(modified_config_short, gamma_con
 
     sr_conf, params_p = modified_config_short(gamma_conf, parallel, 'parallel_conf.conf', 0)
 
-    check_call(f"pyrate conv2tif -f {sr_conf}", shell=True)
-    check_call(f"pyrate prepifg -f {sr_conf}", shell=True)
-    check_call(f"pyrate process -f {sr_conf}", shell=True)
-    check_call(f"pyrate merge -f {sr_conf}", shell=True)
-
+    check_call(f"pyrate workflow -f {sr_conf}", shell=True)
 
     # convert2tif tests, 17 interferograms
     assert_two_dirs_equal(params[cf.OUT_DIR], params_p[cf.OUT_DIR], "*_unw_ifg.tif", 17)
@@ -246,8 +236,7 @@ def test_stack_and_ts_mpi_vs_parallel_vs_serial(modified_config_short, gamma_con
         print("coherence files compared")
 
     # prepifg + process steps that overwrite tifs test
-
-    # 17 ifgs + 1 dem
+    # 17 mlooked ifgs + 1 dem + 17 mlooked coherence files
     if params[cf.COH_MASK]:
         assert_two_dirs_equal(params[cf.OUT_DIR], params_p[cf.OUT_DIR], f"*{params[cf.IFG_CROP_OPT]}cr.tif", 35)
     else:
@@ -262,7 +251,7 @@ def test_stack_and_ts_mpi_vs_parallel_vs_serial(modified_config_short, gamma_con
 
     # compare merge step
     assert_two_dirs_equal(params[cf.OUT_DIR], params_p[cf.OUT_DIR], "stack*.tif", 3)
-    #assert_two_dirs_equal(params[cf.OUT_DIR], params_p[cf.OUT_DIR], "stack*.npy", 3) not saved by default
+    assert_two_dirs_equal(params[cf.OUT_DIR], params_p[cf.OUT_DIR], "stack*.npy", 3)
     assert_two_dirs_equal(params[cf.OUT_DIR], params_p[cf.OUT_DIR], "tscuml*.tif")
 
     print("==========================xxx===========================")

tests/test_prepifg.py

@@ -87,6 +87,7 @@ def test_prepifg_file_types(tempdir, gamma_conf, coh_mask):
     cf.write_config_file(params=params, output_conf_file=output_conf)
     params_s = Configuration(output_conf).__dict__
     conv2tif.main(params_s)
+    # reread params from config
     params_s = Configuration(output_conf).__dict__
     prepifg.main(params_s)
     ifg_files = list(Path(tdir.joinpath(params_s[cf.OUT_DIR])).glob('*_ifg.tif'))
@@ -133,7 +134,7 @@ def test_prepifg_file_types(tempdir, gamma_conf, coh_mask):
     dem.open()
     md = dem.dataset.GetMetadata()
     assert md[ifc.DATA_TYPE] == ifc.MLOOKED_DEM
-    # shutil.rmtree(tdir)
+    shutil.rmtree(tdir)
 
 
 # convenience ifg creation funcs