README.md 45.5 KB
Newer Older
Daniel Scheffler's avatar
Daniel Scheffler committed
1

Daniel Scheffler's avatar
Daniel Scheffler committed
2
[![logo](docs/images/arosics_logo.png)](https://gitext.gfz-potsdam.de/danschef/arosics/)
3

Daniel Scheffler's avatar
Daniel Scheffler committed
4
5
6
### An Automated and Robust Open-Source Image Co-Registration Software for Multi-Sensor Satellite Data


Daniel Scheffler's avatar
Daniel Scheffler committed
7
* Free software: GNU General Public License v3
Daniel Scheffler's avatar
Daniel Scheffler committed
8
* Documentation: http://danschef.gitext.gfz-potsdam.de/arosics/doc/
9
* The (open-access) paper corresponding to this software repository can be found here:
Daniel Scheffler's avatar
Daniel Scheffler committed
10
11
[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)

Daniel Scheffler's avatar
Daniel Scheffler committed
12
13
14
15
16


### Status

[![build status](https://gitext.gfz-potsdam.de/danschef/arosics/badges/master/build.svg)](https://gitext.gfz-potsdam.de/danschef/arosics/commits/master)
Daniel Scheffler's avatar
Daniel Scheffler committed
17
[![coverage report](https://gitext.gfz-potsdam.de/danschef/arosics/badges/master/coverage.svg)](http://danschef.gitext.gfz-potsdam.de/arosics/coverage/)
Daniel Scheffler's avatar
Daniel Scheffler committed
18
[![pypi_status](https://img.shields.io/pypi/v/arosics.svg)](https://pypi.python.org/pypi/arosics)
19
[![license](https://img.shields.io/pypi/l/arosics.svg)](https://gitext.gfz-potsdam.de/danschef/arosics/blob/master/LICENSE)
Daniel Scheffler's avatar
Daniel Scheffler committed
20
[![python versions](https://img.shields.io/pypi/pyversions/arosics.svg)](https://img.shields.io/pypi/pyversions/arosics.svg)
21

Daniel Scheffler's avatar
Daniel Scheffler committed
22

23
See also the latest [coverage](http://danschef.gitext.gfz-potsdam.de/arosics/coverage/) and the [nosetests](http://danschef.gitext.gfz-potsdam.de/arosics/nosetests_reports/nosetests.html) HTML report.
Daniel Scheffler's avatar
Daniel Scheffler committed
24
25
26
27
28



# Description

29
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).
Daniel Scheffler's avatar
Daniel Scheffler committed
30

Daniel Scheffler's avatar
Daniel Scheffler committed
31
AROSICS detects and corrects local as well as global misregistrations between two input images in the subpixel scale, that are often present in satellite imagery.
Daniel Scheffler's avatar
Daniel Scheffler committed
32
33

Prerequisites and hints:
34
35
36
37
38
39
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).
40

41
42
* 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.
43
44


45
46
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).

Daniel Scheffler's avatar
Daniel Scheffler committed
47
For further details regarding the implemented algorithm, example use cases, quality assessment and benchmarks refer to the above mentioned paper ([Scheffler et al. 2017](http://www.mdpi.com/2072-4292/9/7/676)).
Daniel Scheffler's avatar
Daniel Scheffler committed
48
49


50

Daniel Scheffler's avatar
Daniel Scheffler committed
51
52
53
54

Installation
------------

55
56
57
58
59
60
AROSICS depends on some open source packages which are usually installed without problems by the automatic install
routine. However, for some projects, we strongly recommend resolving the dependency before the automatic installer
is run. This approach avoids problems with conflicting versions of the same software.
Using [conda](https://conda.io/docs/), the recommended approach is:

```bash
61
# create virtual environment for arosics, this is optional
62
63
conda create -y -q --name arosics python=3
source activate arosics
64
65
66
67
68
69
conda install -y -q -c conda-forge numpy gdal scikit-image matplotlib pyproj rasterio shapely geopandas

# optional libraries:
conda install -y -q -c conda-forge basemap pykrige
conda install -y -q -c conda-forge pyfftw  # Linux and MacOS
conda install -y -q -c jesserobertson pyfftw  # Windows
70
71
72
```

To install AROSICS, use the pip installer:
Daniel Scheffler's avatar
Daniel Scheffler committed
73

74
```bash
75
pip install arosics
76
```
Daniel Scheffler's avatar
Daniel Scheffler committed
77
78


79
Or clone the repository via GIT and update the PATH environment variable:
Daniel Scheffler's avatar
Daniel Scheffler committed
80

81
82
83
```bash
cd /your/installation/folder
git clone https://gitext.gfz-potsdam.de/danschef/arosics.git
Daniel Scheffler's avatar
Daniel Scheffler committed
84
85
86
git clone https://gitext.gfz-potsdam.de/danschef/geoarray.git
git clone https://gitext.gfz-potsdam.de/danschef/py_tools_ds.git
PATH=$PATH:/path/to/your/installation/folder/arosics:/path/to/your/installation/folder/geoarray:/path/to/your/installation/folder/py_tools_ds
87
```
Daniel Scheffler's avatar
Daniel Scheffler committed
88
89


90
AROSICS has been tested with Python 3.4+ and Python 2.7. It should be fully compatible to all Python versions above 2.7.
Daniel Scheffler's avatar
Daniel Scheffler committed
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
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221


# Modules

## CoReg

This module calculates spatial shifts and performs a global correction (based on a single matching window).

### Python Interface

#### calculate spatial shifts - with input data on disk


```python
from arosics import COREG

#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()
```

    Calculating actual data corner coordinates for reference image...
    Corner coordinates of reference image:
    	[[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, 5790510.0], [352270.0, 5900040.0], [409790.0, 5900040.0], [409790.0, 5790250.0], [319460.0, 5790250.0]]
    Matching window position (X,Y): 354223/5805559
    Detected integer shifts (X/Y):       0/-2
    Detected subpixel shifts (X/Y):      0.357885632465/0.433837319984
    Calculated total shifts in fft pixel units (X/Y):         0.357885632465/-1.56616268002
    Calculated total shifts in reference pixel units (X/Y):   0.357885632465/-1.56616268002
    Calculated total shifts in target pixel units (X/Y):      0.357885632465/-1.56616268002
    Calculated map shifts (X,Y):				  3.578856324660592 15.661626799963415
    Original map info: ['UTM', 1, 1, 300000.0, 5900040.0, 10.0, 10.0, 33, 'North', 'WGS-84']
    Updated map info:  ['UTM', 1, 1, '300003.57885632466', '5900055.6616268', 10.0, 10.0, 33, 'North', 'WGS-84']


#### calculate spatial shifts - without any disk access


```python
from geoarray import GeoArray
from arosics import COREG

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)

ref_ndarray = geoArr[:]            # numpy.ndarray with shape (10980, 10980)
ref_gt      = geoArr.geotransform  # GDAL geotransform: (300000.0, 10.0, 0.0, 5900040.0, 0.0, -10.0)
ref_prj     = geoArr.projection    # projection as WKT string ('PROJCS["WGS 84 / UTM zone 33N....')

# get a sample numpy array with corresponding geoinformation as target image
geoArr  = GeoArray(im_target)

tgt_ndarray = geoArr[:]            # numpy.ndarray with shape (10980, 10980)
tgt_gt      = geoArr.geotransform  # GDAL geotransform: (300000.0, 10.0, 0.0, 5900040.0, 0.0, -10.0)
tgt_prj     = geoArr.projection    # projection as WKT string ('PROJCS["WGS 84 / UTM zone 33N....')

# pass an instance of GeoArray to COREG and calculate spatial shifts
geoArr_reference = GeoArray(ref_ndarray, ref_gt, ref_prj)
geoArr_target    = GeoArray(tgt_ndarray, tgt_gt, tgt_prj)

CR = COREG(geoArr_reference, geoArr_target, wp=(354223, 5805559), ws=(256,256))
CR.calculate_spatial_shifts()
```

    Calculating actual data corner coordinates for reference image...
    Corner coordinates of reference image:
    	[[300000.0, 5848140.0], [409790.0, 5848140.0], [409790.0, 5790250.0], [300000.0, 5790250.0]]
    Calculating actual data corner coordinates for image to be shifted...
    Corner coordinates of image to be shifted:
    	[[300000.0, 5847770.0], [409790.0, 5847770.0], [409790.0, 5790250.0], [300000.0, 5790250.0]]
    Matching window position (X,Y): 354223/5805559
    Detected integer shifts (X/Y):                            0/-2
    Detected subpixel shifts (X/Y):                           0.357885632465/0.433837319984
    Calculated total shifts in fft pixel units (X/Y):         0.357885632465/-1.56616268002
    Calculated total shifts in reference pixel units (X/Y):   0.357885632465/-1.56616268002
    Calculated total shifts in target pixel units (X/Y):      0.357885632465/-1.56616268002
    Calculated map shifts (X,Y):				  3.578856324660592/15.661626799963415
    Calculated absolute shift vector length in map units:     16.065328089207995
    Calculated angle of shift vector in degrees from North:   192.8717191970359
    Original map info: ['UTM', 1, 1, 300000.0, 5900040.0, 10.0, 10.0, 33, 'North', 'WGS-84']
    Updated map info:  ['UTM', 1, 1, '300003.57885632466', '5900055.6616268', 10.0, 10.0, 33, 'North', 'WGS-84']





    'success'



#### correct shifts

CR.correct_shifts() returns an an OrderedDict containing the coregistered numpy array and its corresponding geoinformation.


```python
CR.correct_shifts()
```




    OrderedDict([('band', None),
                 ('is shifted', True),
                 ('is resampled', False),
                 ('updated map info',
                  ['UTM',
                   1,
                   1,
                   300003.57885632466,
                   5900025.6616268,
                   10.0,
                   10.0,
                   33,
                   'North',
                   'WGS-84']),
                 ('updated geotransform',
                  [300000.0, 10.0, 0.0, 5900040.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, ...,  953,  972, 1044],
                         [   0,    0,    0, ..., 1001,  973, 1019],
                         [   0,    0,    0, ...,  953,  985, 1020],
222
                         ...,
Daniel Scheffler's avatar
Daniel Scheffler committed
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
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
267
268
269
270
271
                         [   0,    0,    0, ...,  755,  763,  773],
                         [   0,    0,    0, ...,  760,  763,  749],
                         [9999, 9999, 9999, ..., 9999, 9999, 9999]], dtype=uint16)),
                 ('GeoArray_shifted',
                  <py_tools_ds.io.raster.GeoArray.GeoArray at 0x7f6c5a1cabe0>)])



To write the coregistered image to disk, the COREG class needs to be instanced with a filepath given to keyword 'path_out'. The output raster format can be any format supported by GDAL. Find a list of supported formats here: http://www.gdal.org/formats_list.html

#### 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.

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
from arosics import DESHIFTER

DESHIFTER(im_target1, CR.coreg_info).correct_shifts()
DESHIFTER(im_target2, CR.coreg_info).correct_shifts()
```




    OrderedDict([('band', None),
                 ('is shifted', True),
                 ('is resampled', False),
                 ('updated map info',
                  ['UTM',
                   1,
                   1,
                   300003.57885632466,
                   5900025.6616268,
                   10.0,
                   10.0,
                   33,
                   'North',
                   'WGS-84']),
                 ('updated geotransform',
                  [300000.0, 10.0, 0.0, 5900040.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, ...,  953,  972, 1044],
                         [   0,    0,    0, ..., 1001,  973, 1019],
                         [   0,    0,    0, ...,  953,  985, 1020],
272
                         ...,
Daniel Scheffler's avatar
Daniel Scheffler committed
273
274
275
276
277
278
279
280
                         [   0,    0,    0, ...,  755,  763,  773],
                         [   0,    0,    0, ...,  760,  763,  749],
                         [9999, 9999, 9999, ..., 9999, 9999, 9999]], dtype=uint16)),
                 ('GeoArray_shifted',
                  <py_tools_ds.io.raster.GeoArray.GeoArray at 0x7f6c5a1caa58>)])



281

Daniel Scheffler's avatar
Daniel Scheffler committed
282
283
284
285
286
287
288

### Shell console interface

The help instructions of the console interface can be accessed like this:


```bash
Daniel Scheffler's avatar
Daniel Scheffler committed
289
python arosics_cli.py -h
Daniel Scheffler's avatar
Daniel Scheffler committed
290
291
```

292
Follow these instructions to run AROSICS from a shell console. For example, the most simple call for a global
Daniel Scheffler's avatar
Daniel Scheffler committed
293
294
295
296
co-registration would be like this:


```bash
Daniel Scheffler's avatar
Daniel Scheffler committed
297
python arosics_cli.py global /path/to/your/ref_image.bsq /path/to/your/tgt_image.bsq
Daniel Scheffler's avatar
Daniel Scheffler committed
298
299
300
```


301
302


Daniel Scheffler's avatar
Daniel Scheffler committed
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336

## 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.

### Python interface

#### detect and correct local shifts - with input data on disk


```python
from arosics import COREG_LOCAL

im_reference = '/path/to/your/ref_image.bsq'
im_target    = '/path/to/your/tgt_image.bsq'
kwargs = {
    'grid_res'     : 200,
    'window_size'  : (64,64),
    'path_out'     : 'auto',
    'projectDir'   : 'my_project',
    'q'            : False,
}

CRL = COREG_LOCAL(im_reference,im_target,**kwargs)
CRL.correct_shifts()
```

    Calculating actual data corner coordinates for reference image...
    Corner coordinates of reference image:
    	[[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, 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
337
    Calculating tie point grid (1977 points) in mode 'multiprocessing'...
Daniel Scheffler's avatar
Daniel Scheffler committed
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
    	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, 10979) to /home/gfz-fe/scheffler/jupyter/arosics_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.





    OrderedDict([('band', None),
                 ('is shifted', True),
                 ('is resampled', True),
                 ('updated map info',
                  ['UTM',
                   1,
                   1,
                   300000.0,
                   5900030.0,
                   10.0,
                   10.0,
                   33,
                   'North',
                   'WGS-84']),
                 ('updated geotransform',
                  [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, ..., 1034,  996, 1001],
                         [   0,    0,    0, ..., 1046, 1114, 1124],
                         [   0,    0,    0, ..., 1021, 1126, 1148],
370
                         ...,
Daniel Scheffler's avatar
Daniel Scheffler committed
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
                         [   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.io.raster.GeoArray.GeoArray at 0x7f451ac14a90>)])



#### detect and correct local shifts - without any disk access

All you have to do is to instanciate COREG_LOCAL with two instances of the GeoArray class as described above.


```python
CRL = COREG_LOCAL(GeoArray(ref_ndarray, ref_gt, ref_prj),GeoArray(tgt_ndarray, tgt_gt, tgt_prj),**kwargs)
CRL.correct_shifts()
```

389
#### visualize tie point grid with INITIAL shifts present in your input target image
Daniel Scheffler's avatar
Daniel Scheffler committed
390

391
Use the function COREG_LOCAL.view_CoRegPoints() to visualize the tie point grid with the calculated absolute lenghts of the shift vectors (the unit corresponds to the input projection - UTM in the shown example, thus the unit is 'meters'.).
Daniel Scheffler's avatar
Daniel Scheffler committed
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410

NOTE: a calculation of reliable shifts above cloud covered areas is not possible. In the current version of AROSICS these areas are not masked. A proper masking is planned.


```python
%matplotlib inline

CRL.view_CoRegPoints(figsize=(15,15),backgroundIm='ref')
```

    Note: array has been downsampled to 1000 x 1000 for faster visualization.



![png](docs/images/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.

411
#### visualize tie point grid with shifts present AFTER shift correction
Daniel Scheffler's avatar
Daniel Scheffler committed
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428

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
CRL_after_corr = COREG_LOCAL(im_reference, CRL.path_out, **kwargs)
CRL_after_corr.view_CoRegPoints(figsize=(15,15),backgroundIm='ref')
```

    Calculating actual data corner coordinates for reference image...
    Corner coordinates of reference image:
    	[[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], [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.
429
    Calculating tie point grid (1977 points) in mode 'multiprocessing'...
Daniel Scheffler's avatar
Daniel Scheffler committed
430
431
432
433
434
435
436
437
438
    	progress: |==================================================| 100.0% [1977/1977] Complete 10.78 sek



![png](docs/images/output_44_1.png)


The output figure shows a significant reduction of geometric shifts.

439
#### show the points table of the calculated tie point grid
Daniel Scheffler's avatar
Daniel Scheffler committed
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000

NOTE: Point records where no valid match has been found are filled with -9999.


```python
CRL.CoRegPoints_table
```




<div>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>geometry</th>
      <th>POINT_ID</th>
      <th>X_IM</th>
      <th>Y_IM</th>
      <th>X_UTM</th>
      <th>Y_UTM</th>
      <th>X_WIN_SIZE</th>
      <th>Y_WIN_SIZE</th>
      <th>X_SHIFT_PX</th>
      <th>Y_SHIFT_PX</th>
      <th>X_SHIFT_M</th>
      <th>Y_SHIFT_M</th>
      <th>ABS_SHIFT</th>
      <th>ANGLE</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>0</th>
      <td>POINT (352000 5898040)</td>
      <td>81</td>
      <td>5200</td>
      <td>200</td>
      <td>352000.0</td>
      <td>5898040.0</td>
      <td>64</td>
      <td>64</td>
      <td>-9999.000000</td>
      <td>-9999.000000</td>
      <td>-9999.000000</td>
      <td>-9999.000000</td>
      <td>-9999.000000</td>
      <td>-9999.000000</td>
    </tr>
    <tr>
      <th>1</th>
      <td>POINT (354000 5898040)</td>
      <td>82</td>
      <td>5400</td>
      <td>200</td>
      <td>354000.0</td>
      <td>5898040.0</td>
      <td>64</td>
      <td>64</td>
      <td>0.372470</td>
      <td>-0.285500</td>
      <td>3.724704</td>
      <td>2.855005</td>
      <td>4.693024</td>
      <td>232.529646</td>
    </tr>
    <tr>
      <th>2</th>
      <td>POINT (356000 5898040)</td>
      <td>83</td>
      <td>5600</td>
      <td>200</td>
      <td>356000.0</td>
      <td>5898040.0</td>
      <td>64</td>
      <td>64</td>
      <td>0.260948</td>
      <td>-0.293539</td>
      <td>2.609479</td>
      <td>2.935389</td>
      <td>3.927580</td>
      <td>221.636201</td>
    </tr>
    <tr>
      <th>3</th>
      <td>POINT (358000 5898040)</td>
      <td>84</td>
      <td>5800</td>
      <td>200</td>
      <td>358000.0</td>
      <td>5898040.0</td>
      <td>64</td>
      <td>64</td>
      <td>-9999.000000</td>
      <td>-9999.000000</td>
      <td>-9999.000000</td>
      <td>-9999.000000</td>
      <td>-9999.000000</td>
      <td>-9999.000000</td>
    </tr>
    <tr>
      <th>4</th>
      <td>POINT (360000 5898040)</td>
      <td>85</td>
      <td>6000</td>
      <td>200</td>
      <td>360000.0</td>
      <td>5898040.0</td>
      <td>64</td>
      <td>64</td>
      <td>-9999.000000</td>
      <td>-9999.000000</td>
      <td>-9999.000000</td>
      <td>-9999.000000</td>
      <td>-9999.000000</td>
      <td>-9999.000000</td>
    </tr>
    <tr>
      <th>5</th>
      <td>POINT (362000 5898040)</td>
      <td>86</td>
      <td>6200</td>
      <td>200</td>
      <td>362000.0</td>
      <td>5898040.0</td>
      <td>64</td>
      <td>64</td>
      <td>-9999.000000</td>
      <td>-9999.000000</td>
      <td>-9999.000000</td>
      <td>-9999.000000</td>
      <td>-9999.000000</td>
      <td>-9999.000000</td>
    </tr>
    <tr>
      <th>6</th>
      <td>POINT (364000 5898040)</td>
      <td>87</td>
      <td>6400</td>
      <td>200</td>
      <td>364000.0</td>
      <td>5898040.0</td>
      <td>64</td>
      <td>64</td>
      <td>0.141693</td>
      <td>0.187036</td>
      <td>1.416935</td>
      <td>-1.870360</td>
      <td>2.346476</td>
      <td>322.853405</td>
    </tr>
    <tr>
      <th>7</th>
      <td>POINT (366000 5898040)</td>
      <td>88</td>
      <td>6600</td>
      <td>200</td>
      <td>366000.0</td>
      <td>5898040.0</td>
      <td>64</td>
      <td>64</td>
      <td>-0.230941</td>
      <td>0.121139</td>
      <td>-2.309409</td>
      <td>-1.211389</td>
      <td>2.607841</td>
      <td>62.320969</td>
    </tr>
    <tr>
      <th>8</th>
      <td>POINT (368000 5898040)</td>
      <td>89</td>
      <td>6800</td>
      <td>200</td>
      <td>368000.0</td>
      <td>5898040.0</td>
      <td>64</td>
      <td>64</td>
      <td>-9999.000000</td>
      <td>-9999.000000</td>
      <td>-9999.000000</td>
      <td>-9999.000000</td>
      <td>-9999.000000</td>
      <td>-9999.000000</td>
    </tr>
    <tr>
      <th>9</th>
      <td>POINT (370000 5898040)</td>
      <td>90</td>
      <td>7000</td>
      <td>200</td>
      <td>370000.0</td>
      <td>5898040.0</td>
      <td>64</td>
      <td>64</td>
      <td>-0.035693</td>
      <td>0.084596</td>
      <td>-0.356928</td>
      <td>-0.845957</td>
      <td>0.918172</td>
      <td>22.875994</td>
    </tr>
    <tr>
      <th>10</th>
      <td>POINT (372000 5898040)</td>
      <td>91</td>
      <td>7200</td>
      <td>200</td>
      <td>372000.0</td>
      <td>5898040.0</td>
      <td>64</td>
      <td>64</td>
      <td>-9999.000000</td>
      <td>-9999.000000</td>
      <td>-9999.000000</td>
      <td>-9999.000000</td>
      <td>-9999.000000</td>
      <td>-9999.000000</td>
    </tr>
    <tr>
      <th>11</th>
      <td>POINT (374000 5898040)</td>
      <td>92</td>
      <td>7400</td>
      <td>200</td>
      <td>374000.0</td>
      <td>5898040.0</td>
      <td>64</td>
      <td>64</td>
      <td>-9999.000000</td>
      <td>-9999.000000</td>
      <td>-9999.000000</td>
      <td>-9999.000000</td>
      <td>-9999.000000</td>
      <td>-9999.000000</td>
    </tr>
    <tr>
      <th>12</th>
      <td>POINT (376000 5898040)</td>
      <td>93</td>
      <td>7600</td>
      <td>200</td>
      <td>376000.0</td>
      <td>5898040.0</td>
      <td>64</td>
      <td>64</td>
      <td>-9999.000000</td>
      <td>-9999.000000</td>
      <td>-9999.000000</td>
      <td>-9999.000000</td>
      <td>-9999.000000</td>
      <td>-9999.000000</td>
    </tr>
    <tr>
      <th>13</th>
      <td>POINT (378000 5898040)</td>
      <td>94</td>
      <td>7800</td>
      <td>200</td>
      <td>378000.0</td>
      <td>5898040.0</td>
      <td>64</td>
      <td>64</td>
      <td>-9999.000000</td>
      <td>-9999.000000</td>
      <td>-9999.000000</td>
      <td>-9999.000000</td>
      <td>-9999.000000</td>
      <td>-9999.000000</td>
    </tr>
    <tr>
      <th>14</th>
      <td>POINT (380000 5898040)</td>
      <td>95</td>
      <td>8000</td>
      <td>200</td>
      <td>380000.0</td>
      <td>5898040.0</td>
      <td>64</td>
      <td>64</td>
      <td>-9999.000000</td>
      <td>-9999.000000</td>
      <td>-9999.000000</td>
      <td>-9999.000000</td>
      <td>-9999.000000</td>
      <td>-9999.000000</td>
    </tr>
    <tr>
      <th>15</th>
      <td>POINT (382000 5898040)</td>
      <td>96</td>
      <td>8200</td>
      <td>200</td>
      <td>382000.0</td>
      <td>5898040.0</td>
      <td>64</td>
      <td>64</td>
      <td>-9999.000000</td>
      <td>-9999.000000</td>
      <td>-9999.000000</td>
      <td>-9999.000000</td>
      <td>-9999.000000</td>
      <td>-9999.000000</td>
    </tr>
    <tr>
      <th>16</th>
      <td>POINT (384000 5898040)</td>
      <td>97</td>
      <td>8400</td>
      <td>200</td>
      <td>384000.0</td>
      <td>5898040.0</td>
      <td>64</td>
      <td>64</td>
      <td>-9999.000000</td>
      <td>-9999.000000</td>
      <td>-9999.000000</td>
      <td>-9999.000000</td>
      <td>-9999.000000</td>
      <td>-9999.000000</td>
    </tr>
    <tr>
      <th>17</th>
      <td>POINT (386000 5898040)</td>
      <td>98</td>
      <td>8600</td>
      <td>200</td>
      <td>386000.0</td>
      <td>5898040.0</td>
      <td>64</td>
      <td>64</td>
      <td>-9999.000000</td>
      <td>-9999.000000</td>
      <td>-9999.000000</td>
      <td>-9999.000000</td>
      <td>-9999.000000</td>
      <td>-9999.000000</td>
    </tr>
    <tr>
      <th>18</th>
      <td>POINT (388000 5898040)</td>
      <td>99</td>
      <td>8800</td>
      <td>200</td>
      <td>388000.0</td>
      <td>5898040.0</td>
      <td>64</td>
      <td>64</td>
      <td>0.656098</td>
      <td>2.533985</td>
      <td>6.560977</td>
      <td>-25.339852</td>
      <td>26.175457</td>
      <td>345.483797</td>
    </tr>
    <tr>
      <th>19</th>
      <td>POINT (390000 5898040)</td>
      <td>100</td>
      <td>9000</td>
      <td>200</td>
      <td>390000.0</td>
      <td>5898040.0</td>
      <td>64</td>
      <td>64</td>
      <td>-9999.000000</td>
      <td>-9999.000000</td>
      <td>-9999.000000</td>
      <td>-9999.000000</td>
      <td>-9999.000000</td>
      <td>-9999.000000</td>
    </tr>
    <tr>
      <th>20</th>
      <td>POINT (392000 5898040)</td>
      <td>101</td>
      <td>9200</td>
      <td>200</td>
      <td>392000.0</td>
      <td>5898040.0</td>
      <td>64</td>
      <td>64</td>
      <td>-9999.000000</td>
      <td>-9999.000000</td>
      <td>-9999.000000</td>
      <td>-9999.000000</td>
      <td>-9999.000000</td>
      <td>-9999.000000</td>
    </tr>
    <tr>
      <th>21</th>
      <td>POINT (394000 5898040)</td>
      <td>102</td>
      <td>9400</td>
      <td>200</td>
      <td>394000.0</td>
      <td>5898040.0</td>
      <td>64</td>
      <td>64</td>
      <td>-9999.000000</td>
      <td>-9999.000000</td>
      <td>-9999.000000</td>
      <td>-9999.000000</td>
      <td>-9999.000000</td>
      <td>-9999.000000</td>
    </tr>
    <tr>
      <th>22</th>
      <td>POINT (396000 5898040)</td>
      <td>103</td>
      <td>9600</td>
      <td>200</td>
      <td>396000.0</td>
      <td>5898040.0</td>
      <td>64</td>
      <td>64</td>
      <td>-9999.000000</td>
      <td>-9999.000000</td>
      <td>-9999.000000</td>
      <td>-9999.000000</td>
      <td>-9999.000000</td>
      <td>-9999.000000</td>
    </tr>
    <tr>
      <th>23</th>
      <td>POINT (398000 5898040)</td>
      <td>104</td>
      <td>9800</td>
      <td>200</td>
      <td>398000.0</td>
      <td>5898040.0</td>
      <td>64</td>
      <td>64</td>
      <td>-9999.000000</td>
      <td>-9999.000000</td>
      <td>-9999.000000</td>
      <td>-9999.000000</td>
      <td>-9999.000000</td>
      <td>-9999.000000</td>
    </tr>
    <tr>
      <th>24</th>
      <td>POINT (400000 5898040)</td>
      <td>105</td>
      <td>10000</td>
      <td>200</td>
      <td>400000.0</td>
      <td>5898040.0</td>
      <td>64</td>
      <td>64</td>
      <td>-0.147210</td>
      <td>-0.223871</td>
      <td>-1.472098</td>
      <td>2.238708</td>
      <td>2.679344</td>
      <td>146.672433</td>
    </tr>
    <tr>
      <th>25</th>
      <td>POINT (402000 5898040)</td>
      <td>106</td>
      <td>10200</td>
      <td>200</td>
      <td>402000.0</td>
      <td>5898040.0</td>
      <td>64</td>
      <td>64</td>
      <td>-9999.000000</td>
      <td>-9999.000000</td>
      <td>-9999.000000</td>
      <td>-9999.000000</td>
      <td>-9999.000000</td>
      <td>-9999.000000</td>
    </tr>
    <tr>
      <th>26</th>
      <td>POINT (404000 5898040)</td>
      <td>107</td>
      <td>10400</td>
      <td>200</td>
      <td>404000.0</td>
      <td>5898040.0</td>
      <td>64</td>
      <td>64</td>
      <td>-9999.000000</td>
      <td>-9999.000000</td>
      <td>-9999.000000</td>
      <td>-9999.000000</td>
      <td>-9999.000000</td>
      <td>-9999.000000</td>
    </tr>
    <tr>
      <th>27</th>
      <td>POINT (406000 5898040)</td>
      <td>108</td>
      <td>10600</td>
      <td>200</td>
      <td>406000.0</td>
      <td>5898040.0</td>
      <td>64</td>
      <td>64</td>
      <td>0.249318</td>
      <td>0.214416</td>
      <td>2.493182</td>
      <td>-2.144158</td>
      <td>3.288369</td>
      <td>310.695805</td>
    </tr>
    <tr>
      <th>28</th>
      <td>POINT (408000 5898040)</td>
      <td>109</td>
      <td>10800</td>
      <td>200</td>
      <td>408000.0</td>
      <td>5898040.0</td>
      <td>64</td>
      <td>64</td>
      <td>0.372511</td>
      <td>-1.410450</td>
      <td>3.725107</td>
      <td>14.104504</td>
      <td>14.588127</td>
      <td>194.794441</td>
    </tr>
    <tr>
      <th>29</th>
      <td>POINT (352000 5896040)</td>
      <td>136</td>
      <td>5200</td>
      <td>400</td>
      <td>352000.0</td>
      <td>5896040.0</td>
      <td>64</td>
      <td>64</td>
      <td>-9999.000000</td>
      <td>-9999.000000</td>
      <td>-9999.000000</td>
      <td>-9999.000000</td>
      <td>-9999.000000</td>
      <td>-9999.000000</td>
    </tr>
    <tr>
      <th>...</th>
      <td>...</td>
      <td>...</td>
      <td>...</td>
      <td>...</td>
      <td>...</td>
      <td>...</td>
      <td>...</td>
      <td>...</td>
      <td>...</td>
      <td>...</td>
      <td>...</td>
      <td>...</td>
      <td>...</td>
      <td>...</td>
    </tr>
    <tr>
For faster browsing, not all history is shown. View entire blame