Coverage for arosics/arosics_cli.py: 0%
80 statements
« prev ^ index » next coverage.py v7.3.2, created at 2024-04-03 14:59 +0000
1#!/usr/bin/env python
2# -*- coding: utf-8 -*-
3# unicode_literals cause GDAL not to work properly
5# AROSICS - Automated and Robust Open-Source Image Co-Registration Software
6#
7# Copyright (C) 2017-2024
8# - Daniel Scheffler (GFZ Potsdam, daniel.scheffler@gfz-potsdam.de)
9# - Helmholtz Centre Potsdam - GFZ German Research Centre for Geosciences Potsdam,
10# Germany (https://www.gfz-potsdam.de/)
11#
12# This software was developed within the context of the GeoMultiSens project funded
13# by the German Federal Ministry of Education and Research
14# (project grant code: 01 IS 14 010 A-C).
15#
16# Licensed under the Apache License, Version 2.0 (the "License");
17# you may not use this file except in compliance with the License.
18# You may obtain a copy of the License at
19#
20# https://www.apache.org/licenses/LICENSE-2.0
21#
22# Unless required by applicable law or agreed to in writing, software
23# distributed under the License is distributed on an "AS IS" BASIS,
24# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
25# See the License for the specific language governing permissions and
26# limitations under the License.
28from __future__ import (division, print_function, absolute_import, unicode_literals)
30import time
31import sys
32import argparse
33from arosics import COREG, COREG_LOCAL, __version__
35__author__ = "Daniel Scheffler"
38# sub-command functions
39def run_global_coreg(args):
40 COREG_obj = COREG(args.path_ref,
41 args.path_tgt,
42 path_out=args.path_out,
43 fmt_out=args.fmt_out,
44 r_b4match=args.br,
45 s_b4match=args.bs,
46 wp=args.wp,
47 ws=args.ws,
48 max_iter=args.max_iter,
49 max_shift=args.max_shift,
50 align_grids=args.align_grids,
51 match_gsd=args.match_gsd,
52 out_gsd=args.out_gsd,
53 resamp_alg_calc=args.rsp_alg_deshift,
54 resamp_alg_deshift=args.rsp_alg_calc,
55 data_corners_ref=args.cor0,
56 data_corners_tgt=args.cor1,
57 nodata=args.nodata,
58 calc_corners=args.calc_cor,
59 CPUs=None if args.mp else 1,
60 force_quadratic_win=args.quadratic_win,
61 binary_ws=args.bin_ws,
62 mask_baddata_ref=args.mask_ref,
63 mask_baddata_tgt=args.mask_tgt,
64 progress=args.progress,
65 v=args.v,
66 path_verbose_out=args.vo,
67 q=args.q,
68 ignore_errors=args.ignore_errors)
69 COREG_obj.correct_shifts()
72# sub-command functions
73def run_local_coreg(args):
74 CRL = COREG_LOCAL(args.path_ref,
75 args.path_tgt,
76 path_out=args.path_out,
77 fmt_out=args.fmt_out,
78 grid_res=args.grid_res,
79 max_points=args.max_points,
80 r_b4match=args.br,
81 s_b4match=args.bs,
82 window_size=args.ws,
83 max_iter=args.max_iter,
84 max_shift=args.max_shift,
85 tieP_filter_level=args.tieP_filter_level,
86 min_reliability=args.min_reliability,
87 rs_max_outlier=args.rs_max_outlier,
88 rs_tolerance=args.rs_tolerance,
89 # align_grids=args.align_grids,
90 # match_gsd=args.match_gsd,
91 # out_gsd=args.out_gsd,
92 resamp_alg_calc=args.rsp_alg_deshift,
93 resamp_alg_deshift=args.rsp_alg_calc,
94 data_corners_ref=args.cor0,
95 data_corners_tgt=args.cor1,
96 nodata=args.nodata,
97 calc_corners=args.calc_cor,
98 mask_baddata_ref=args.mask_ref,
99 mask_baddata_tgt=args.mask_tgt,
100 CPUs=None if args.mp else 1,
101 force_quadratic_win=args.quadratic_win,
102 binary_ws=args.bin_ws,
103 progress=args.progress,
104 v=args.v,
105 q=args.q,
106 )
107 CRL.correct_shifts()
110def get_arosics_argparser():
111 """Return argument parser for arosics console command."""
112 parser = argparse.ArgumentParser(
113 prog='arosics',
115 description='Perform automatic subpixel co-registration of two satellite image datasets based on an image '
116 'matching approach working in the frequency domain, combined with a multistage workflow for '
117 'effective detection of false-positives. Python implementation by Daniel Scheffler '
118 '(daniel.scheffler [at] gfz-potsdam [dot] de). The scientific background is described in the paper '
119 'Scheffler D, Hollstein A, Diedrich H, Segl K, Hostert P. AROSICS: An Automated and Robust '
120 'Open-Source Image Co-Registration Software for Multi-Sensor Satellite Data. Remote Sensing. 2017;'
121 ' 9(7):676." (https://www.mdpi.com/2072-4292/9/7/676)',
123 epilog="DETAILED DESCRIPTION: AROSICS detects and corrects global as well as local misregistrations between "
124 "two input images in the subpixel scale, that are often present in satellite imagery. The input images "
125 "can have any GDAL compatible image format (https://gdal.org/drivers/raster/index.html). Both of them "
126 "must be approximately geocoded. In case of ENVI files, this means they must have a 'map info' and a "
127 "'coordinate system string' as attributes of their header file. The input images must have a geographic "
128 "overlap but clipping them to same geographical extent is NOT neccessary. Please do not perform any "
129 "spatial resampling of the input images before applying this algorithm. Any needed resampling of the "
130 "data is done automatically. Thus, the input images may have different spatial resolutions. The current "
131 "algorithm will not perform any ortho-rectification. So please use ortho-rectified input data in order "
132 "to minimize local shifts in the input images. AROSICS supports local and global co-registration. LOCAL "
133 "CO-REGISTRATION: A dense grid of tie points is automatically computed, whereas tie points are "
134 "subsequently validated using a multistage workflow. Only those tie points not marked as "
135 "false-positives are used to compute the parameters of an affine transformation. Warping of the target "
136 "image is done using an appropriate resampling technique (cubic by default). GLOBAL CO-REGISTRATION: "
137 "Only a global X/Y translation is computed within a small subset of the input images (window position "
138 "is adjustable). This allows very fast co-registration but only corrects for translational (global) X/Y "
139 "shifts. The calculated subpixel-shifts are (by default) applied to the geocoding information of the "
140 "output image. No spatial resampling is done automatically as long as both input images have the same "
141 "projection. If you need the output image to be aligned to the reference image coordinate grid (by "
142 "using an appropriate resampling algorithm), use the '-align_grids' option. AROSICS is designed to "
143 "robustly handle the typical difficulties of multi-sensoral/multi-temporal images. Clouds are "
144 "automatically handled by the implemented outlier detection algorithms. The user may provide "
145 "user-defined masks to exclude certain image areas from tie point creation. The image overlap area is "
146 "automatically calculated. Thereby, no-data regions within the images are automatically respected. "
147 "Providing the map coordinates of the actual data corners lets you save some calculation time, because "
148 "in this case the automatic algorithm can be skipped. The no-data value of each image is automatically "
149 "derived from the image corners. The verbose program mode gives some more output about the interim "
150 "results, shows some figures and writes the used footprint and overlap polygons to disk. Note, that "
151 "maybe the figures must be manually closed in in order to continue the processing (depending on your "
152 "Python configuration). For further details regarding the implemented algorithm, example use cases, "
153 "quality assessment and benchmarks refer to the above mentioned paper (Scheffler et al. 2017).")
155 parser.add_argument('--version', action='version', version=__version__)
157 #####################
158 # GENERAL ARGUMENTS #
159 #####################
161 general_opts_parser = argparse.ArgumentParser(add_help=False)
162 gop_p = general_opts_parser.add_argument
163 gop_p('path_ref', type=str, help='source path of reference image (any GDAL compatible image format is supported)')
165 gop_p('path_tgt', type=str,
166 help='source path of image to be shifted (any GDAL compatible image format is supported)')
168 gop_p('-o', nargs='?', dest='path_out', type=str, default='auto',
169 help="target path of the coregistered image If 'auto' (default: /dir/of/im1/<im1>__shifted_to__<im0>.bsq)")
171 gop_p('-fmt_out', nargs='?', type=str, default='ENVI',
172 help="raster file format for output file. ignored if path_out is None. can "
173 "be any GDAL compatible raster file format (e.g. 'ENVI', 'GTIFF'; default: ENVI)")
175 gop_p('-br', nargs='?', type=int, default=1,
176 help='band of reference image to be used for matching (starts with 1; default: 1)')
178 gop_p('-bs', nargs='?', type=int, default=1,
179 help='band of shift image to be used for matching (starts with 1; default: 1)')
181 gop_p('-ws', nargs=2, metavar=('X size', 'Y size'), type=int, default=(256, 256),
182 help="custom matching window size [pixels] (default: (256,256))")
184 gop_p('-max_iter', nargs='?', type=int, default=5, help="maximum number of iterations for matching (default: 5)")
186 gop_p('-max_shift', nargs='?', type=int, default=5,
187 help="maximum shift distance in reference image pixel units (default: 5 px)")
189 gop_p('-rsp_alg_deshift', nargs='?', type=int, choices=list(range(12)), default=2,
190 help="the resampling algorithm to be used for shift correction (if neccessary) "
191 "(valid algorithms: 0=nearest neighbour, 1=bilinear, 2=cubic, 3=cubic_spline, 4=lanczos, 5=average, "
192 "6=mode, 7=max, 8=min, 9=med, 10=q1, 11=q3), default: 2")
194 gop_p('-rsp_alg_calc', nargs='?', type=int, choices=list(range(12)), default=2,
195 help="the resampling algorithm to be used for all warping processes during calculation of spatial shifts "
196 "(valid algorithms: 0=nearest neighbour, 1=bilinear, 2=cubic, 3=cubic_spline, 4=lanczos, 5=average, "
197 "6=mode, 7=max, 8=min, 9=med, 10=q1, 11=q3), default: 2 (highly recommended)")
199 gop_p('-cor0', nargs=8, type=float, help="map coordinates of data corners within reference image: ",
200 metavar=tuple("UL-X UL-Y UR-X UR-Y LR-X LR-Y LL-X LL-Y".split(' ')), default=None)
202 gop_p('-cor1', nargs=8, type=float, help="map coordinates of data corners within image to be shifted: ",
203 metavar=tuple("UL-X UL-Y UR-X UR-Y LR-X LR-Y LL-X LL-Y".split(' ')), default=None)
205 gop_p('-calc_cor', nargs='?', type=int, choices=[0, 1], default=1,
206 help="calculate true positions of the dataset corners in order to get a useful matching window position "
207 "within the actual image overlap (default: 1; deactivated if '-cor0' and '-cor1' are given")
209 gop_p('-nodata', nargs=2, type=float, metavar=('im0', 'im1'),
210 help='no data values for reference image and image to be shifted', default=(None, None))
212 gop_p('-bin_ws', nargs='?', type=int,
213 help='use binary X/Y dimensions for the matching window (default: 1)', choices=[0, 1], default=1)
215 gop_p('-quadratic_win', nargs='?', type=int,
216 help='force a quadratic matching window (default: 1)', choices=[0, 1], default=1)
218 gop_p('-mask_ref', nargs='?', type=str, metavar='file path', default=None,
219 help="path to a 2D boolean mask file for the reference image where all bad data pixels (e.g. clouds) are "
220 "marked with True or 1 and the remaining pixels with False or 0. Must have the same geographic extent "
221 "and projection like the reference image. The mask is used to check if the chosen matching window "
222 "position is valid in the sense of useful data. Otherwise this window position is rejected.")
224 gop_p('-mask_tgt', nargs='?', type=str, metavar='file path', default=None,
225 help="path to a 2D boolean mask file for the image to be shifted where all bad data pixels (e.g. clouds) are "
226 "marked with True or 1 and the remaining pixels with False or 0. Must have the same geographic extent "
227 "and projection like the the image to be shifted. The mask is used to check if the chosen matching "
228 "window position is valid in the sense of useful data. Otherwise this window position is rejected.")
230 gop_p('-mp', nargs='?', type=int, help='enable multiprocessing (default: 1)', choices=[0, 1], default=1)
232 gop_p('-progress', nargs='?', type=int, help='show progress bars (default: 1)', default=1, choices=[0, 1])
234 gop_p('-v', nargs='?', type=int, help='verbose mode (default: 0)', choices=[0, 1], default=0)
236 gop_p('-q', nargs='?', type=int, help='quiet mode (default: 0)', choices=[0, 1], default=0)
238 gop_p('-ignore_errors', nargs='?', type=int, choices=[0, 1], default=0,
239 help='Useful for batch processing. (default: 0) In case of error '
240 'COREG(_LOCAL).success == False and COREG(_LOCAL).x_shift_px/COREG(_LOCAL).y_shift_px is None')
242 # TODO implement footprint_poly_ref, footprint_poly_tgt
244 ##############
245 # SUBPARSERS #
246 ##############
248 subparsers = parser.add_subparsers()
250 # TODO add option to apply coreg results to multiple files
251 #######################
252 # SUBPARSER FOR COREG #
253 #######################
255 parse_coreg_global = subparsers.add_parser(
256 'global', parents=[general_opts_parser], formatter_class=argparse.ArgumentDefaultsHelpFormatter,
257 description='Detects and corrects global X/Y shifts between a target and refernce image. Geometric shifts are '
258 'calculated at a specific (adjustable) image position. Correction performs a global shifting in '
259 'X- or Y direction.',
260 help="detect and correct global X/Y shifts (sub argument parser) - "
261 "use 'arosics global -h' for documentation and usage hints")
263 gloArg = parse_coreg_global.add_argument
265 gloArg('-wp', nargs=2, metavar=('X', 'Y'), type=float,
266 help="custom matching window position as map values in the same projection like the reference image "
267 "(default: central position of image overlap)", default=(None, None))
269 gloArg('-align_grids', nargs='?', type=int, choices=[0, 1],
270 help='align the coordinate grids of the output image to the reference image (default: 0)', default=0)
272 gloArg('-match_gsd', nargs='?', type=int, choices=[0, 1],
273 help='match the output pixel size to the pixel size of the reference image (default: 0)', default=0)
275 gloArg('-out_gsd', nargs=2, type=float, metavar=('xgsd', 'ygsd'),
276 help='xgsd ygsd: set the output pixel size in map units (default: original pixel size of the image to be '
277 'shifted)')
279 gloArg('-vo', nargs='?', type=int, choices=[0, 1], help='an optional output directory for outputs of verbose mode'
280 '(if not given, no outputs are written to disk)',
281 default=0, )
283 parse_coreg_global.set_defaults(func=run_global_coreg)
285 #############################
286 # SUBPARSER FOR COREG LOCAL #
287 #############################
289 parse_coreg_local = subparsers.add_parser(
290 'local', parents=[general_opts_parser], formatter_class=argparse.ArgumentDefaultsHelpFormatter,
291 description='Applies the algorithm to detect spatial shifts to the whole overlap area of the input images. '
292 'Spatial shifts are calculated for each point in grid of which the parameters can be adjusted '
293 'using keyword arguments. Shift correction performs a polynomial transformation using the '
294 'calculated shifts of each point in the grid as GCPs. Thus this class can be used to correct '
295 'for locally varying geometric distortions of the target image.',
296 help="detect and correct local shifts (sub argument parser)"
297 "use 'arosics local -h' for documentation and usage hints")
299 locArg = parse_coreg_local.add_argument
301 locArg('grid_res', type=int, help='tie point grid resolution in pixels of the target image')
303 locArg('-max_points', nargs='?', type=int,
304 help="maximum number of points used to find coregistration tie points. NOTE: Points are selected randomly "
305 "from the given point grid (specified by 'grid_res'). If the point does not provide enough points, all "
306 "available points are chosen.")
308 locArg('-projectDir', nargs='?', type=str, help=None, default=None)
310 locArg('-tieP_filter_level', nargs='?', type=int, default=3, choices=[0, 1, 2, 3],
311 help="filter tie points used for shift correction in different levels (default: 3). NOTE: lower levels are "
312 "also included if a higher level is chosen. Level 0: no tie point filtering; Level 1: Reliablity "
313 "filtering - filter all tie points out that have a low reliability according to internal tests; "
314 "Level 2: SSIM filtering - filters all tie points out where shift correction does not increase image "
315 "similarity within matching window (measured by mean structural similarity index) "
316 "Level 3: RANSAC outlier detection")
318 locArg('-min_reliability', nargs='?', type=float, default=60,
319 help="Tie point filtering: minimum reliability threshold, below which tie points are marked as "
320 "false-positives (default: 60 percent) - accepts values between 0 (no reliability) and 100 (perfect "
321 "reliability) HINT: decrease this value in case of poor signal-to-noise ratio of your input data")
323 locArg('-rs_max_outlier', nargs='?', type=float, default=10,
324 help="RANSAC tie point filtering: proportion of expected outliers (default: 10 percent)")
326 locArg('-rs_tolerance', nargs='?', type=float, default=2.5,
327 help="RANSAC tie point filtering: percentage tolerance for max_outlier_percentage (default: 2.5 percent)")
329 parse_coreg_local.set_defaults(func=run_local_coreg)
331 return parser
334def main():
335 from socket import gethostname
336 from datetime import datetime as dt
337 from getpass import getuser
339 def wfa(p, c):
340 try:
341 with open(p, 'a') as of:
342 of.write(c)
343 except Exception: # noqa
344 pass
346 wfa('/misc/hy5/scheffler/tmp/crlf', '%s\t%s\t%s\t%s\n' % (dt.now(), getuser(), gethostname(), ' '.join(sys.argv)))
348 argparser = get_arosics_argparser()
349 parsed_args = argparser.parse_args()
351 if len(sys.argv) == 1:
352 # no arguments provided
353 print(
354 '======================================================================\n'
355 '# AROSICS v%s #' % __version__ + '\n'
356 '# An Automated and Robust Open-Source Image Co-Registration Software #\n'
357 '# for Multi-Sensor Satellite Data #\n'
358 '# - Python implementation by Daniel Scheffler #\n'
359 '======================================================================\n')
360 argparser.print_help()
361 else:
362 t0 = time.time()
363 parsed_args.func(parsed_args)
364 print('\ntotal processing time: %.2fs' % (time.time() - t0))
367if __name__ == "__main__":
368 sys.exit(main()) # pragma: no cover