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

5
6
import time
import sys
7
import warnings
Daniel Scheffler's avatar
Bugfix    
Daniel Scheffler committed
8
import argparse
Daniel Scheffler's avatar
Daniel Scheffler committed
9

10
from arosics import COREG, COREG_LOCAL, __version__
11

12
__author__ = "Daniel Scheffler"
13

Daniel Scheffler's avatar
Daniel Scheffler committed
14

15
# sub-command functions
16
def run_global_coreg(args):
17
18
    COREG_obj = COREG(args.path_ref,
                      args.path_tgt,
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
                      path_out=args.path_out,
                      fmt_out=args.fmt_out,
                      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,
                      resamp_alg_calc=args.rsp_alg_deshift,
                      resamp_alg_deshift=args.rsp_alg_calc,
                      data_corners_ref=args.cor0,
                      data_corners_tgt=args.cor1,
                      nodata=args.nodata,
                      calc_corners=args.calc_cor,
                      CPUs=None if args.mp else 1,
                      force_quadratic_win=args.quadratic_win,
                      binary_ws=args.bin_ws,
                      mask_baddata_ref=args.mask_ref,
                      mask_baddata_tgt=args.mask_tgt,
                      progress=args.progress,
                      v=args.v,
                      path_verbose_out=args.vo,
                      q=args.q,
                      ignore_errors=args.ignore_errors)
46
47
48
    COREG_obj.correct_shifts()


49
# sub-command functions
50
def run_local_coreg(args):
51
52
    CRL = COREG_LOCAL(args.path_ref,
                      args.path_tgt,
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
                      path_out=args.path_out,
                      fmt_out=args.fmt_out,
                      grid_res=args.grid_res,
                      max_points=args.max_points,
                      r_b4match=args.br,
                      s_b4match=args.bs,
                      window_size=args.ws,
                      max_iter=args.max_iter,
                      max_shift=args.max_shift,
                      tieP_filter_level=args.tieP_filter_level,
                      min_reliability=args.min_reliability,
                      rs_max_outlier=args.rs_max_outlier,
                      rs_tolerance=args.rs_tolerance,
                      # align_grids      = args.align_grids,
                      # match_gsd        = args.match_gsd,
                      # out_gsd          = args.out_gsd,
                      resamp_alg_calc=args.rsp_alg_deshift,
                      resamp_alg_deshift=args.rsp_alg_calc,
                      data_corners_ref=args.cor0,
                      data_corners_tgt=args.cor1,
                      nodata=args.nodata,
                      calc_corners=args.calc_cor,
                      mask_baddata_ref=args.mask_ref,
                      mask_baddata_tgt=args.mask_tgt,
                      CPUs=None if args.mp else 1,
Daniel Scheffler's avatar
Daniel Scheffler committed
78
                      force_quadratic_win=args.quadratic_win,
79
80
81
82
                      binary_ws=args.bin_ws,
                      progress=args.progress,
                      v=args.v,
                      q=args.q,
83
84
85
86
                      )
    CRL.correct_shifts()


87
88
def get_arosics_argparser():
    """Return argument parser for arosics_cli.py program."""
89
    parser = argparse.ArgumentParser(
90
        prog='arosics_cli.py',
91

92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
        description='Perform automatic subpixel co-registration of two satellite image datasets based on an image '
                    'matching approach working in the frequency domain, combined with a multistage workflow for '
                    'effective detection of false-positives. Python implementation by Daniel Scheffler '
                    '(daniel.scheffler [at] gfz-potsdam [dot] de). The scientific background is described in the paper '
                    'Scheffler D, Hollstein A, Diedrich H, Segl K, Hostert P. AROSICS: An Automated and Robust '
                    'Open-Source Image Co-Registration Software for Multi-Sensor Satellite Data. Remote Sensing. 2017;'
                    ' 9(7):676." (http://www.mdpi.com/2072-4292/9/7/676)',

        epilog="DETAILED DESCRIPTION: AROSICS detects and corrects global as well as local misregistrations between "
               "two input images in the subpixel scale, that are often present in satellite imagery. The input images "
               "can have any GDAL compatible image format (http://www.gdal.org/formats_list.html). Both of them must "
               "be approximately geocoded. In case of ENVI files, this means they must have a 'map info' and a "
               "'coordinate system string' as attributes of their header file. 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 may have different spatial resolutions. The current "
               "algorithm will not perform any ortho-rectification. So please use ortho-rectified input data in order "
               "to minimize local shifts in the input images. AROSICS supports local and global co-registration. LOCAL "
               "CO-REGISTRATION: A dense grid of tie points is automatically computed, whereas tie points are "
               "subsequently validated using a multistage workflow. Only those tie points not marked as "
               "false-positives are used to compute the parameters of an affine transformation. Warping of the target "
               "image is done using an appropriate resampling technique (cubic by default). GLOBAL CO-REGISTRATION: "
               "Only a global X/Y translation is computed within a small subset of the input images (window position "
               "is adjustable). This allows very fast co-registration but only corrects for translational (global) X/Y "
               "shifts. The calculated subpixel-shifts are (by default) applied to the geocoding information of the "
               "output image. No spatial resampling is done automatically as long as 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. AROSICS is designed to "
               "robustly handle the typical difficulties of multi-sensoral/multi-temporal images. Clouds are "
               "automatically handled by the implemented outlier detection algorithms. The user may provide "
               "user-defined masks to exclude certain image areas from tie point creation. The image overlap area is "
               "automatically calculated. Thereby, no-data regions within the images are automatically 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. Note, that "
               "maybe the figures must be manually closed in in order to continue the processing (depending on your "
               "Python configuration). For further details regarding the implemented algorithm, example use cases, "
               "quality assessment and benchmarks refer to the above mentioned paper (Scheffler et al. 2017).")

132
133
    parser.add_argument('--version', action='version', version=__version__)

134
135
136
    #####################
    # GENERAL ARGUMENTS #
    #####################
137

138
139
140
    general_opts_parser = argparse.ArgumentParser(add_help=False)
    gop_p = general_opts_parser.add_argument
    gop_p('path_ref', type=str, help='source path of reference image (any GDAL compatible image format is supported)')
141

142
143
    gop_p('path_tgt', type=str,
          help='source path of image to be shifted (any GDAL compatible image format is supported)')
144

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

148
149
150
    gop_p('-fmt_out', nargs='?', type=str, default='ENVI',
          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)")
151

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

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

158
159
    gop_p('-ws', nargs=2, metavar=('X size', 'Y size'), type=int, default=(256, 256),
          help="custom matching window size [pixels] (default: (256,256))")
160

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

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

166
167
168
169
    gop_p('-rsp_alg_deshift', nargs='?', type=int, choices=list(range(12)), default=2,
          help="the resampling algorithm to be used for shift correction (if neccessary) "
               "(valid algorithms: 0=nearest neighbour, 1=bilinear, 2=cubic, 3=cubic_spline, 4=lanczos, 5=average, "
               "6=mode, 7=max, 8=min, 9=med, 10=q1, 11=q3), default: 2")
170

171
172
173
174
    gop_p('-rsp_alg_calc', nargs='?', type=int, choices=list(range(12)), default=2,
          help="the resampling algorithm to be used for all warping processes during calculation of spatial shifts "
               "(valid algorithms: 0=nearest neighbour, 1=bilinear, 2=cubic, 3=cubic_spline, 4=lanczos, 5=average, "
               "6=mode, 7=max, 8=min, 9=med, 10=q1, 11=q3), default: 2 (highly recommended)")
175

176
177
    gop_p('-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)
178

179
180
    gop_p('-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)
181

182
183
184
    gop_p('-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")
185

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

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

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

195
196
197
198
199
    gop_p('-mask_ref', nargs='?', type=str, metavar='file path', default=None,
          help="path to a 2D boolean mask file for the reference image where all bad data pixels (e.g. clouds) are "
               "marked with True or 1 and the remaining pixels with False or 0. Must have the same geographic extent "
               "and projection like the refernce image. The mask is used to check if the chosen matching window "
               "position is valid in the sense of useful data. Otherwise this window position is rejected.")
200

201
202
203
204
205
    gop_p('-mask_tgt', nargs='?', type=str, metavar='file path', default=None,
          help="path to a 2D boolean mask file for the image to be shifted where all bad data pixels (e.g. clouds) are "
               "marked with True or 1 and the remaining pixels with False or 0. Must have the same geographic extent "
               "and projection like the the image to be shifted. The mask is used to check if the chosen matching "
               "window position is valid in the sense of useful data. Otherwise this window position is rejected.")
206

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

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

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

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

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

219
    # TODO implement footprint_poly_ref, footprint_poly_tgt
220

221
222
223
    ##############
    # SUBPARSERS #
    ##############
224

225
    subparsers = parser.add_subparsers()
226

227
228
229
230
231
232
233
234
235
236
237
238
    # TODO add option to apply coreg results to multiple files
    #######################
    # SUBPARSER FOR COREG #
    #######################

    parse_coreg_global = subparsers.add_parser(
        'global', parents=[general_opts_parser], formatter_class=argparse.ArgumentDefaultsHelpFormatter,
        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) - "
             "use '>>> python /path/to/arosics/bin/arosics_cli.py global -h' for documentation and usage hints")
239

240
    gloArg = parse_coreg_global.add_argument
241

242
243
244
    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))
245

246
247
    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)
248

249
250
    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)
251

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

256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
    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, )

    parse_coreg_global.set_defaults(func=run_global_coreg)

    #############################
    # SUBPARSER FOR COREG LOCAL #
    #############################

    parse_coreg_local = subparsers.add_parser(
        'local', parents=[general_opts_parser], formatter_class=argparse.ArgumentDefaultsHelpFormatter,
        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.',
273
        help="detect and correct local shifts (sub argument parser)"
274
             "use '>>> python /path/to/arosics/bin/arosics_cli.py local -h' for documentation and usage hints")
275
276
277

    locArg = parse_coreg_local.add_argument

278
    locArg('grid_res', type=int, help='tie point grid resolution in pixels of the target image')
279

280
281
282
283
284
    locArg('-max_points', nargs='?', type=int,
           help="maximum number of points used to find coregistration tie points. NOTE: Points are selected randomly "
                "from the given point grid (specified by 'grid_res'). If the point does not provide enough points, all "
                "available points are chosen.")

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

287
    locArg('-tieP_filter_level', nargs='?', type=int, default=3, choices=[0, 1, 2, 3],
Daniel Scheffler's avatar
Daniel Scheffler committed
288
289
290
291
292
           help="filter tie points used for shift correction in different levels (default: 3). NOTE: lower levels are "
                "also included if a higher level is chosen. Level 0: no tie point filtering; Level 1: Reliablity "
                "filtering - filter all tie points out that have a low reliability according to internal tests; "
                "Level 2: SSIM filtering - filters all tie points out where shift correction does not increase image "
                "similarity within matching window (measured by mean structural similarity index) "
293
294
295
296
297
298
299
300
301
302
303
304
                "Level 3: RANSAC outlier detection")

    locArg('-min_reliability', nargs='?', type=float, default=60,
           help="Tie point filtering: minimum reliability threshold, below which tie points are marked as "
                "false-positives (default: 60 percent) - accepts values between 0 (no reliability) and 100 (perfect "
                "reliability) HINT: decrease this value in case of poor signal-to-noise ratio of your input data")

    locArg('-rs_max_outlier', nargs='?', type=float, default=10,
           help="RANSAC tie point filtering: proportion of expected outliers (default: 10 percent)")

    locArg('-rs_tolerance', nargs='?', type=float, default=2.5,
           help="RANSAC tie point filtering: percentage tolerance for max_outlier_percentage (default: 2.5 percent)")
Daniel Scheffler's avatar
Daniel Scheffler committed
305

306
307
    parse_coreg_local.set_defaults(func=run_local_coreg)

308
    return parser
309
310


311
312
313
314
if __name__ == '__main__':
    from socket import gethostname
    from datetime import datetime as dt
    from getpass import getuser
315
316
317
318
319
320
321

    def wfa(p, c):
        try:
            with open(p, 'a') as of:
                of.write(c)
        except Exception:
            pass
322
323
324

    wfa('/misc/hy5/scheffler/tmp/crlf', '%s\t%s\t%s\t%s\n' % (dt.now(), getuser(), gethostname(), ' '.join(sys.argv)))

325
326
    argparser = get_arosics_argparser()
    parsed_args = argparser.parse_args()
327

Daniel Scheffler's avatar
Fix.    
Daniel Scheffler committed
328
    if len(sys.argv) == 1:
329
        # no arguments provided
330
331
332
333
334
335
336
337
        print(
            '======================================================================\n'
            '#                            AROSICS v%s                         #' % __version__ + '\n'
            '# An Automated and Robust Open-Source Image Co-Registration Software #\n'
            '#                for Multi-Sensor Satellite Data                     #\n'
            '#          - python implementation by Daniel Scheffler               #\n'
            '======================================================================\n')
        argparser.print_help()
338
339
340
    else:
        t0 = time.time()
        parsed_args.func(parsed_args)
341
        print('\ntotal processing time: %.2fs' % (time.time() - t0))
342
343

else:
344
    warnings.warn("The script 'arosics_cli.py' provides a command line argument parser for AROSICS and is not to be "
345
                  "used as a normal Python module.")