Daniel Scheffler
arosics

Repository

- gdal, osr, ogr
- geopandas
- pykrige
- argparse
- shapely
- pyfftw is optional but will speed up calculation
- folium and geojson for some visualization functions
https://gitext.gfz-potsdam.de/danschef/py_tools_ds
from CoReg_Sat 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']
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'

# 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'
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>)])
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),
             ('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 0x7f6c5a1caa58>)])
cd /path/to/CoReg_Sat/
python ./dsc__CoReg_Sat_FourierShiftTheorem.py -h
python ./dsc__CoReg_Sat_FourierShiftTheorem.py /path/to/your/ref_image.bsq /path/to/your/tgt_image.bsq
from CoReg_Sat 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
Calculating geometric quality grid (1977 points) in mode 'multiprocessing'...
	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/CoReg_Sat_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],
                     ...,
                     [   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.ptds.io.raster.GeoArray.GeoArray at 0x7f451ac14a90>)])
CRL = COREG_LOCAL(GeoArray(ref_ndarray, ref_gt, ref_prj),GeoArray(tgt_ndarray, tgt_gt, tgt_prj),**kwargs)
CRL.correct_shifts()
%matplotlib inline

CRL.view_CoRegPoints(figsize=(15,15),backgroundIm='ref')
Note: array has been downsampled to 1000 x 1000 for faster visualization.
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.
Calculating geometric quality grid (1977 points) in mode 'multiprocessing'...
	progress: |==================================================| 100.0% [1977/1977] Complete 10.78 sek
CRL.CoRegPoints_table
CRL.quality_grid.to_PointShapefile(path_out='/path/to/your/output_shapefile.shp')