updated README; renamed dsc__CoReg_Sat_FourierShiftTheorem.py to coreg_cmd.py

components.CoReg.COREG:
- edited docstring

- renamed dsc__CoReg_Sat_FourierShiftTheorem.py to coreg_cmd.py
- updated __version__
- updated README and associated figures

README.md

@@ -9,6 +9,8 @@ Prerequisites and hints:
 The input images can have any GDAL compatible image format (http://www.gdal.org/formats_list.html). Both of them must be georeferenced. In case of ENVI files, this means they must have a 'map info' and a 'coordinate system string' as attributes of their header file. Different projection systems are currently not supported. The input images must have a geographic overlap but clipping them to same geographical extent is NOT neccessary. Please do not perform any spatial resampling of the input images before applying this algorithm. Any needed resampling of the data is done automatically. Thus the input images can have different spatial resolutions. The current algorithm will not perform any ortho-rectification. So please use ortho-rectified input data in order to prevent local shifts in the output image. By default the calculated subpixel-shifts are applied to the header file of the output image. No spatial resampling is done automatically as long as the both input images have the same projection. If you need the output image to be aligned to the reference image coordinate grid (by using an appropriate resampling algorithm), use the '-align_grids' option. The image overlap area is automatically calculated. Thereby no-data regions within the images are standardly respected. Providing the map coordinates of the actual data corners lets you save some calculation time, because in this case the automatic algorithm can be skipped. The no-data value of each image is automatically derived from the image corners. The verbose program mode gives some more output about the interim results, shows some figures and writes the used footprint and overlap polygons to disk. The figures must be manually closed in in order to continue the processing.
 
 
+ 
+
 # Install
 
 For now, there is no automatic install script. Just clone the repository, install the dependencies and add the root directory of CoReg_Sat to your PATH environment variable.
@@ -31,6 +33,8 @@ In addition clone the repository "py_tools_ds" and add its root directory to you
     
 Since py_tools_ds is not a public repository right now, contact Daniel Scheffler if you can not access it.
 
+ 
+
 # Modules
 
 ## CoReg
@@ -45,8 +49,8 @@ This module calculates spatial shifts and performs a global correction (based on
 ```python
 from CoReg_Sat import COREG
 
-im_reference = '/path/to/your/ref_image.bsq'
-im_target    = '/path/to/your/tgt_image.bsq'
+#im_reference = '/path/to/your/ref_image.bsq'
+#im_target    = '/path/to/your/tgt_image.bsq'
 
 CR = COREG(im_reference, im_target, wp=(354223, 5805559), ws=(256,256))
 CR.calculate_spatial_shifts()
@@ -76,8 +80,8 @@ CR.calculate_spatial_shifts()
 from py_tools_ds.ptds import GeoArray
 from CoReg_Sat import COREG
 
-#im_reference = '/path/to/your/ref_image.bsq'
-#im_target    = '/path/to/your/tgt_image.bsq'
+im_reference = '/path/to/your/ref_image.bsq'
+im_target    = '/path/to/your/tgt_image.bsq'
 
 # get a sample numpy array with corresponding geoinformation as reference image
 geoArr  = GeoArray(im_reference)
@@ -174,7 +178,9 @@ To write the coregistered image to disk, the COREG class needs to be instanced w
 #### apply detected shifts to multiple images
 
 Sometimes it can be useful to apply the same shifts to multiple images - e.g. to different mask images derived from the same satellite dataset.
-For this purpose you can calculate spatial shifts using the COREG class (see above) and then apply the calculated shifts to mulitple images using the DESHIFTER class:
+For this purpose you can calculate spatial shifts using the COREG class (see above) and then apply the calculated shifts to mulitple images using the DESHIFTER class.
+
+Take a look at the keyword arguments of the DESHIFTER class when you need further adjustments (e.g. output paths for the corrected images; aligned output grid, ...).
 
 
 ```python
@@ -217,7 +223,7 @@ DESHIFTER(im_target2, CR.coreg_info).correct_shifts()
 
 
 
-Take a look at the keyword arguments of the DESHIFTER class when you need further adjustments (e.g. output paths for the corrected images; aligned output grid, ...).
+ 
 
 ### Shell console interface
 
@@ -236,6 +242,10 @@ Follow these instructions to run CoReg_Sat from a shell console. For example, th
 python ./dsc__CoReg_Sat_FourierShiftTheorem.py /path/to/your/ref_image.bsq /path/to/your/tgt_image.bsq
 ```
 
+ 
+
+ 
+
 ## CoReg_local
 
 This module has been designed to detect and correct geometric shifts present locally in your input image. The class COREG_LOCAL calculates a grid of spatial shifts with points spread over the whole overlap area of the input images. Based on this grid a correction of local shifts can be performed.
@@ -270,9 +280,12 @@ CRL.correct_shifts()
     	[[319460.0, 5790510.0], [352270.0, 5900040.0], [409790.0, 5900040.0], [409790.0, 5790250.0], [319460.0, 5790250.0]]
     Matching window position (X,Y): 372220.10753674706/5841066.947109019
     Calculating geometric quality grid (1977 points) in mode 'multiprocessing'...
+    	progress: |==================================================| 100.0% [1977/1977] Complete 9.75 sek
+    Found 1144 valid GCPs.
+    Correcting geometric shifts...
     Translating progress |==================================================| 100.0% Complete
     Warping progress     |==================================================| 100.0% Complete
-    Writing GeoArray of size (10979, 9033) to /home/gfz-fe/scheffler/jupyter/CoReg_Sat_jupyter/my_project/S2A_OPER_MSI_L1C_TL_SGS__20160608T153121_A005024_T33UUU_B03__shifted_to__S2A_OPER_MSI_L1C_TL_SGS__20160529T153631_A004881_T33UUU_B03.bsq.
+    Writing GeoArray of size (10979, 10979) to /home/gfz-fe/scheffler/jupyter/CoReg_Sat_jupyter/my_project/S2A_OPER_MSI_L1C_TL_SGS__20160608T153121_A005024_T33UUU_B03__shifted_to__S2A_OPER_MSI_L1C_TL_SGS__20160529T153631_A004881_T33UUU_B03.bsq.
 
 
 
@@ -285,7 +298,7 @@ CRL.correct_shifts()
                   ['UTM',
                    1,
                    1,
-                   319460.0,
+                   300000.0,
                    5900030.0,
                    10.0,
                    10.0,
@@ -293,18 +306,18 @@ CRL.correct_shifts()
                    'North',
                    'WGS-84']),
                  ('updated geotransform',
-                  [319460.0, 10.0, 0.0, 5900030.0, 0.0, -10.0]),
+                  [300000.0, 10.0, 0.0, 5900030.0, 0.0, -10.0]),
                  ('updated projection',
                   'PROJCS["WGS 84 / UTM zone 33N",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"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",15],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32633"]]'),
-                 ('arr_shifted', array([[   0,    0,    0, ...,  952,  984, 1058],
-                         [   0,    0,    0, ...,  997,  976, 1037],
-                         [   0,    0,    0, ...,  960,  991, 1030],
+                 ('arr_shifted', array([[   0,    0,    0, ..., 1034,  996, 1001],
+                         [   0,    0,    0, ..., 1046, 1114, 1124],
+                         [   0,    0,    0, ..., 1021, 1126, 1148],
                          ..., 
-                         [   0,    0,    0, ...,  760,  769,  804],
-                         [   0,    0,    0, ...,  762,  754,  765],
+                         [   0,    0,    0, ...,  760,  769,  805],
+                         [   0,    0,    0, ...,  762,  755,  765],
                          [   0,    0,    0, ...,    0,    0,    0]], dtype=uint16)),
                  ('GeoArray_shifted',
-                  <py_tools_ds.ptds.io.raster.GeoArray.GeoArray at 0x7f0592b297f0>)])
+                  <py_tools_ds.ptds.io.raster.GeoArray.GeoArray at 0x7f451ac14a90>)])
 
 
 
@@ -335,12 +348,14 @@ CRL.view_CoRegPoints(figsize=(15,15),backgroundIm='ref')
 
 
 
-![png](output_36_1.png)
+![png](output_40_1.png)
 
 
+The output figure shows the calculated absolute lenghts of the shift vectors - in this case with shifts up to ~25 meters.
+
 #### visualize geometric quality grid with shifts present AFTER shift correction
 
-The remaining shifts after local correction can be visualized by instanciating COREG_LOCAL with the output path of the above instance of COREG_LOCAL.
+The remaining shifts after local correction can be calculated and visualized by instanciating COREG_LOCAL with the output path of the above instance of COREG_LOCAL.
 
 
 ```python
@@ -353,16 +368,19 @@ CRL_after_corr.view_CoRegPoints(figsize=(15,15),backgroundIm='ref')
     	[[319090.0, 5790510.0], [351800.0, 5899940.0], [409790.0, 5900040.0], [409790.0, 5790250.0], [319090.0, 5790250.0]]
     Calculating actual data corner coordinates for image to be shifted...
     Corner coordinates of image to be shifted:
-    	[[319460.0, 5790540.0], [352280.0, 5900030.0], [409780.0, 5900030.0], [409780.0, 5790260.0], [322940.0, 5790250.0], [319460.0, 5790280.0]]
-    Matching window position (X,Y): 372215.77905147866/5841065.066319922
+    	[[319460.0, 5790540.0], [352270.0, 5900030.0], [409780.0, 5900030.0], [409780.0, 5790260.0], [322970.0, 5790250.0], [319460.0, 5790280.0]]
+    Matching window position (X,Y): 372216.38593955856/5841068.390957352
     Note: array has been downsampled to 1000 x 1000 for faster visualization.
-    Calculating geometric quality grid (2014 points) in mode 'multiprocessing'...
+    Calculating geometric quality grid (1977 points) in mode 'multiprocessing'...
+    	progress: |==================================================| 100.0% [1977/1977] Complete 10.78 sek
 
 
 
-![png](output_39_1.png)
+![png](output_44_1.png)
 
 
+The output figure shows a significant reduction of geometric shifts.
+
 #### show the points table of the calculated geometric quality grid
 
 NOTE: Point records where no valid match has been found are filled with -9999.
@@ -1448,6 +1466,8 @@ CRL.CoRegPoints_table
 CRL.quality_grid.to_PointShapefile(path_out='/path/to/your/output_shapefile.shp')
 ```
 
+ 
+
 ### Shell console interface
 
 By far, there is no shell console interface for this module.

__init__.py

@@ -9,7 +9,7 @@ from .components import utilities
 from .components import geometry
 
 __author__ = 'Daniel Scheffler'
-__version__= '2016-11-03_01'
+__version__= '2016-11-03_02'
 
 __all__=['COREG',
          'COREG_LOCAL',

components/CoReg.py

@@ -110,12 +110,17 @@ class imParamObj(object):
 
 
 class COREG(object):
+    """See help(COREG) for documentation!"""
+
     def __init__(self, im_ref, im_tgt, path_out=None, fmt_out='ENVI', r_b4match=1, s_b4match=1, wp=(None,None),
                  ws=(512, 512), max_iter=5, max_shift=5, align_grids=False, match_gsd=False, out_gsd=None,
                  resamp_alg_deshift='cubic', resamp_alg_calc='cubic', data_corners_im0=None, data_corners_im1=None,
                  nodata=(None,None), calc_corners=True, multiproc=True, binary_ws=True, force_quadratic_win=True,
                  progress=True, v=False, path_verbose_out=None, q=False, ignore_errors=False):
-        """
+
+        """Detects and corrects global X/Y shifts between a target and refernce image. Geometric shifts are calculated
+        at a specific (adjustable) image position. Correction performs a global shifting in X- or Y direction.
+
         :param im_ref(str, GeoArray):   source path (any GDAL compatible image format is supported) or GeoArray instance
                                         of reference image
         :param im_tgt(str, GeoArray):   source path (any GDAL compatible image format is supported) or GeoArray instance

components/CoReg_local.py

@@ -22,6 +22,8 @@ from py_tools_ds.ptds                 import GeoArray
 
 
 class COREG_LOCAL(object):
+    """See help(COREG_LOCAL) for documentation!"""
+
     def __init__(self, im_ref, im_tgt, grid_res, window_size=(256,256), path_out=None, fmt_out='ENVI', projectDir=None,
                  r_b4match=1, s_b4match=1, max_iter=5, max_shift=5, data_corners_im0=None,
                  data_corners_im1=None, outFillVal=-9999, nodata=(None,None), calc_corners=True, binary_ws=True,

components/DeShifter.py

@@ -23,12 +23,13 @@ _dict_rspAlg_rsp_Int = {'nearest': 0, 'bilinear': 1, 'cubic': 2, 'cubic_spline':
                         'mode': 6, 'max': 7, 'min': 8 , 'med': 9, 'q1':10, 'q2':11}
 
 class DESHIFTER(object):
+    """See help(DESHIFTER) for documentation!"""
     def __init__(self, im2shift, coreg_results, **kwargs):
         """
         Deshift an image array or one of its products by applying the coregistration info calculated by COREG class.
 
         :param im2shift:            <path,GeoArray> path of an image to be de-shifted or alternatively a GeoArray object
-        :param coreg_results:       <dict> the results of the co-registration as given by COREG.coreg_info  or
+        :param coreg_results:       <dict> the results of the co-registration as given by COREG.coreg_info or
                                     COREG_LOCAL.coreg_info respectively
 
         :Keyword Arguments:
@@ -313,4 +314,28 @@ class DESHIFTER(object):
         deshift_results.update({'updated projection'  : self.updated_projection})
         deshift_results.update({'arr_shifted'         : self.arr_shifted})
         deshift_results.update({'GeoArray_shifted'    : self.GeoArray_shifted})
-        return deshift_results
\ No newline at end of file
+        return deshift_results
+
+
+
+def deshift_image_using_coreg_info(im2shift, coreg_results, path_out=None, fmt_out='ENVI', q=False):
+    """Corrects a geometrically distorted image using previously calculated coregistration info. This function can be
+    used for example to corrects spatial shifts of mask files using the same transformation parameters that have been
+    used to correct their source images.
+
+    :param im2shift:      <path,GeoArray> path of an image to be de-shifted or alternatively a GeoArray object
+    :param coreg_results: <dict> the results of the co-registration as given by COREG.coreg_info or
+                          COREG_LOCAL.coreg_info respectively
+    :param path_out:      /output/directory/filename for coregistered results. If None, no output is written - only
+                          the shift corrected results are returned.
+    :param fmt_out:       raster file format for output file. ignored if path_out is None. can be any GDAL
+                                        compatible raster file format (e.g. 'ENVI', 'GeoTIFF'; default: ENVI)
+    :param q:             quiet mode (default: False)
+    :return:
+    """
+    deshift_results = DESHIFTER(im2shift, coreg_results).correct_shifts()
+
+    if path_out:
+        deshift_results['GeoArray_shifted'].save(path_out, fmt_out=fmt_out, q=q)
+
+    return deshift_results

components/Geom_Quality_Grid.py

@@ -31,6 +31,8 @@ global_shared_im2shift = None
 
 
 class Geom_Quality_Grid(object):
+    """See help(Geom_Quality_Grid) for documentation!"""
+
     def __init__(self, COREG_obj, grid_res, outFillVal=-9999, dir_out=None, CPUs=None, progress=True, v=False, q=False):
 
         """Applies the algorithm to detect spatial shifts to the whole overlap area of the input images. Spatial shifts

dsc__CoReg_Sat_FourierShiftTheorem.py → coreg_cmd.py

@@ -6,6 +6,7 @@ __author__ = "Daniel Scheffler"
 import time
 import sys
 import os
+import warnings
 
 # custom
 try:
@@ -101,7 +102,7 @@ if __name__ == '__main__':
 
     wfa('/misc/hy5/scheffler/tmp/crlf', '%s\t%s\t%s\t%s\n' % (dt.now(), getuser(), gethostname(), ' '.join(sys.argv)))
     parser = argparse.ArgumentParser(
-        prog='dsc__CoReg_Sat_FourierShiftTheorem.py',
+        prog='coreg_cmd.py',
 
         description='Perform subpixel coregistration of two satellite image datasets ' \
                     'using Fourier Shift Theorem proposed by Foroosh et al. 2002: ' \
@@ -133,7 +134,7 @@ if __name__ == '__main__':
         "continue the processing."
         "The following non-standard Python libraries are required: gdal, osr, ogr, geopandas, rasterio, pykrige, "\
         "argparse and shapely. pyfftw is optional but will speed up calculation.")
-
+    # TODO update epilog
     parser.add_argument('--version', action='version', version=__version__)
 
     subparsers = parser.add_subparsers()
@@ -141,7 +142,9 @@ if __name__ == '__main__':
     # TODO add option to apply coreg results to multiple files
     ### SUBPARSER FOR COREG
     parse_coreg_global = subparsers.add_parser('global',
-        description= None,
+        description= 'Detects and corrects global X/Y shifts between a target and refernce image. Geometric shifts are '
+                     'calculated at a specific (adjustable) image position. Correction performs a global shifting in '
+                     'X- or Y direction.',
         help='detect and correct global X/Y shifts (sub argument parser)')
 
     gloArg = parse_coreg_global.add_argument
@@ -295,4 +298,5 @@ if __name__ == '__main__':
     print('\ntotal processing time: %.2fs' %(time.time()-t0))
 
 else:
-    pass # TODO
\ No newline at end of file
+    warnings.warn("The script 'coreg_cmd.py' provides a command line argument parser for CoReg_Sat and is not to be "
+                  "used as a normal Python module.")
\ No newline at end of file