Figure 6 script¶

This notebook requires additional steps to be reproduced:

Make sure you have installed all necessary packages, including vmap and others specified in the cell below.
Update the input file paths so they point to the actual files on your local machine. These files include the source Landsat 8 image (LC08_L1TP_061018_20180304_20180319_01_T1_B8_s.TIF, downloadable at https://doi.org/10.17605/OSF.IO/HE7YR) and the RGI glacier outline (01_rgi60_Alaska.shp, downloadable at https://doi.org/10.7265/4m1f-gd79).

import numpy as np
import matplotlib.pyplot as plt
import os,sys,glob
from pygeotools.lib import geolib,iolib,warplib,malib
from imview import pltlib
import rasterio
import geopandas as gpd

%matplotlib inline

%cd /nobackup/sbhusha1/feature_tracking_wg

img_fn = 'Cropped/LC08_L1TP_061018_20180304_20180319_01_T1_B8_s.TIF'

def create_synthetic_offset(imgfile, mode='subpixel', block_size=500):
    """
    imgfile: str, geotiff file path
    mode: 'subpixel' or 'multipixel'
    block_size: int, increment block size
    ----
    returns:
    shift_arx: np.ndarray, offset field (x), in pixels
    shift_ary: np.ndarray, offset field (y), in pixels
    ----
    """  
    with rasterio.open(imgfile) as src:
        data_shape = (src.height, src.width)
    idxy, idxx = np.indices(data_shape)   
    # for Numpy array, first is row element (-> geotiff's y direction, height) 
    # and second is column element (-> geotiff's x direction, width)
    
    if mode == 'subpixel':
        shift_arx = idxx // block_size
        shift_arx = 0.1 * shift_arx + 0.1
        shift_ary = idxy // block_size
        shift_ary = -0.1 * shift_ary - 0.1
    elif mode == 'multipixel':
        shift_arx = 1 + idxx // block_size
        shift_ary = -1 - idxy // block_size
    else:
        raise ValueError('Mode is not defined.')
        
    return shift_arx, shift_ary

def apply_synthetic_offset(imgfile, shift_arx, shift_ary, spline_order=1):
    """
    imgfile: str, geotiff file path
    shift_arx: np.ndarray, offset field (x) from gftt.create_synthetic_offset
    shift_ary: np.ndarray, offset field (y) from gftt.create_synthetic_offset
    ----
    returns:
    ----
    """ 
    import rasterio
    from scipy.ndimage import map_coordinates
    with rasterio.open(imgfile) as src:
        data_shape = (src.height, src.width)
        data = src.read(1)
    idxy, idxx = np.indices(data_shape)
    shifted_y = idxy + shift_ary
    shifted_x = idxx + shift_arx
    shifted_yx = np.vstack((shifted_y.flatten(), shifted_x.flatten()))
    
    shifted_val = map_coordinates(data, shifted_yx, order=spline_order, mode='nearest')
    shifted_val = np.reshape(shifted_val, data_shape)
    
    return shifted_val

data_shape = img.shape
shift_arx, shift_ary = create_synthetic_offset(img_fn)
shift_arx

array([[0.1, 0.1, 0.1, ..., 0.8, 0.8, 0.8],
       [0.1, 0.1, 0.1, ..., 0.8, 0.8, 0.8],
       [0.1, 0.1, 0.1, ..., 0.8, 0.8, 0.8],
       ...,
       [0.1, 0.1, 0.1, ..., 0.8, 0.8, 0.8],
       [0.1, 0.1, 0.1, ..., 0.8, 0.8, 0.8],
       [0.1, 0.1, 0.1, ..., 0.8, 0.8, 0.8]])

shifted_image = apply_synthetic_offset(img_fn, shift_arx, shift_ary)
outfn = os.path.splitext(img_fn)[0]+'_subpixel_shift.tif'
iolib.writeGTiff(shifted_image,outfn,src_ds=ds_list[2],ndv=0)

f,ax = plt.subplots()
pltlib.iv(iolib.fn_getma(outfn),ax=ax,cmap='gray')

<Axes: >

../../_images/feature_tracking_figure6_8_1.png

!vmap.py $img_fn $outfn -dt none

2023-07-08 02:11:05.369578
2023-07-08 09:11:05.369647 UTC


Warping all inputs to the following:
Resolution: 15.0
Extent: [584872.5, 6717982.5, 641572.5, 6755182.5]
Projection: '+proj=utm +zone=7 +datum=WGS84 +units=m +no_defs'
Resampling alg: average

1 of 2: Cropped/LC08_L1TP_061018_20180304_20180319_01_T1_B8_s.TIF
2 of 2: Cropped/LC08_L1TP_061018_20180304_20180319_01_T1_B8_s_subpixel_shift.tif
/nobackupp16/swbuild3/sbhusha1/StereoPipeline-3.2.1-alpha-2023-04-19-x86_64-Linux/bin/stereo_pprc -t pinhole --threads 56 --alignment-method None --corr-kernel 35 35 --corr-max-levels 5 --corr-timeout 1200 --subpixel-mode 1 --subpixel-kernel 35 35 --stereo-algorithm asp_bm --erode-max-size 1024 LC08_L1TP_061018_20180304_20180319_01_T1_B8_s__LC08_L1TP_061018_20180304_20180319_01_T1_B8_s_subpixel_shift_vmap_minm_35px_spm1/LC08_L1TP_061018_20180304_20180319_01_T1_B8_s_warp.tif LC08_L1TP_061018_20180304_20180319_01_T1_B8_s__LC08_L1TP_061018_20180304_20180319_01_T1_B8_s_subpixel_shift_vmap_minm_35px_spm1/LC08_L1TP_061018_20180304_20180319_01_T1_B8_s_subpixel_shift_warp.tif /nobackupp16/swbuild3/sbhusha1/pip_git_sw/vmap/vmap/dummy.tsai /nobackupp16/swbuild3/sbhusha1/pip_git_sw/vmap/vmap/dummy2.tsai LC08_L1TP_061018_20180304_20180319_01_T1_B8_s__LC08_L1TP_061018_20180304_20180319_01_T1_B8_s_subpixel_shift_vmap_minm_35px_spm1/vmap

[ 2023-Jul-08 02:11:24 ] : Stage 0 --> PREPROCESSING 
	--> Setting number of processing threads to: 56
Stereo file ./stereo.default could not be found. Will use default settings and command line options only.
Writing log info to: LC08_L1TP_061018_20180304_20180319_01_T1_B8_s__LC08_L1TP_061018_20180304_20180319_01_T1_B8_s_subpixel_shift_vmap_minm_35px_spm1/vmap-log-stereo_pprc-07-08-0211-78852.txt
Using session: pinhole
Loading camera model: LC08_L1TP_061018_20180304_20180319_01_T1_B8_s__LC08_L1TP_061018_20180304_20180319_01_T1_B8_s_subpixel_shift_vmap_minm_35px_spm1/LC08_L1TP_061018_20180304_20180319_01_T1_B8_s_warp.tif /nobackupp16/swbuild3/sbhusha1/pip_git_sw/vmap/vmap/dummy.tsai
Loading camera model: LC08_L1TP_061018_20180304_20180319_01_T1_B8_s__LC08_L1TP_061018_20180304_20180319_01_T1_B8_s_subpixel_shift_vmap_minm_35px_spm1/LC08_L1TP_061018_20180304_20180319_01_T1_B8_s_subpixel_shift_warp.tif /nobackupp16/swbuild3/sbhusha1/pip_git_sw/vmap/vmap/dummy2.tsai
Georef 1: -- Proj.4 Geospatial Reference Object --
	PROJCS name: WGS 84 / UTM zone 7N
	Transform  : Matrix3x3((15,0,584872)(0,-15,6.75518e+06)(0,0,1))
	Geodetic Datum --> Name: WGS_1984  Spheroid: WGS 84  Semi-major axis: 6378137  Semi-minor axis: 6356752.3142451793  Meridian: Greenwich at 0  Proj4 Str: +datum=WGS84
	Proj.4 String: +proj=utm +zone=7 +units=m
	Pixel Interpretation: pixel as point
longitude range: [-180, 180]

Georef 2: -- Proj.4 Geospatial Reference Object --
	PROJCS name: WGS 84 / UTM zone 7N
	Transform  : Matrix3x3((15,0,584872)(0,-15,6.75518e+06)(0,0,1))
	Geodetic Datum --> Name: WGS_1984  Spheroid: WGS 84  Semi-major axis: 6378137  Semi-minor axis: 6356752.3142451793  Meridian: Greenwich at 0  Proj4 Str: +datum=WGS84
	Proj.4 String: +proj=utm +zone=7 +units=m
	Pixel Interpretation: pixel as point
longitude range: [-180, 180]

Warning: It appears that the input images are map-projected. In that case a DEM needs to be provided for stereo to give correct results.
Warning: Your cameras appear to be in the same location!
	You should double check your given camera
	models as most likely stereo won't be able
	to triangulate or perform epipolar rectification.
Distance between camera centers in meters: 0.
Using image files:  LC08_L1TP_061018_20180304_20180319_01_T1_B8_s__LC08_L1TP_061018_20180304_20180319_01_T1_B8_s_subpixel_shift_vmap_minm_35px_spm1/LC08_L1TP_061018_20180304_20180319_01_T1_B8_s_warp.tif, LC08_L1TP_061018_20180304_20180319_01_T1_B8_s__LC08_L1TP_061018_20180304_20180319_01_T1_B8_s_subpixel_shift_vmap_minm_35px_spm1/LC08_L1TP_061018_20180304_20180319_01_T1_B8_s_subpixel_shift_warp.tif
Using camera files: /nobackupp16/swbuild3/sbhusha1/pip_git_sw/vmap/vmap/dummy.tsai, /nobackupp16/swbuild3/sbhusha1/pip_git_sw/vmap/vmap/dummy2.tsai
Computing statistics for left
Using downsample scale: 4
	    Writing stats file: LC08_L1TP_061018_20180304_20180319_01_T1_B8_s__LC08_L1TP_061018_20180304_20180319_01_T1_B8_s_subpixel_shift_vmap_minm_35px_spm1/vmap-LC08_L1TP_061018_20180304_20180319_01_T1_B8_s_warp-stats.tif
	    left: [ lo: 6430 hi: 46285 mean: 18455.6 std_dev: 8616.76 ]
Computing statistics for right
Using downsample scale: 4
	    Writing stats file: LC08_L1TP_061018_20180304_20180319_01_T1_B8_s__LC08_L1TP_061018_20180304_20180319_01_T1_B8_s_subpixel_shift_vmap_minm_35px_spm1/vmap-LC08_L1TP_061018_20180304_20180319_01_T1_B8_s_subpixel_shift_warp-stats.tif
	    right: [ lo: 6444 hi: 45544 mean: 18432.6 std_dev: 8477.22 ]
	--> Applying alignment method: none
	--> Normalizing globally to: [1222.04 35689.1]
	--> Writing pre-aligned images.
	--> Writing: LC08_L1TP_061018_20180304_20180319_01_T1_B8_s__LC08_L1TP_061018_20180304_20180319_01_T1_B8_s_subpixel_shift_vmap_minm_35px_spm1/vmap-L.tif.
          L:  [******************************************************] Complete!
	--> Writing: LC08_L1TP_061018_20180304_20180319_01_T1_B8_s__LC08_L1TP_061018_20180304_20180319_01_T1_B8_s_subpixel_shift_vmap_minm_35px_spm1/vmap-R.tif.
          R:  [******************************************************] Complete!
	--> Generating image masks... 
Writing masks: LC08_L1TP_061018_20180304_20180319_01_T1_B8_s__LC08_L1TP_061018_20180304_20180319_01_T1_B8_s_subpixel_shift_vmap_minm_35px_spm1/vmap-lMask.tif LC08_L1TP_061018_20180304_20180319_01_T1_B8_s__LC08_L1TP_061018_20180304_20180319_01_T1_B8_s_subpixel_shift_vmap_minm_35px_spm1/vmap-rMask.tif.
         Mask L: [***************************************************] Complete!
         Mask R: [***************************************************] Complete!
	--> Creating previews. Subsampling by 0.489914 by using a tile of size 256 and 56 threads.
            Sub L: [*************************************************] Complete!
            Sub R: [*************************************************] Complete!
            Sub L Mask: [********************************************] Complete!
            Sub R Mask: [********************************************] Complete!

[ 2023-Jul-08 02:11:32 ] : PREPROCESSING FINISHED 
/nobackupp16/swbuild3/sbhusha1/StereoPipeline-3.2.1-alpha-2023-04-19-x86_64-Linux/bin/stereo_corr --compute-low-res-disparity-only -t pinhole --threads 56 --alignment-method None --corr-kernel 35 35 --corr-max-levels 5 --corr-timeout 1200 --subpixel-mode 1 --subpixel-kernel 35 35 --stereo-algorithm asp_bm --erode-max-size 1024 LC08_L1TP_061018_20180304_20180319_01_T1_B8_s__LC08_L1TP_061018_20180304_20180319_01_T1_B8_s_subpixel_shift_vmap_minm_35px_spm1/LC08_L1TP_061018_20180304_20180319_01_T1_B8_s_warp.tif LC08_L1TP_061018_20180304_20180319_01_T1_B8_s__LC08_L1TP_061018_20180304_20180319_01_T1_B8_s_subpixel_shift_vmap_minm_35px_spm1/LC08_L1TP_061018_20180304_20180319_01_T1_B8_s_subpixel_shift_warp.tif /nobackupp16/swbuild3/sbhusha1/pip_git_sw/vmap/vmap/dummy.tsai /nobackupp16/swbuild3/sbhusha1/pip_git_sw/vmap/vmap/dummy2.tsai LC08_L1TP_061018_20180304_20180319_01_T1_B8_s__LC08_L1TP_061018_20180304_20180319_01_T1_B8_s_subpixel_shift_vmap_minm_35px_spm1/vmap
	--> Setting number of processing threads to: 56
Stereo file ./stereo.default could not be found. Will use default settings and command line options only.
Writing log info to: LC08_L1TP_061018_20180304_20180319_01_T1_B8_s__LC08_L1TP_061018_20180304_20180319_01_T1_B8_s_subpixel_shift_vmap_minm_35px_spm1/vmap-log-stereo_corr-07-08-0211-79772.txt
Using session: pinhole
Loading camera model: LC08_L1TP_061018_20180304_20180319_01_T1_B8_s__LC08_L1TP_061018_20180304_20180319_01_T1_B8_s_subpixel_shift_vmap_minm_35px_spm1/LC08_L1TP_061018_20180304_20180319_01_T1_B8_s_warp.tif /nobackupp16/swbuild3/sbhusha1/pip_git_sw/vmap/vmap/dummy.tsai
Loading camera model: LC08_L1TP_061018_20180304_20180319_01_T1_B8_s__LC08_L1TP_061018_20180304_20180319_01_T1_B8_s_subpixel_shift_vmap_minm_35px_spm1/LC08_L1TP_061018_20180304_20180319_01_T1_B8_s_subpixel_shift_warp.tif /nobackupp16/swbuild3/sbhusha1/pip_git_sw/vmap/vmap/dummy2.tsai
Georef 1: -- Proj.4 Geospatial Reference Object --
	PROJCS name: WGS 84 / UTM zone 7N
	Transform  : Matrix3x3((15,0,584872)(0,-15,6.75518e+06)(0,0,1))
	Geodetic Datum --> Name: WGS_1984  Spheroid: WGS 84  Semi-major axis: 6378137  Semi-minor axis: 6356752.3142451793  Meridian: Greenwich at 0  Proj4 Str: +datum=WGS84
	Proj.4 String: +proj=utm +zone=7 +units=m
	Pixel Interpretation: pixel as point
longitude range: [-180, 180]

Georef 2: -- Proj.4 Geospatial Reference Object --
	PROJCS name: WGS 84 / UTM zone 7N
	Transform  : Matrix3x3((15,0,584872)(0,-15,6.75518e+06)(0,0,1))
	Geodetic Datum --> Name: WGS_1984  Spheroid: WGS 84  Semi-major axis: 6378137  Semi-minor axis: 6356752.3142451793  Meridian: Greenwich at 0  Proj4 Str: +datum=WGS84
	Proj.4 String: +proj=utm +zone=7 +units=m
	Pixel Interpretation: pixel as point
longitude range: [-180, 180]

Warning: It appears that the input images are map-projected. In that case a DEM needs to be provided for stereo to give correct results.
Warning: Your cameras appear to be in the same location!
	You should double check your given camera
	models as most likely stereo won't be able
	to triangulate or perform epipolar rectification.
Distance between camera centers in meters: 0.
	--> Using LOG pre-processing filter with 1.5 sigma blur.

[ 2023-Jul-08 02:11:35 ] : Stage 1 --> CORRELATION

[ 2023-Jul-08 02:11:35 ] : Stage 1 --> LOW-RESOLUTION CORRELATION
No IP file found, computing IP now.
	    * Detecting interest points.
	--> Matching interest points using homography.
	    Looking for IP in left image...
	    Using 559 interest points per tile (1024^2 px).
	    Detecting IP
	    Removing IP near nodata with radius 4
	    Building descriptors
	    Found interest points: 5000
	    Recording interest points to file: LC08_L1TP_061018_20180304_20180319_01_T1_B8_s__LC08_L1TP_061018_20180304_20180319_01_T1_B8_s_subpixel_shift_vmap_minm_35px_spm1/vmap-L.vwip
	    Looking for IP in right image...
	    Using 559 interest points per tile (1024^2 px).
	    Detecting IP
	    Removing IP near nodata with radius 4
	    Building descriptors
	    Found interest points: 5000
	    Recording interest points to file: LC08_L1TP_061018_20180304_20180319_01_T1_B8_s__LC08_L1TP_061018_20180304_20180319_01_T1_B8_s_subpixel_shift_vmap_minm_35px_spm1/vmap-R.vwip
	--> Uniqueness threshold: 0.8
           Matching: [***************************************************.] 100%
	    Matched points: 4103
	    Inlier threshold:                     200
	    RANSAC iterations:                    100
	    * Writing match file: LC08_L1TP_061018_20180304_20180319_01_T1_B8_s__LC08_L1TP_061018_20180304_20180319_01_T1_B8_s_subpixel_shift_vmap_minm_35px_spm1/vmap-L__R.match
	--> Using interest points to determine search window.
	    * Loading match file: LC08_L1TP_061018_20180304_20180319_01_T1_B8_s__LC08_L1TP_061018_20180304_20180319_01_T1_B8_s_subpixel_shift_vmap_minm_35px_spm1/vmap-L__R.match
Removed 21 outliers based on percentiles of differences of interest points with --outlier-removal-params.
D_sub search range: (Origin: (-8.75, -2.625) width: 17.5 height: 6.25) px
Writing: LC08_L1TP_061018_20180304_20180319_01_T1_B8_s__LC08_L1TP_061018_20180304_20180319_01_T1_B8_s_subpixel_shift_vmap_minm_35px_spm1/vmap-D_sub.tif
        --> Low-resolution disparity:[*******************************] Complete!
	--> Full-res search range based on D_sub: (Origin: (-3, -3) width: 6 height: 6)

[ 2023-Jul-08 02:11:53 ] : LOW-RESOLUTION CORRELATION FINISHED

[ 2023-Jul-08 02:11:53 ] : CORRELATION FINISHED
/nobackupp16/swbuild3/sbhusha1/StereoPipeline-3.2.1-alpha-2023-04-19-x86_64-Linux/bin/stereo_corr -t pinhole --threads 56 --alignment-method None --corr-kernel 35 35 --corr-max-levels 5 --corr-timeout 1200 --subpixel-mode 1 --subpixel-kernel 35 35 --stereo-algorithm asp_bm --erode-max-size 1024 LC08_L1TP_061018_20180304_20180319_01_T1_B8_s__LC08_L1TP_061018_20180304_20180319_01_T1_B8_s_subpixel_shift_vmap_minm_35px_spm1/LC08_L1TP_061018_20180304_20180319_01_T1_B8_s_warp.tif LC08_L1TP_061018_20180304_20180319_01_T1_B8_s__LC08_L1TP_061018_20180304_20180319_01_T1_B8_s_subpixel_shift_vmap_minm_35px_spm1/LC08_L1TP_061018_20180304_20180319_01_T1_B8_s_subpixel_shift_warp.tif /nobackupp16/swbuild3/sbhusha1/pip_git_sw/vmap/vmap/dummy.tsai /nobackupp16/swbuild3/sbhusha1/pip_git_sw/vmap/vmap/dummy2.tsai LC08_L1TP_061018_20180304_20180319_01_T1_B8_s__LC08_L1TP_061018_20180304_20180319_01_T1_B8_s_subpixel_shift_vmap_minm_35px_spm1/vmap
	--> Setting number of processing threads to: 56
Stereo file ./stereo.default could not be found. Will use default settings and command line options only.
Writing log info to: LC08_L1TP_061018_20180304_20180319_01_T1_B8_s__LC08_L1TP_061018_20180304_20180319_01_T1_B8_s_subpixel_shift_vmap_minm_35px_spm1/vmap-log-stereo_corr-07-08-0211-79938.txt
Using session: pinhole
Loading camera model: LC08_L1TP_061018_20180304_20180319_01_T1_B8_s__LC08_L1TP_061018_20180304_20180319_01_T1_B8_s_subpixel_shift_vmap_minm_35px_spm1/LC08_L1TP_061018_20180304_20180319_01_T1_B8_s_warp.tif /nobackupp16/swbuild3/sbhusha1/pip_git_sw/vmap/vmap/dummy.tsai
Loading camera model: LC08_L1TP_061018_20180304_20180319_01_T1_B8_s__LC08_L1TP_061018_20180304_20180319_01_T1_B8_s_subpixel_shift_vmap_minm_35px_spm1/LC08_L1TP_061018_20180304_20180319_01_T1_B8_s_subpixel_shift_warp.tif /nobackupp16/swbuild3/sbhusha1/pip_git_sw/vmap/vmap/dummy2.tsai
Georef 1: -- Proj.4 Geospatial Reference Object --
	PROJCS name: WGS 84 / UTM zone 7N
	Transform  : Matrix3x3((15,0,584872)(0,-15,6.75518e+06)(0,0,1))
	Geodetic Datum --> Name: WGS_1984  Spheroid: WGS 84  Semi-major axis: 6378137  Semi-minor axis: 6356752.3142451793  Meridian: Greenwich at 0  Proj4 Str: +datum=WGS84
	Proj.4 String: +proj=utm +zone=7 +units=m
	Pixel Interpretation: pixel as point
longitude range: [-180, 180]

Georef 2: -- Proj.4 Geospatial Reference Object --
	PROJCS name: WGS 84 / UTM zone 7N
	Transform  : Matrix3x3((15,0,584872)(0,-15,6.75518e+06)(0,0,1))
	Geodetic Datum --> Name: WGS_1984  Spheroid: WGS 84  Semi-major axis: 6378137  Semi-minor axis: 6356752.3142451793  Meridian: Greenwich at 0  Proj4 Str: +datum=WGS84
	Proj.4 String: +proj=utm +zone=7 +units=m
	Pixel Interpretation: pixel as point
longitude range: [-180, 180]

Warning: It appears that the input images are map-projected. In that case a DEM needs to be provided for stereo to give correct results.
Warning: Your cameras appear to be in the same location!
	You should double check your given camera
	models as most likely stereo won't be able
	to triangulate or perform epipolar rectification.
Distance between camera centers in meters: 0.
	--> Using LOG pre-processing filter with 1.5 sigma blur.

[ 2023-Jul-08 02:11:55 ] : Stage 1 --> CORRELATION

[ 2023-Jul-08 02:11:55 ] : Stage 1 --> LOW-RESOLUTION CORRELATION
Cached IP match file found: LC08_L1TP_061018_20180304_20180319_01_T1_B8_s__LC08_L1TP_061018_20180304_20180319_01_T1_B8_s_subpixel_shift_vmap_minm_35px_spm1/vmap-L__R.match
	--> Using interest points to determine search window.
	    * Loading match file: LC08_L1TP_061018_20180304_20180319_01_T1_B8_s__LC08_L1TP_061018_20180304_20180319_01_T1_B8_s_subpixel_shift_vmap_minm_35px_spm1/vmap-L__R.match
Removed 21 outliers based on percentiles of differences of interest points with --outlier-removal-params.
	--> Using cached low-resolution disparity: LC08_L1TP_061018_20180304_20180319_01_T1_B8_s__LC08_L1TP_061018_20180304_20180319_01_T1_B8_s_subpixel_shift_vmap_minm_35px_spm1/vmap-D_sub.tif

[ 2023-Jul-08 02:11:55 ] : LOW-RESOLUTION CORRELATION FINISHED
	--> Full-res search range based on D_sub: (Origin: (-3, -3) width: 6 height: 6)
	--------------------------------------------------
	   Kernel size:    Vector2(35,35)
	   Search range:   (Origin: (-3, -3) width: 6 height: 6)
	   Cost mode:      2
	--------------------------------------------------
Writing: LC08_L1TP_061018_20180304_20180319_01_T1_B8_s__LC08_L1TP_061018_20180304_20180319_01_T1_B8_s_subpixel_shift_vmap_minm_35px_spm1/vmap-D.tif
        --> Correlation :[*******************************************] Complete!

[ 2023-Jul-08 02:12:22 ] : CORRELATION FINISHED
/nobackupp16/swbuild3/sbhusha1/StereoPipeline-3.2.1-alpha-2023-04-19-x86_64-Linux/bin/stereo_rfne -t pinhole --threads 56 --alignment-method None --corr-kernel 35 35 --corr-max-levels 5 --corr-timeout 1200 --subpixel-mode 1 --subpixel-kernel 35 35 --stereo-algorithm asp_bm --erode-max-size 1024 LC08_L1TP_061018_20180304_20180319_01_T1_B8_s__LC08_L1TP_061018_20180304_20180319_01_T1_B8_s_subpixel_shift_vmap_minm_35px_spm1/LC08_L1TP_061018_20180304_20180319_01_T1_B8_s_warp.tif LC08_L1TP_061018_20180304_20180319_01_T1_B8_s__LC08_L1TP_061018_20180304_20180319_01_T1_B8_s_subpixel_shift_vmap_minm_35px_spm1/LC08_L1TP_061018_20180304_20180319_01_T1_B8_s_subpixel_shift_warp.tif /nobackupp16/swbuild3/sbhusha1/pip_git_sw/vmap/vmap/dummy.tsai /nobackupp16/swbuild3/sbhusha1/pip_git_sw/vmap/vmap/dummy2.tsai LC08_L1TP_061018_20180304_20180319_01_T1_B8_s__LC08_L1TP_061018_20180304_20180319_01_T1_B8_s_subpixel_shift_vmap_minm_35px_spm1/vmap

[ 2023-Jul-08 02:12:23 ] : Stage 2 --> REFINEMENT 
	--> Setting number of processing threads to: 56
Stereo file ./stereo.default could not be found. Will use default settings and command line options only.
Writing log info to: LC08_L1TP_061018_20180304_20180319_01_T1_B8_s__LC08_L1TP_061018_20180304_20180319_01_T1_B8_s_subpixel_shift_vmap_minm_35px_spm1/vmap-log-stereo_rfne-07-08-0212-80059.txt
Using session: pinhole
Loading camera model: LC08_L1TP_061018_20180304_20180319_01_T1_B8_s__LC08_L1TP_061018_20180304_20180319_01_T1_B8_s_subpixel_shift_vmap_minm_35px_spm1/LC08_L1TP_061018_20180304_20180319_01_T1_B8_s_warp.tif /nobackupp16/swbuild3/sbhusha1/pip_git_sw/vmap/vmap/dummy.tsai
Loading camera model: LC08_L1TP_061018_20180304_20180319_01_T1_B8_s__LC08_L1TP_061018_20180304_20180319_01_T1_B8_s_subpixel_shift_vmap_minm_35px_spm1/LC08_L1TP_061018_20180304_20180319_01_T1_B8_s_subpixel_shift_warp.tif /nobackupp16/swbuild3/sbhusha1/pip_git_sw/vmap/vmap/dummy2.tsai
Georef 1: -- Proj.4 Geospatial Reference Object --
	PROJCS name: WGS 84 / UTM zone 7N
	Transform  : Matrix3x3((15,0,584872)(0,-15,6.75518e+06)(0,0,1))
	Geodetic Datum --> Name: WGS_1984  Spheroid: WGS 84  Semi-major axis: 6378137  Semi-minor axis: 6356752.3142451793  Meridian: Greenwich at 0  Proj4 Str: +datum=WGS84
	Proj.4 String: +proj=utm +zone=7 +units=m
	Pixel Interpretation: pixel as point
longitude range: [-180, 180]

Georef 2: -- Proj.4 Geospatial Reference Object --
	PROJCS name: WGS 84 / UTM zone 7N
	Transform  : Matrix3x3((15,0,584872)(0,-15,6.75518e+06)(0,0,1))
	Geodetic Datum --> Name: WGS_1984  Spheroid: WGS 84  Semi-major axis: 6378137  Semi-minor axis: 6356752.3142451793  Meridian: Greenwich at 0  Proj4 Str: +datum=WGS84
	Proj.4 String: +proj=utm +zone=7 +units=m
	Pixel Interpretation: pixel as point
longitude range: [-180, 180]

Warning: It appears that the input images are map-projected. In that case a DEM needs to be provided for stereo to give correct results.
Warning: Your cameras appear to be in the same location!
	You should double check your given camera
	models as most likely stereo won't be able
	to triangulate or perform epipolar rectification.
Distance between camera centers in meters: 0.
Left image nodata: -32768
Right image nodata: -32768
	--> Using LOG pre-processing filter with 1.5 sigma blur.
	--> Using parabola subpixel mode.
Writing: LC08_L1TP_061018_20180304_20180319_01_T1_B8_s__LC08_L1TP_061018_20180304_20180319_01_T1_B8_s_subpixel_shift_vmap_minm_35px_spm1/vmap-RD.tif
        --> Refinement :[********************************************] Complete!

[ 2023-Jul-08 02:12:29 ] : REFINEMENT FINISHED 
/nobackupp16/swbuild3/sbhusha1/StereoPipeline-3.2.1-alpha-2023-04-19-x86_64-Linux/bin/stereo_fltr -t pinhole --threads 56 --alignment-method None --corr-kernel 35 35 --corr-max-levels 5 --corr-timeout 1200 --subpixel-mode 1 --subpixel-kernel 35 35 --stereo-algorithm asp_bm --erode-max-size 1024 LC08_L1TP_061018_20180304_20180319_01_T1_B8_s__LC08_L1TP_061018_20180304_20180319_01_T1_B8_s_subpixel_shift_vmap_minm_35px_spm1/LC08_L1TP_061018_20180304_20180319_01_T1_B8_s_warp.tif LC08_L1TP_061018_20180304_20180319_01_T1_B8_s__LC08_L1TP_061018_20180304_20180319_01_T1_B8_s_subpixel_shift_vmap_minm_35px_spm1/LC08_L1TP_061018_20180304_20180319_01_T1_B8_s_subpixel_shift_warp.tif /nobackupp16/swbuild3/sbhusha1/pip_git_sw/vmap/vmap/dummy.tsai /nobackupp16/swbuild3/sbhusha1/pip_git_sw/vmap/vmap/dummy2.tsai LC08_L1TP_061018_20180304_20180319_01_T1_B8_s__LC08_L1TP_061018_20180304_20180319_01_T1_B8_s_subpixel_shift_vmap_minm_35px_spm1/vmap

[ 2023-Jul-08 02:12:30 ] : Stage 3 --> FILTERING 
Warning: Hole-filling is disabled by default in stereo_fltr. It is suggested to use instead point2dem's analogous functionality. It can be re-enabled using --enable-fill-holes.
	--> Setting number of processing threads to: 56
Stereo file ./stereo.default could not be found. Will use default settings and command line options only.
Writing log info to: LC08_L1TP_061018_20180304_20180319_01_T1_B8_s__LC08_L1TP_061018_20180304_20180319_01_T1_B8_s_subpixel_shift_vmap_minm_35px_spm1/vmap-log-stereo_fltr-07-08-0212-80235.txt
Using session: pinhole
Loading camera model: LC08_L1TP_061018_20180304_20180319_01_T1_B8_s__LC08_L1TP_061018_20180304_20180319_01_T1_B8_s_subpixel_shift_vmap_minm_35px_spm1/LC08_L1TP_061018_20180304_20180319_01_T1_B8_s_warp.tif /nobackupp16/swbuild3/sbhusha1/pip_git_sw/vmap/vmap/dummy.tsai
Loading camera model: LC08_L1TP_061018_20180304_20180319_01_T1_B8_s__LC08_L1TP_061018_20180304_20180319_01_T1_B8_s_subpixel_shift_vmap_minm_35px_spm1/LC08_L1TP_061018_20180304_20180319_01_T1_B8_s_subpixel_shift_warp.tif /nobackupp16/swbuild3/sbhusha1/pip_git_sw/vmap/vmap/dummy2.tsai
Georef 1: -- Proj.4 Geospatial Reference Object --
	PROJCS name: WGS 84 / UTM zone 7N
	Transform  : Matrix3x3((15,0,584872)(0,-15,6.75518e+06)(0,0,1))
	Geodetic Datum --> Name: WGS_1984  Spheroid: WGS 84  Semi-major axis: 6378137  Semi-minor axis: 6356752.3142451793  Meridian: Greenwich at 0  Proj4 Str: +datum=WGS84
	Proj.4 String: +proj=utm +zone=7 +units=m
	Pixel Interpretation: pixel as point
longitude range: [-180, 180]

Georef 2: -- Proj.4 Geospatial Reference Object --
	PROJCS name: WGS 84 / UTM zone 7N
	Transform  : Matrix3x3((15,0,584872)(0,-15,6.75518e+06)(0,0,1))
	Geodetic Datum --> Name: WGS_1984  Spheroid: WGS 84  Semi-major axis: 6378137  Semi-minor axis: 6356752.3142451793  Meridian: Greenwich at 0  Proj4 Str: +datum=WGS84
	Proj.4 String: +proj=utm +zone=7 +units=m
	Pixel Interpretation: pixel as point
longitude range: [-180, 180]

Warning: It appears that the input images are map-projected. In that case a DEM needs to be provided for stereo to give correct results.
Warning: Your cameras appear to be in the same location!
	You should double check your given camera
	models as most likely stereo won't be able
	to triangulate or perform epipolar rectification.
Distance between camera centers in meters: 0.
	--> Cleaning up disparity map prior to filtering processes (1 pass).
Writing: LC08_L1TP_061018_20180304_20180319_01_T1_B8_s__LC08_L1TP_061018_20180304_20180319_01_T1_B8_s_subpixel_shift_vmap_minm_35px_spm1/vmap-GoodPixelMap.tif
        --> Good pixel map: [****************************************] Complete!
	--> Removing small blobs.
Writing: LC08_L1TP_061018_20180304_20180319_01_T1_B8_s__LC08_L1TP_061018_20180304_20180319_01_T1_B8_s_subpixel_shift_vmap_minm_35px_spm1/vmap-F.tif
        --> Filtering: [*********************************************] Complete!

[ 2023-Jul-08 02:12:44 ] : FILTERING FINISHED 

2023-07-08 02:12:44.657561
2023-07-08 09:12:44.657648 UTC

swig/python detected a memory leak of type 'OSRSpatialReferenceShadow *', no destructor found.
swig/python detected a memory leak of type 'OSRSpatialReferenceShadow *', no destructor found.
swig/python detected a memory leak of type 'OSRSpatialReferenceShadow *', no destructor found.
swig/python detected a memory leak of type 'OSRSpatialReferenceShadow *', no destructor found.
swig/python detected a memory leak of type 'OSRSpatialReferenceShadow *', no destructor found.
swig/python detected a memory leak of type 'OSRSpatialReferenceShadow *', no destructor found.
swig/python detected a memory leak of type 'OSRSpatialReferenceShadow *', no destructor found.
swig/python detected a memory leak of type 'OSRSpatialReferenceShadow *', no destructor found.
swig/python detected a memory leak of type 'OSRSpatialReferenceShadow *', no destructor found.
swig/python detected a memory leak of type 'OSRSpatialReferenceShadow *', no destructor found.
swig/python detected a memory leak of type 'OSRSpatialReferenceShadow *', no destructor found.
swig/python detected a memory leak of type 'OSRSpatialReferenceShadow *', no destructor found.
swig/python detected a memory leak of type 'OSRSpatialReferenceShadow *', no destructor found.
swig/python detected a memory leak of type 'OSRSpatialReferenceShadow *', no destructor found.
swig/python detected a memory leak of type 'OSRSpatialReferenceShadow *', no destructor found.
swig/python detected a memory leak of type 'OSRSpatialReferenceShadow *', no destructor found.

out_disp_fn = 'LC08_L1TP_061018_20180304_20180319_01_T1_B8_s__LC08_L1TP_061018_20180304_20180319_01_T1_B8_s_subpixel_shift_vmap_minm_35px_spm1/vmap-F.tif'
out_dx,out_dy = [iolib.fn_getma(out_disp_fn,b) for b in [1,2]]

f,ax = plt.subplots(1,2)
clim_x = (0.1,0.8)
clim_y = (-0.6,-0.1)
cb1 = ax[0].imshow(shift_arx,cmap='inferno',clim=clim_x)
cb2 = ax[1].imshow(shift_ary,cmap='inferno',clim=clim_y)

../../_images/feature_tracking_figure6_11_0.png

ds = iolib.fn_getds(out_disp_fn)
extent = geolib.ds_extent(ds)
fig_extent = [extent[0],extent[2],extent[1],extent[3]]

glac_shp = gpd.read_file('/nobackup/sbhusha1/reference_data/rgi60/regions/01_rgi60_Alaska.shp')

kaskwulsh_mask = glac_shp['RGIId'] == 'RGI60-01.16201'
kaskwulsh_shp = (glac_shp[kaskwulsh_mask]).to_crs("EPSG:32607")

kaskwulsh_shp.plot()

<Axes: >

../../_images/feature_tracking_figure6_15_1.png

f,axa = plt.subplots(2,2,figsize=(10,8),sharex=True,sharey=True)
ax = axa.ravel()
clim_x = (0.1,0.8)
clim_y = (-0.6,-0.1)
pltlib.iv(shift_arx,ax=ax[0],cmap='inferno',clim=clim_x,title='input shift in x (px)',extent=fig_extent)
pltlib.iv(shift_ary,ax=ax[1],cmap='inferno',clim=clim_y,title='input shift in y (px)',extent=fig_extent)
kaskwulsh_shp.plot(ax=ax[2],facecolor="None",edgecolor='green',linewidth=0.95)
kaskwulsh_shp.plot(ax=ax[3],facecolor="None",edgecolor='green',linewidth=0.95)
pltlib.iv(-1*out_dx,ax=ax[2],cmap='inferno',clim=clim_x,cbar=False,title='measured shift in x (px)',extent=fig_extent)
pltlib.add_cbar(ax[2],mappable=cb1)
pltlib.iv(-1*out_dy,ax=ax[3],cmap='inferno',clim=clim_y,cbar=False,title='measured shift in y (px)',extent=fig_extent)

pltlib.add_cbar(ax[3],cb2)
plt.tight_layout()
#f.savefig('figure6_manuscript.png',dpi=300,bbox_inches='tight', pad_inches=0.1)
f.savefig('/nobackup/sbhusha1/notebooks/velocity/figure6_manuscript_revised_kaskawulsh_only.png',dpi=300,bbox_inches='tight', pad_inches=0.1)

../../_images/feature_tracking_figure6_16_0.png

GLAFT manual and supporting material

Figure 6 script

Figure 6 script¶