Commit afb6b894 authored by Daniel Scheffler's avatar Daniel Scheffler
Browse files

updated README and added documentation how to apply the same shifts to multiple images

parent 265b2fc2
......@@ -74,8 +74,8 @@ CR.calculate_spatial_shifts()
from py_tools_ds.ptds import GeoArray
from CoReg_Sat import COREG
im_reference = '/path/to/your/ref_image.bsq'
im_target = '/path/to/your/tgt_image.bsq'
#im_reference = '/path/to/your/ref_image.bsq'
#im_target = '/path/to/your/tgt_image.bsq'
# get a sample numpy array with corresponding geoinformation as reference image
geoArr = GeoArray(im_reference)
......@@ -106,16 +106,25 @@ CR.calculate_spatial_shifts()
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
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 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.
......@@ -125,6 +134,57 @@ CR.correct_shifts() returns an an OrderedDict containing the coregistered numpy
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],
...,
[ 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.ptds.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:
```python
from CoReg_Sat 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),
......@@ -139,6 +199,8 @@ CR.correct_shifts()
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],
......@@ -147,11 +209,13 @@ CR.correct_shifts()
...,
[ 0, 0, 0, ..., 755, 763, 773],
[ 0, 0, 0, ..., 760, 763, 749],
[9999, 9999, 9999, ..., 9999, 9999, 9999]], dtype=uint16))])
[9999, 9999, 9999, ..., 9999, 9999, 9999]], dtype=uint16)),
('GeoArray_shifted',
<py_tools_ds.ptds.io.raster.GeoArray.GeoArray at 0x7f6c5a1caa58>)])
To write the coregistered image to disk, the COREG class needs to be instanced with a filepath given to keyword 'path_out'.
Take a look at the keyword arguments of the DESHIFTER class when you need further adjustments (e.g. output paths for the corrected images; aligned output grid, ...).
### Shell console interface
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment