coreg_cmd.py 16.5 KB
Newer Older
Daniel Scheffler's avatar
Daniel Scheffler committed
1
# -*- coding: utf-8 -*-
2
from __future__ import (division, print_function, absolute_import, unicode_literals) #unicode_literals cause GDAL not to work properly
3

Daniel Scheffler's avatar
Daniel Scheffler committed
4
5
__author__ = "Daniel Scheffler"

6
7
import time
import sys
8
import os
9
import warnings
Daniel Scheffler's avatar
Daniel Scheffler committed
10

11
# custom
Daniel Scheffler's avatar
Daniel Scheffler committed
12
13
14
15
16
17
18
19
20
21
22
23
24
try:
    import pyfftw
except ImportError:
    print('PYFFTW library is missing. However, coregistration works. But in some cases it can be much slower.')
try:
    import gdal
    import osr
    import ogr
except ImportError:
    from osgeo import gdal
    from osgeo import osr
    from osgeo import ogr

25
# internal modules
26
try:
27
    from CoReg_Sat import COREG, COREG_LOCAL, __version__
28
29
except ImportError:
    sys.path.append(os.path.abspath(os.path.join(os.path.dirname( __file__ ), '..')))
30
    from CoReg_Sat import COREG, COREG_LOCAL, __version__
31
32
33
34
35
36

try:
    import py_tools_ds
except ImportError:
    sys.path.append(os.path.abspath(os.path.dirname(__file__))) # append CoReg_Sat root directory
    import py_tools_ds
37
38


Daniel Scheffler's avatar
Daniel Scheffler committed
39

40
41
42
43
44
#sub-command functions
def run_global_coreg(args):
    COREG_obj = COREG(args.path_im0,
                      args.path_im1,
                      path_out         = args.path_out,
45
                      fmt_out          = args.fmt_out,
46
47
48
49
50
51
52
53
54
                      r_b4match        = args.br,
                      s_b4match        = args.bs,
                      wp               = args.wp,
                      ws               = args.ws,
                      max_iter         = args.max_iter,
                      max_shift        = args.max_shift,
                      align_grids      = args.align_grids,
                      match_gsd        = args.match_gsd,
                      out_gsd          = args.out_gsd,
55
56
                      data_corners_ref = args.cor0,
                      data_corners_tgt = args.cor1,
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
                      nodata           = args.nodata,
                      calc_corners     = args.calc_cor,
                      multiproc        = args.mp,
                      binary_ws        = args.bin_ws,
                      v                = args.v,
                      path_verbose_out = args.vo,
                      q                = args.q,
                      ignore_errors    = args.ignore_errors)
    COREG_obj.correct_shifts()


#sub-command functions
def run_local_coreg(args):
    CRL = COREG_LOCAL(args.path_im0,
                      args.path_im1,
                      path_out         = args.path_out,
                      fmt_out          = args.fmt_out,
                      grid_res         = args.grid_res,
                      r_b4match        = args.br,
                      s_b4match        = args.bs,
                      window_size      = args.ws,
                      max_iter         = args.max_iter,
                      max_shift        = args.max_shift,
                      #align_grids      = args.align_grids,
                      #match_gsd        = args.match_gsd,
                      #out_gsd          = args.out_gsd,
83
84
                      data_corners_ref = args.cor0,
                      data_corners_tgt = args.cor1,
85
86
87
88
89
90
91
92
93
94
95
96
                      nodata           = args.nodata,
                      calc_corners     = args.calc_cor,
                      CPUs             = None if args.mp else 1,
                      binary_ws        = args.bin_ws,
                      progress         = args.progress,
                      v                = args.v,
                      q                = args.q,
                      )
    CRL.correct_shifts()



Daniel Scheffler's avatar
Daniel Scheffler committed
97
98
if __name__ == '__main__':
    import argparse
99
100
101
102
    from socket import gethostname
    from datetime import datetime as dt
    from getpass import getuser
    from components.io import wfa
103

104
    wfa('/misc/hy5/scheffler/tmp/crlf', '%s\t%s\t%s\t%s\n' % (dt.now(), getuser(), gethostname(), ' '.join(sys.argv)))
105
    parser = argparse.ArgumentParser(
106
        prog='coreg_cmd.py',
107

Daniel Scheffler's avatar
Daniel Scheffler committed
108
        description='Perform subpixel coregistration of two satellite image datasets ' \
109
110
111
112
113
                    'using Fourier Shift Theorem proposed by Foroosh et al. 2002: ' \
                    'Foroosh, H., Zerubia, J. B., & Berthod, M. (2002). Extension of phase correlation to subpixel ' \
                    'registration. IEEE Transactions on Image Processing, 11(3), 188-199. doi:10.1109/83.988953); '\
                    'Python implementation by Daniel Scheffler (daniel.scheffler@gfz-potsdam.de)',

Daniel Scheffler's avatar
Daniel Scheffler committed
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
        epilog="DETAILED DESCRIPTION: The program detects and corrects global X/Y-shifts between two input images in "\
        "the subpixel scale, that are often present in satellite imagery. It does not correct scaling or rotation "\
        "issues and will not apply any higher grade transformation. Therefore it will also not correct for shifts "\
        "that are locally present in the input images. "
        "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."
136
        "The following non-standard Python libraries are required: gdal, osr, ogr, geopandas, rasterio, pykrige, "\
Daniel Scheffler's avatar
Daniel Scheffler committed
137
        "argparse and shapely. pyfftw is optional but will speed up calculation.")
138
    # TODO update epilog
139
140
141
    parser.add_argument('--version', action='version', version=__version__)

    subparsers = parser.add_subparsers()
142

143
144
145
    # TODO add option to apply coreg results to multiple files
    ### SUBPARSER FOR COREG
    parse_coreg_global = subparsers.add_parser('global',
146
147
148
        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.',
149
150
        help="detect and correct global X/Y shifts (sub argument parser) - "
             "use '>>> python /path/to/CoReg_Sat/coreg_cmd.py global -h' for documentation and usage hints")
151

152
153
    gloArg = parse_coreg_global.add_argument
    gloArg('path_im0', type=str, help='source path of reference image (any GDAL compatible image format is supported)')
154

155
    gloArg('path_im1', type=str, help='source path of image to be shifted (any GDAL compatible image format is supported)')
156

157
158
    gloArg('-o', nargs='?', dest='path_out', type=str, default='auto',
           help="target path of the coregistered image (default: /dir/of/im1/<im1>__shifted_to__<im0>.bsq)")
159

160
161
162
163
164
    gloArg('-fmt_out', nargs='?', type=str, help="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)", default='ENVI')



165
166
    gloArg('-br', nargs='?', type=int,
           help='band of reference image to be used for matching (starts with 1; default: 1)', default=1)
167

168
169
    gloArg('-bs', nargs='?', type=int,
           help='band of shift image to be used for matching (starts with 1; default: 1)', default=1)
170

171
172
173
    gloArg('-wp', nargs=2, metavar=('X', 'Y'), type=float,
           help="custom matching window position as map values in the same projection like the reference image "
                "(default: central position of image overlap)", default=(None,None))
174

175
176
    gloArg('-ws', nargs=2, metavar=('X size', 'Y size'), type=float,
           help="custom matching window size [pixels] (default: (512,512))", default=(512,512))
177

178
    gloArg('-max_iter', nargs='?', type=int, help="maximum number of iterations for matching (default: 5)", default=5)
179

180
181
    gloArg('-max_shift', nargs='?', type=int,
           help="maximum shift distance in reference image pixel units (default: 5 px)", default=5)
182

183
184
    gloArg('-align_grids', nargs='?', type=int, choices=[0, 1],
           help='align the coordinate grids of the output image to the reference image (default: 0)', default=0)
185

186
187
    gloArg('-match_gsd', nargs='?', type=int, choices=[0, 1],
           help='match the output pixel size to the pixel size of the reference image (default: 0)', default=0)
188

189
190
    gloArg('-out_gsd', nargs=2, type=float, help='xgsd ygsd: set the output pixel size in map units'\
           '(default: original pixel size of the image to be shifted)', metavar=('xgsd','ygsd'))
191

192
193
    # TODO implement footprint_poly_ref, footprint_poly_tgt

194
195
    gloArg('-cor0', nargs=8, type=float, help="map coordinates of data corners within reference image: ",
           metavar=tuple("UL-X UL-Y UR-X UR-Y LR-X LR-Y LL-X LL-Y".split(' ')), default=None)
196

197
198
    gloArg('-cor1', nargs=8, type=float, help="map coordinates of data corners within image to be shifted: ",
           metavar=tuple("UL-X UL-Y UR-X UR-Y LR-X LR-Y LL-X LL-Y".split(' ')), default=None)
199

200
201
    gloArg('-nodata', nargs=2, type=float, metavar=('im0', 'im1'),
           help='no data values for reference image and image to be shifted', default=(None,None))
202

203
204
205
    gloArg('-calc_cor', nargs='?', type=int, choices=[0, 1], default=1,
           help="calculate true positions of the dataset corners in order to get a useful matching window position "
                "within the actual image overlap (default: 1; deactivated if '-cor0' and '-cor1' are given")
206

207
    gloArg('-mp', nargs='?', type=int, help='enable multiprocessing (default: 1)', default=1, choices=[0, 1])
208

209
210
    gloArg('-bin_ws', nargs='?', type=int,
           help='use binary X/Y dimensions for the matching window (default: 1)', default=1, choices=[0, 1])
211

212
213
    gloArg('-quadratic_win', nargs='?', type=int,
           help='force a quadratic matching window (default: 1)', default=1, choices=[0, 1])
214

215
216
217
218
219
220
221
222
223
224
225
    gloArg('-v', nargs='?', type=int, help='verbose mode (default: 0)', default=0, choices=[0, 1])

    gloArg('-vo', nargs='?', type=str, choices=[0, 1], help='an optional output directory for outputs of verbose mode'
           '(if not given, no outputs are written to disk)', default=0, )

    gloArg('-q', nargs='?', type=int, help='quiet mode (default: 0)', default=0, choices=[0, 1])

    gloArg('-ignore_errors', nargs='?', type=int, help='Useful for batch processing. (default: 0) In case of error '
           'COREG.success == False and COREG.x_shift_px/COREG.y_shift_px is None', default=0, choices=[0,1])

    parse_coreg_global.set_defaults(func=run_global_coreg)
226
227
228



229
230
231
232
233
234
235
    ### SUBPARSER FOR COREG LOCAL
    parse_coreg_local = subparsers.add_parser('local',
        description= 'Applies the algorithm to detect spatial shifts to the whole overlap area of the input images. '
                     'Spatial shifts are calculated for each point in grid of which the parameters can be adjusted '
                     'using keyword arguments. Shift correction performs a polynomial transformation using the '
                     'calculated shifts of each point in the grid as GCPs. Thus this class can be used to correct ' \
                     'for locally varying geometric distortions of the target image.',
236
237
        help="detect and correct local shifts (sub argument parser)"
             "use '>>> python /path/to/CoReg_Sat/coreg_cmd.py local -h' for documentation and usage hints")
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266

    locArg = parse_coreg_local.add_argument
    locArg('path_im0', type=str, help='source path of reference image (any GDAL compatible image format is supported)')

    locArg('path_im1', type=str, help='source path of image to be shifted (any GDAL compatible image format is supported)')

    locArg('grid_res', type=int, help='quality grid resolution in pixels of the target image')

    locArg('-o', nargs='?', type=str, dest='path_out', default='auto',
           help="target path of the coregistered image. If 'auto' (default): /dir/of/im1/<im1>__shifted_to__<im0>.bsq")

    locArg('-fmt_out', nargs='?', type=str, help="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)", default='ENVI')

    locArg('-projectDir', nargs='?', type=str, help=None, default=None)

    locArg('-ws', nargs=2, type=int, help='custom matching window size [pixels] (default: (256,256))')

    locArg('-br', nargs='?', type=int,
           help='band of reference image to be used for matching (starts with 1; default: 1)', default=1)

    locArg('-bs', nargs='?', type=int,
           help='band of shift image to be used for matching (starts with 1; default: 1)', default=1)

    locArg('-max_iter', nargs='?', type=int, help="maximum number of iterations for matching (default: 5)", default=5)

    locArg('-max_shift', nargs='?', type=int,
           help="maximum shift distance in reference image pixel units (default: 5 px)", default=5)

267
268
    # TODO implement footprint_poly_ref, footprint_poly_tgt

269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
    locArg('-cor0', nargs=8, type=float, help="map coordinates of data corners within reference image: ",
           metavar=tuple("UL-X UL-Y UR-X UR-Y LR-X LR-Y LL-X LL-Y".split(' ')), default=None)

    locArg('-cor1', nargs=8, type=float, help="map coordinates of data corners within image to be shifted: ",
           metavar=tuple("UL-X UL-Y UR-X UR-Y LR-X LR-Y LL-X LL-Y".split(' ')), default=None)

    locArg('-nodata', nargs=2, type=float, metavar=('im0', 'im1'),
           help='no data values for reference image and image to be shifted', default=(None,None))

    locArg('-calc_cor', nargs='?', type=int, choices=[0, 1], default=1,
           help="calculate true positions of the dataset corners in order to get a useful matching window position "
                "within the actual image overlap (default: 1; deactivated if '-cor0' and '-cor1' are given")

    locArg('-mp', nargs='?', type=int, help='enable multiprocessing (default: 1)', default=1, choices=[0, 1])

    locArg('-bin_ws', nargs='?', type=int,
           help='use binary X/Y dimensions for the matching window (default: 1)', default=1, choices=[0, 1])

    locArg('-quadratic_win', nargs='?', type=int,
           help='force a quadratic matching window (default: 1)', default=1, choices=[0, 1])

    locArg('-progress', nargs='?', type=int, help='show progress bars (default: 1)', default=1, choices=[0, 1])

    locArg('-v', nargs='?', type=int, help='verbose mode (default: 0)', default=0, choices=[0, 1])

    locArg('-q', nargs='?', type=int, help='quiet mode (default: 0)', default=0, choices=[0, 1])

    parse_coreg_local.set_defaults(func=run_local_coreg)



    parsed_args = parser.parse_args()
Daniel Scheffler's avatar
Daniel Scheffler committed
301
302

    print('==================================================================\n'
303
          '#                     CoReg_Sat v%s                   #'%__version__+'\n'
Daniel Scheffler's avatar
Daniel Scheffler committed
304
305
306
307
308
309
          '#          SUBPIXEL COREGISTRATION FOR SATELLITE IMAGERY         #\n'
          '#          - algorithm proposed by Foroosh et al. 2002           #\n'
          '#          - python implementation by Daniel Scheffler           #\n'
          '==================================================================\n')

    t0 = time.time()
310
    parsed_args.func(parsed_args)
311
312
313
    print('\ntotal processing time: %.2fs' %(time.time()-t0))

else:
314
315
    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.")