Bump autoRIFT 1.0.7 -> 1.1.0

contrib/geo_autoRIFT/README.md

@@ -3,9 +3,9 @@
 
 
 [![Language](https://img.shields.io/badge/python-3.6%2B-blue.svg)](https://www.python.org/)
-[![Latest version](https://img.shields.io/badge/latest%20version-v1.0.7-yellowgreen.svg)](https://github.com/leiyangleon/autoRIFT/releases)
+[![Latest version](https://img.shields.io/badge/latest%20version-v1.0.8-yellowgreen.svg)](https://github.com/leiyangleon/autoRIFT/releases)
 [![License](https://img.shields.io/badge/License-Apache%202.0-blue.svg)](https://github.com/leiyangleon/autoRIFT/blob/master/LICENSE)
-[![Citation](https://img.shields.io/badge/DOI-10.5281/zenodo.4025445-blue)](https://doi.org/10.5281/zenodo.4025445)
+[![Citation](https://img.shields.io/badge/DOI-10.5281/zenodo.4211013-blue)](https://doi.org/10.5281/zenodo.4211013)
 
 
 
@@ -22,7 +22,7 @@ Copyright (C) 2019 California Institute of Technology.  Government Sponsorship A
 
 Link: https://github.com/leiyangleon/autoRIFT
 
-Citation: https://doi.org/10.5281/zenodo.4025445
+Citation: https://doi.org/10.5281/zenodo.4211013
 
 
 ## 1. Authors
@@ -58,6 +58,7 @@ This effort was funded by the NASA MEaSUREs program in contribution to the Inter
 * the current version can be installed with the ISCE software (that supports both Cartesian and radar coordinates) or as a standalone Python module (Cartesian coordinates only)
 * when used in combination with the Geogrid Python module (https://github.com/leiyangleon/Geogrid), autoRIFT can be used for feature tracking between image pair over a grid defined in an arbitrary geographic Cartesian (northing/easting) coordinate projection
 * when the grid is provided in geographic Cartesian (northing/easting) coordinates, outputs are returned in geocoded GeoTIFF image file format with the same EPSG projection code as input search grid
+* **[NEW]** the program now supports fetching optical images (Landsat-8 GeoTIFF and Sentinel-2 COG formats are included) as well as other inputs (e.g. DEM, slope, etc; all in GeoTIFF format) from either local machine or URL links. See the changes on the Geogrid [commands](https://github.com/leiyangleon/Geogrid). When using the autoRIFT commands below, users need to append a url flag: "-urlflag 1" for using URL links and performing coregistration, and "-urlflag 0" (default; can be omitted) for using coregistered files on local machine. 
 
 ## 4. Possible Future Development
 

contrib/geo_autoRIFT/autoRIFT/autoRIFT.py

@@ -56,8 +56,7 @@ class autoRIFT:
 #        import scipy.io as sio
         
         
-        if self.zeroMask is not None:
+        self.zeroMask = (self.I1 == 0)
-            self.zeroMask = (self.I1 == 0)
         
         kernel = np.ones((self.WallisFilterWidth,self.WallisFilterWidth), dtype=np.float32)
         
@@ -115,8 +114,6 @@ class autoRIFT:
         import numpy as np
         
 
-        if self.zeroMask is not None:
-            self.zeroMask = (self.I1 == 0)
 
         kernel = -np.ones((self.WallisFilterWidth,self.WallisFilterWidth), dtype=np.float32)
 
@@ -140,8 +137,8 @@ class autoRIFT:
         import numpy as np
         
         
-        if self.zeroMask is not None:
+
-            self.zeroMask = (self.I1 == 0)
+        self.zeroMask = (self.I1 == 0)
         
 #        pdb.set_trace()
 
@@ -161,8 +158,6 @@ class autoRIFT:
         
         
         
-        if self.zeroMask is not None:
-            self.zeroMask = (self.I1 == 0)
 
         sobelx = cv2.getDerivKernels(1,0,self.WallisFilterWidth)
         
@@ -194,8 +189,8 @@ class autoRIFT:
         import numpy as np
         
         
-        if self.zeroMask is not None:
+
-            self.zeroMask = (self.I1 == 0)
+        self.zeroMask = (self.I1 == 0)
 
         self.I1 = 20.0 * np.log10(self.I1)
         self.I1 = cv2.Laplacian(self.I1,-1,ksize=self.WallisFilterWidth,borderType=cv2.BORDER_CONSTANT)
@@ -216,19 +211,27 @@ class autoRIFT:
         import numpy as np
         
         if self.DataType == 0:
-#            validData = np.logical_not(np.isnan(self.I1))
+            if self.zeroMask is not None:
-            validData = np.isfinite(self.I1)
+    #            validData = np.logical_not(np.isnan(self.I1))
-            S1 = np.std(self.I1[validData])*np.sqrt(self.I1[validData].size/(self.I1[validData].size-1.0))
+                validData = np.isfinite(self.I1)
-            M1 = np.mean(self.I1[validData])
+                temp = self.I1[validData]
+            else:
+                temp = self.I1
+            S1 = np.std(temp)*np.sqrt(temp.size/(temp.size-1.0))
+            M1 = np.mean(temp)
             self.I1 = (self.I1 - (M1 - 3*S1)) / (6*S1) * (2**8 - 0)
             
 #            self.I1[np.logical_not(np.isfinite(self.I1))] = 0
             self.I1 = np.round(np.clip(self.I1, 0, 255)).astype(np.uint8)
 
-#            validData = np.logical_not(np.isnan(self.I2))
+            if self.zeroMask is not None:
-            validData = np.isfinite(self.I2)
+    #            validData = np.logical_not(np.isnan(self.I2))
-            S2 = np.std(self.I2[validData])*np.sqrt(self.I2[validData].size/(self.I2[validData].size-1.0))
+                validData = np.isfinite(self.I2)
-            M2 = np.mean(self.I2[validData])
+                temp = self.I2[validData]
+            else:
+                temp = self.I2
+            S2 = np.std(temp)*np.sqrt(temp.size/(temp.size-1.0))
+            M2 = np.mean(temp)
             self.I2 = (self.I2 - (M2 - 3*S2)) / (6*S2) * (2**8 - 0)
             
 #            self.I2[np.logical_not(np.isfinite(self.I2))] = 0
@@ -240,8 +243,11 @@ class autoRIFT:
                 self.zeroMask = None
         
         elif self.DataType == 1:
-            self.I1[np.logical_not(np.isfinite(self.I1))] = 0
+            
-            self.I2[np.logical_not(np.isfinite(self.I2))] = 0
+            if self.zeroMask is not None:
+                self.I1[np.logical_not(np.isfinite(self.I1))] = 0
+                self.I2[np.logical_not(np.isfinite(self.I2))] = 0
+            
             self.I1 = self.I1.astype(np.float32)
             self.I2 = self.I2.astype(np.float32)
             
@@ -409,9 +415,9 @@ class autoRIFT:
 #            pdb.set_trace()
 
             if self.I1.dtype == np.uint8:
-                DxC, DyC = arImgDisp_u(self.I2.copy(), self.I1.copy(), xGrid0C.copy(), yGrid0C.copy(), ChipSizeXC, ChipSizeYC, SearchLimitX0C.copy(), SearchLimitY0C.copy(), Dx0C.copy(), Dy0C.copy(), SubPixFlag, overSampleRatio)
+                DxC, DyC = arImgDisp_u(self.I2.copy(), self.I1.copy(), xGrid0C.copy(), yGrid0C.copy(), ChipSizeXC, ChipSizeYC, SearchLimitX0C.copy(), SearchLimitY0C.copy(), Dx0C.copy(), Dy0C.copy(), SubPixFlag, overSampleRatio, self.MultiThread)
             elif self.I1.dtype == np.float32:
-                DxC, DyC = arImgDisp_s(self.I2.copy(), self.I1.copy(), xGrid0C.copy(), yGrid0C.copy(), ChipSizeXC, ChipSizeYC, SearchLimitX0C.copy(), SearchLimitY0C.copy(), Dx0C.copy(), Dy0C.copy(), SubPixFlag, overSampleRatio)
+                DxC, DyC = arImgDisp_s(self.I2.copy(), self.I1.copy(), xGrid0C.copy(), yGrid0C.copy(), ChipSizeXC, ChipSizeYC, SearchLimitX0C.copy(), SearchLimitY0C.copy(), Dx0C.copy(), Dy0C.copy(), SubPixFlag, overSampleRatio, self.MultiThread)
             else:
                 sys.exit('invalid data type for the image pair which must be unsigned integer 8 or 32-bit float')
             
@@ -459,9 +465,9 @@ class autoRIFT:
             ChipSizeYF = np.float32(np.round(ChipSizeXF*self.ScaleChipSizeY/2)*2)
 #            pdb.set_trace()
             if self.I1.dtype == np.uint8:
-                DxF, DyF = arImgDisp_u(self.I2.copy(), self.I1.copy(), xGrid0.copy(), yGrid0.copy(), ChipSizeXF, ChipSizeYF, SearchLimitX0.copy(), SearchLimitY0.copy(), Dx00.copy(), Dy00.copy(), SubPixFlag, overSampleRatio)
+                DxF, DyF = arImgDisp_u(self.I2.copy(), self.I1.copy(), xGrid0.copy(), yGrid0.copy(), ChipSizeXF, ChipSizeYF, SearchLimitX0.copy(), SearchLimitY0.copy(), Dx00.copy(), Dy00.copy(), SubPixFlag, overSampleRatio, self.MultiThread)
             elif self.I1.dtype == np.float32:
-                DxF, DyF = arImgDisp_s(self.I2.copy(), self.I1.copy(), xGrid0.copy(), yGrid0.copy(), ChipSizeXF, ChipSizeYF, SearchLimitX0.copy(), SearchLimitY0.copy(), Dx00.copy(), Dy00.copy(), SubPixFlag, overSampleRatio)
+                DxF, DyF = arImgDisp_s(self.I2.copy(), self.I1.copy(), xGrid0.copy(), yGrid0.copy(), ChipSizeXF, ChipSizeYF, SearchLimitX0.copy(), SearchLimitY0.copy(), Dx00.copy(), Dy00.copy(), SubPixFlag, overSampleRatio, self.MultiThread)
             else:
                 sys.exit('invalid data type for the image pair which must be unsigned integer 8 or 32-bit float')
 
@@ -647,6 +653,7 @@ class autoRIFT:
         self.CoarseCorCutoff = 0.01
         self.OverSampleRatio = 16
         self.DataType = 0
+        self.MultiThread = 0
 
 
 
@@ -659,8 +666,23 @@ class AUTO_RIFT_CORE:
         self._autoriftcore = None
 
 
+var_dict = {}
 
-def arImgDisp_u(I1, I2, xGrid, yGrid, ChipSizeX, ChipSizeY, SearchLimitX, SearchLimitY, Dx0, Dy0, SubPixFlag, oversample):
+def initializer(I1, I2, xGrid, yGrid, SearchLimitX, SearchLimitY, ChipSizeX, ChipSizeY, Dx0, Dy0):
+    var_dict['I1'] = I1
+    var_dict['I2'] = I2
+    var_dict['xGrid'] = xGrid
+    var_dict['yGrid'] = yGrid
+    var_dict['SearchLimitX'] = SearchLimitX
+    var_dict['SearchLimitY'] = SearchLimitY
+    var_dict['ChipSizeX'] = ChipSizeX
+    var_dict['ChipSizeY'] = ChipSizeY
+    var_dict['Dx0'] = Dx0
+    var_dict['Dy0'] = Dy0
+
+
+
+def unpacking_loop_u(tup):
     
     import numpy as np
     from . import autoriftcore
@@ -671,6 +693,168 @@ def arImgDisp_u(I1, I2, xGrid, yGrid, ChipSizeX, ChipSizeY, SearchLimitX, Search
     
     core._autoriftcore = autoriftcore.createAutoRiftCore_Py()
     
+    k, chunkInds, SubPixFlag, oversample, in_shape, I_shape = tup
+
+    I1 = np.frombuffer(var_dict['I1'],dtype=np.uint8).reshape(I_shape)
+    I2 = np.frombuffer(var_dict['I2'],dtype=np.uint8).reshape(I_shape)
+    xGrid = np.frombuffer(var_dict['xGrid'],dtype=np.float32).reshape(in_shape)
+    yGrid = np.frombuffer(var_dict['yGrid'],dtype=np.float32).reshape(in_shape)
+    SearchLimitX = np.frombuffer(var_dict['SearchLimitX'],dtype=np.float32).reshape(in_shape)
+    SearchLimitY = np.frombuffer(var_dict['SearchLimitY'],dtype=np.float32).reshape(in_shape)
+    ChipSizeX = np.frombuffer(var_dict['ChipSizeX'],dtype=np.float32).reshape(in_shape)
+    ChipSizeY = np.frombuffer(var_dict['ChipSizeY'],dtype=np.float32).reshape(in_shape)
+    Dx0 = np.frombuffer(var_dict['Dx0'],dtype=np.float32).reshape(in_shape)
+    Dy0 = np.frombuffer(var_dict['Dy0'],dtype=np.float32).reshape(in_shape)
+    
+    Dx = np.empty(chunkInds.shape,dtype=np.float32)
+    Dx.fill(np.nan)
+    Dy = Dx.copy()
+
+
+#    print(k)
+#    print(np.min(chunkInds),np.max(chunkInds))
+#    print(chunkInds.shape)
+
+    for ind in chunkInds:
+        
+        ind1 = np.where(chunkInds == ind)[0][0]
+        
+        ii, jj = [v[0] for v in np.unravel_index([ind], in_shape)]
+
+        if (SearchLimitX[ii,jj] == 0) & (SearchLimitY[ii,jj] == 0):
+            continue
+    
+        # remember motion terms Dx and Dy correspond to I1 relative to I2 (reference)
+        clx = np.floor(ChipSizeX[ii,jj]/2)
+        ChipRangeX = slice(int(-clx - Dx0[ii,jj] + xGrid[ii,jj]) , int(clx - Dx0[ii,jj] + xGrid[ii,jj]))
+        cly = np.floor(ChipSizeY[ii,jj]/2)
+        ChipRangeY = slice(int(-cly - Dy0[ii,jj] + yGrid[ii,jj]) , int(cly - Dy0[ii,jj] + yGrid[ii,jj]))
+        ChipI = I2[ChipRangeY,ChipRangeX]
+        
+        SearchRangeX = slice(int(-clx - SearchLimitX[ii,jj] + xGrid[ii,jj]) , int(clx + SearchLimitX[ii,jj] - 1 + xGrid[ii,jj]))
+        SearchRangeY = slice(int(-cly - SearchLimitY[ii,jj] + yGrid[ii,jj]) , int(cly + SearchLimitY[ii,jj] - 1 + yGrid[ii,jj]))
+        RefI = I1[SearchRangeY,SearchRangeX]
+        
+        minChipI = np.min(ChipI)
+        if minChipI < 0:
+            ChipI = ChipI - minChipI
+        if np.all(ChipI == ChipI[0,0]):
+            continue
+        
+        minRefI = np.min(RefI)
+        if minRefI < 0:
+            RefI = RefI - minRefI
+        if np.all(RefI == RefI[0,0]):
+            continue
+
+
+        if SubPixFlag:
+            # call C++
+            Dx[ind1], Dy[ind1] = np.float32(autoriftcore.arSubPixDisp_u_Py(core._autoriftcore,ChipI.shape[1],ChipI.shape[0],ChipI.ravel(),RefI.shape[1],RefI.shape[0],RefI.ravel(),oversample))
+#                   # call Python
+#                   Dx1[ii], Dy1[ii] = arSubPixDisp(ChipI,RefI)
+        else:
+            # call C++
+            Dx[ind1], Dy[ind1] = np.float32(autoriftcore.arPixDisp_u_Py(core._autoriftcore,ChipI.shape[1],ChipI.shape[0],ChipI.ravel(),RefI.shape[1],RefI.shape[0],RefI.ravel()))
+#                   # call Python
+#                   Dx1[ii], Dy1[ii] = arPixDisp(ChipI,RefI)
+    return Dx, Dy
+
+
+def unpacking_loop_s(tup):
+    
+    import numpy as np
+    from . import autoriftcore
+    
+    core = AUTO_RIFT_CORE()
+    if core._autoriftcore is not None:
+        autoriftcore.destroyAutoRiftCore_Py(core._autoriftcore)
+    
+    core._autoriftcore = autoriftcore.createAutoRiftCore_Py()
+
+
+    k, chunkInds, SubPixFlag, oversample, in_shape, I_shape = tup
+
+    I1 = np.frombuffer(var_dict['I1'],dtype=np.float32).reshape(I_shape)
+    I2 = np.frombuffer(var_dict['I2'],dtype=np.float32).reshape(I_shape)
+    xGrid = np.frombuffer(var_dict['xGrid'],dtype=np.float32).reshape(in_shape)
+    yGrid = np.frombuffer(var_dict['yGrid'],dtype=np.float32).reshape(in_shape)
+    SearchLimitX = np.frombuffer(var_dict['SearchLimitX'],dtype=np.float32).reshape(in_shape)
+    SearchLimitY = np.frombuffer(var_dict['SearchLimitY'],dtype=np.float32).reshape(in_shape)
+    ChipSizeX = np.frombuffer(var_dict['ChipSizeX'],dtype=np.float32).reshape(in_shape)
+    ChipSizeY = np.frombuffer(var_dict['ChipSizeY'],dtype=np.float32).reshape(in_shape)
+    Dx0 = np.frombuffer(var_dict['Dx0'],dtype=np.float32).reshape(in_shape)
+    Dy0 = np.frombuffer(var_dict['Dy0'],dtype=np.float32).reshape(in_shape)
+
+    
+    Dx = np.empty(chunkInds.shape,dtype=np.float32)
+    Dx.fill(np.nan)
+    Dy = Dx.copy()
+
+#    print(k)
+#    print(np.min(chunkInds),np.max(chunkInds))
+#    print(chunkInds.shape)
+    
+    for ind in chunkInds:
+        
+        ind1 = np.where(chunkInds == ind)[0][0]
+        
+        ii, jj = [v[0] for v in np.unravel_index([ind], in_shape)]
+        
+        if (SearchLimitX[ii,jj] == 0) & (SearchLimitY[ii,jj] == 0):
+            continue
+        
+        # remember motion terms Dx and Dy correspond to I1 relative to I2 (reference)
+        clx = np.floor(ChipSizeX[ii,jj]/2)
+        ChipRangeX = slice(int(-clx - Dx0[ii,jj] + xGrid[ii,jj]) , int(clx - Dx0[ii,jj] + xGrid[ii,jj]))
+        cly = np.floor(ChipSizeY[ii,jj]/2)
+        ChipRangeY = slice(int(-cly - Dy0[ii,jj] + yGrid[ii,jj]) , int(cly - Dy0[ii,jj] + yGrid[ii,jj]))
+        ChipI = I2[ChipRangeY,ChipRangeX]
+        
+        SearchRangeX = slice(int(-clx - SearchLimitX[ii,jj] + xGrid[ii,jj]) , int(clx + SearchLimitX[ii,jj] - 1 + xGrid[ii,jj]))
+        SearchRangeY = slice(int(-cly - SearchLimitY[ii,jj] + yGrid[ii,jj]) , int(cly + SearchLimitY[ii,jj] - 1 + yGrid[ii,jj]))
+        RefI = I1[SearchRangeY,SearchRangeX]
+        
+        minChipI = np.min(ChipI)
+        if minChipI < 0:
+            ChipI = ChipI - minChipI
+        if np.all(ChipI == ChipI[0,0]):
+            continue
+        
+        minRefI = np.min(RefI)
+        if minRefI < 0:
+            RefI = RefI - minRefI
+        if np.all(RefI == RefI[0,0]):
+            continue
+        
+        
+        if SubPixFlag:
+            # call C++
+            Dx[ind1], Dy[ind1] = np.float32(autoriftcore.arSubPixDisp_s_Py(core._autoriftcore,ChipI.shape[1],ChipI.shape[0],ChipI.ravel(),RefI.shape[1],RefI.shape[0],RefI.ravel(),oversample))
+#                   # call Python
+#                   Dx1[ii], Dy1[ii] = arSubPixDisp(ChipI,RefI)
+        else:
+            # call C++
+            Dx[ind1], Dy[ind1] = np.float32(autoriftcore.arPixDisp_s_Py(core._autoriftcore,ChipI.shape[1],ChipI.shape[0],ChipI.ravel(),RefI.shape[1],RefI.shape[0],RefI.ravel()))
+#                   # call Python
+#                   Dx1[ii], Dy1[ii] = arPixDisp(ChipI,RefI)
+    return Dx, Dy
+
+
+
+
+def arImgDisp_u(I1, I2, xGrid, yGrid, ChipSizeX, ChipSizeY, SearchLimitX, SearchLimitY, Dx0, Dy0, SubPixFlag, oversample, MultiThread):
+
+    import numpy as np
+    from . import autoriftcore
+    import multiprocessing as mp
+    
+    core = AUTO_RIFT_CORE()
+    if core._autoriftcore is not None:
+        autoriftcore.destroyAutoRiftCore_Py(core._autoriftcore)
+
+    core._autoriftcore = autoriftcore.createAutoRiftCore_Py()
+    
     
 #    if np.size(I1) == 1:
 #        if np.logical_not(isinstance(I1,np.uint8) & isinstance(I2,np.uint8)):
@@ -756,50 +940,136 @@ def arImgDisp_u(I1, I2, xGrid, yGrid, ChipSizeX, ChipSizeY, SearchLimitX, Search
     Dx.fill(np.nan)
     Dy = Dx.copy()
 
-
+    if MultiThread == 0:
-    for jj in range(xGrid.shape[1]):
+        for jj in range(xGrid.shape[1]):
-        if np.all(SearchLimitX[:,jj] == 0) & np.all(SearchLimitY[:,jj] == 0):
+            if np.all(SearchLimitX[:,jj] == 0) & np.all(SearchLimitY[:,jj] == 0):
-            continue
-        Dx1 = Dx[:,jj]
-        Dy1 = Dy[:,jj]
-        for ii in range(xGrid.shape[0]):
-            if (SearchLimitX[ii,jj] == 0) & (SearchLimitY[ii,jj] == 0):
                 continue
+            Dx1 = Dx[:,jj]
+            Dy1 = Dy[:,jj]
+            for ii in range(xGrid.shape[0]):
+                if (SearchLimitX[ii,jj] == 0) & (SearchLimitY[ii,jj] == 0):
+                    continue
                 
-            # remember motion terms Dx and Dy correspond to I1 relative to I2 (reference)
+                # remember motion terms Dx and Dy correspond to I1 relative to I2 (reference)
-            clx = np.floor(ChipSizeX[ii,jj]/2)
+                clx = np.floor(ChipSizeX[ii,jj]/2)
-            ChipRangeX = slice(int(-clx - Dx0[ii,jj] + xGrid[ii,jj]) , int(clx - Dx0[ii,jj] + xGrid[ii,jj]))
+                ChipRangeX = slice(int(-clx - Dx0[ii,jj] + xGrid[ii,jj]) , int(clx - Dx0[ii,jj] + xGrid[ii,jj]))
-            cly = np.floor(ChipSizeY[ii,jj]/2)
+                cly = np.floor(ChipSizeY[ii,jj]/2)
-            ChipRangeY = slice(int(-cly - Dy0[ii,jj] + yGrid[ii,jj]) , int(cly - Dy0[ii,jj] + yGrid[ii,jj]))
+                ChipRangeY = slice(int(-cly - Dy0[ii,jj] + yGrid[ii,jj]) , int(cly - Dy0[ii,jj] + yGrid[ii,jj]))
-            ChipI = I2[ChipRangeY,ChipRangeX]
+                ChipI = I2[ChipRangeY,ChipRangeX]
 
-            SearchRangeX = slice(int(-clx - SearchLimitX[ii,jj] + xGrid[ii,jj]) , int(clx + SearchLimitX[ii,jj] - 1 + xGrid[ii,jj]))
+                SearchRangeX = slice(int(-clx - SearchLimitX[ii,jj] + xGrid[ii,jj]) , int(clx + SearchLimitX[ii,jj] - 1 + xGrid[ii,jj]))
-            SearchRangeY = slice(int(-cly - SearchLimitY[ii,jj] + yGrid[ii,jj]) , int(cly + SearchLimitY[ii,jj] - 1 + yGrid[ii,jj]))
+                SearchRangeY = slice(int(-cly - SearchLimitY[ii,jj] + yGrid[ii,jj]) , int(cly + SearchLimitY[ii,jj] - 1 + yGrid[ii,jj]))
-            RefI = I1[SearchRangeY,SearchRangeX]
+                RefI = I1[SearchRangeY,SearchRangeX]
                 
-            minChipI = np.min(ChipI)
+                minChipI = np.min(ChipI)
-            if minChipI < 0:
+                if minChipI < 0:
-                ChipI = ChipI - minChipI
+                    ChipI = ChipI - minChipI
-            if np.all(ChipI == ChipI[0,0]):
+                if np.all(ChipI == ChipI[0,0]):
-                continue
+                    continue
                 
-            minRefI = np.min(RefI)
+                minRefI = np.min(RefI)
-            if minRefI < 0:
+                if minRefI < 0:
-                RefI = RefI - minRefI
+                    RefI = RefI - minRefI
-            if np.all(RefI == RefI[0,0]):
+                if np.all(RefI == RefI[0,0]):
-                continue
+                    continue
                 
 
-            if SubPixFlag:
+                if SubPixFlag:
-                # call C++
+                    # call C++
-                Dx1[ii], Dy1[ii] = np.float32(autoriftcore.arSubPixDisp_u_Py(core._autoriftcore,ChipI.shape[1],ChipI.shape[0],ChipI.ravel(),RefI.shape[1],RefI.shape[0],RefI.ravel(),oversample))
+                    Dx1[ii], Dy1[ii] = np.float32(autoriftcore.arSubPixDisp_u_Py(core._autoriftcore,ChipI.shape[1],ChipI.shape[0],ChipI.ravel(),RefI.shape[1],RefI.shape[0],RefI.ravel(),oversample))
 #                   # call Python
 #                   Dx1[ii], Dy1[ii] = arSubPixDisp(ChipI,RefI)
-            else:
+                else:
-                # call C++
+                    # call C++
-                Dx1[ii], Dy1[ii] = np.float32(autoriftcore.arPixDisp_u_Py(core._autoriftcore,ChipI.shape[1],ChipI.shape[0],ChipI.ravel(),RefI.shape[1],RefI.shape[0],RefI.ravel()))
+                    Dx1[ii], Dy1[ii] = np.float32(autoriftcore.arPixDisp_u_Py(core._autoriftcore,ChipI.shape[1],ChipI.shape[0],ChipI.ravel(),RefI.shape[1],RefI.shape[0],RefI.ravel()))
 #                   # call Python
 #                   Dx1[ii], Dy1[ii] = arPixDisp(ChipI,RefI)
+    else:
+        #   Preparation for parallel
+        in_shape = xGrid.shape
+        I_shape = I1.shape
+        shape_prod = np.asscalar(np.prod(in_shape))
+        
+        #        import pdb
+        #        pdb.set_trace()
+        XI1 = mp.RawArray('b', np.asscalar(np.prod(I_shape)))
+        XI1_np = np.frombuffer(XI1,dtype=np.uint8).reshape(I_shape)
+        np.copyto(XI1_np,I1)
+        del I1
+        
+        XI2 = mp.RawArray('b', np.asscalar(np.prod(I_shape)))
+        XI2_np = np.frombuffer(XI2,dtype=np.uint8).reshape(I_shape)
+        np.copyto(XI2_np,I2)
+        del I2
+        
+        XxGrid = mp.RawArray('f', shape_prod)
+        XxGrid_np = np.frombuffer(XxGrid,dtype=np.float32).reshape(in_shape)
+        np.copyto(XxGrid_np,xGrid)
+        del xGrid
+        
+        XyGrid = mp.RawArray('f', shape_prod)
+        XyGrid_np = np.frombuffer(XyGrid,dtype=np.float32).reshape(in_shape)
+        np.copyto(XyGrid_np,yGrid)
+        del yGrid
+        
+        XSearchLimitX = mp.RawArray('f', shape_prod)
+        XSearchLimitX_np = np.frombuffer(XSearchLimitX,dtype=np.float32).reshape(in_shape)
+        np.copyto(XSearchLimitX_np,SearchLimitX)
+        
+        XSearchLimitY = mp.RawArray('f', shape_prod)
+        XSearchLimitY_np = np.frombuffer(XSearchLimitY,dtype=np.float32).reshape(in_shape)
+        np.copyto(XSearchLimitY_np,SearchLimitY)
+        
+        XChipSizeX = mp.RawArray('f', shape_prod)
+        XChipSizeX_np = np.frombuffer(XChipSizeX,dtype=np.float32).reshape(in_shape)
+        np.copyto(XChipSizeX_np,ChipSizeX)
+        del ChipSizeX
+        
+        XChipSizeY = mp.RawArray('f', shape_prod)
+        XChipSizeY_np = np.frombuffer(XChipSizeY,dtype=np.float32).reshape(in_shape)
+        np.copyto(XChipSizeY_np,ChipSizeY)
+        del ChipSizeY
+        
+        XDx0 = mp.RawArray('f', shape_prod)
+        XDx0_np = np.frombuffer(XDx0,dtype=np.float32).reshape(in_shape)
+        np.copyto(XDx0_np,Dx0)
+        
+        XDy0 = mp.RawArray('f', shape_prod)
+        XDy0_np = np.frombuffer(XDy0,dtype=np.float32).reshape(in_shape)
+        np.copyto(XDy0_np,Dy0)
+        #        import pdb
+        #        pdb.set_trace()
+        
+
+#        Nchunks = mp.cpu_count() // 8 * MultiThread
+        Nchunks = MultiThread
+        chunkSize = int(np.floor(shape_prod / Nchunks))
+        chunkRem = shape_prod - chunkSize * Nchunks
+
+        CHUNKS = []
+
+        for k in range(Nchunks):
+#            print(k)
+            if k == (Nchunks-1):
+                chunkInds = np.arange(k*chunkSize, (k+1)*chunkSize+chunkRem)
+            else:
+                chunkInds = np.arange(k*chunkSize, (k+1)*chunkSize)
+            CHUNKS.append(chunkInds)
+#            print(CHUNKS)
+
+        chunk_inputs = [(kk, CHUNKS[kk], SubPixFlag, oversample, in_shape, I_shape)
+                        for kk in range(Nchunks)]
+
+        with mp.Pool(initializer=initializer, initargs=(XI1, XI2, XxGrid, XyGrid, XSearchLimitX, XSearchLimitY, XChipSizeX, XChipSizeY, XDx0, XDy0)) as pool:
+            Dx, Dy = zip(*pool.map(unpacking_loop_u, chunk_inputs))
+
+        Dx = np.concatenate(Dx)
+        Dy = np.concatenate(Dy)
+
+        Dx = np.reshape(Dx, in_shape)
+        Dy = np.reshape(Dy, in_shape)
+
+
 
 
     # add back 1) I1 (RefI) relative to I2 (ChipI) initial offset Dx0 and Dy0, and
@@ -821,7 +1091,7 @@ def arImgDisp_u(I1, I2, xGrid, yGrid, ChipSizeX, ChipSizeY, SearchLimitX, Search
 
 
 
-def arImgDisp_s(I1, I2, xGrid, yGrid, ChipSizeX, ChipSizeY, SearchLimitX, SearchLimitY, Dx0, Dy0, SubPixFlag, oversample):
+def arImgDisp_s(I1, I2, xGrid, yGrid, ChipSizeX, ChipSizeY, SearchLimitX, SearchLimitY, Dx0, Dy0, SubPixFlag, oversample, MultiThread):
     
     import numpy as np
     from . import autoriftcore
@@ -917,51 +1187,134 @@ def arImgDisp_s(I1, I2, xGrid, yGrid, ChipSizeX, ChipSizeY, SearchLimitX, Search
     Dx.fill(np.nan)
     Dy = Dx.copy()
     
-    
+    if MultiThread == 0:
-    for jj in range(xGrid.shape[1]):
+        for jj in range(xGrid.shape[1]):
-        if np.all(SearchLimitX[:,jj] == 0) & np.all(SearchLimitY[:,jj] == 0):
+            if np.all(SearchLimitX[:,jj] == 0) & np.all(SearchLimitY[:,jj] == 0):
-            continue
-        Dx1 = Dx[:,jj]
-        Dy1 = Dy[:,jj]
-        for ii in range(xGrid.shape[0]):
-            if (SearchLimitX[ii,jj] == 0) & (SearchLimitY[ii,jj] == 0):
                 continue
+            Dx1 = Dx[:,jj]
+            Dy1 = Dy[:,jj]
+            for ii in range(xGrid.shape[0]):
+                if (SearchLimitX[ii,jj] == 0) & (SearchLimitY[ii,jj] == 0):
+                    continue
             
-            # remember motion terms Dx and Dy correspond to I1 relative to I2 (reference)
+                # remember motion terms Dx and Dy correspond to I1 relative to I2 (reference)
-            clx = np.floor(ChipSizeX[ii,jj]/2)
+                clx = np.floor(ChipSizeX[ii,jj]/2)
-            ChipRangeX = slice(int(-clx - Dx0[ii,jj] + xGrid[ii,jj]) , int(clx - Dx0[ii,jj] + xGrid[ii,jj]))
+                ChipRangeX = slice(int(-clx - Dx0[ii,jj] + xGrid[ii,jj]) , int(clx - Dx0[ii,jj] + xGrid[ii,jj]))
-            cly = np.floor(ChipSizeY[ii,jj]/2)
+                cly = np.floor(ChipSizeY[ii,jj]/2)
-            ChipRangeY = slice(int(-cly - Dy0[ii,jj] + yGrid[ii,jj]) , int(cly - Dy0[ii,jj] + yGrid[ii,jj]))
+                ChipRangeY = slice(int(-cly - Dy0[ii,jj] + yGrid[ii,jj]) , int(cly - Dy0[ii,jj] + yGrid[ii,jj]))
-            ChipI = I2[ChipRangeY,ChipRangeX]
+                ChipI = I2[ChipRangeY,ChipRangeX]
                 
-            SearchRangeX = slice(int(-clx - SearchLimitX[ii,jj] + xGrid[ii,jj]) , int(clx + SearchLimitX[ii,jj] - 1 + xGrid[ii,jj]))
+                SearchRangeX = slice(int(-clx - SearchLimitX[ii,jj] + xGrid[ii,jj]) , int(clx + SearchLimitX[ii,jj] - 1 + xGrid[ii,jj]))
-            SearchRangeY = slice(int(-cly - SearchLimitY[ii,jj] + yGrid[ii,jj]) , int(cly + SearchLimitY[ii,jj] - 1 + yGrid[ii,jj]))
+                SearchRangeY = slice(int(-cly - SearchLimitY[ii,jj] + yGrid[ii,jj]) , int(cly + SearchLimitY[ii,jj] - 1 + yGrid[ii,jj]))
-            RefI = I1[SearchRangeY,SearchRangeX]
+                RefI = I1[SearchRangeY,SearchRangeX]
                 
-            minChipI = np.min(ChipI)
+                minChipI = np.min(ChipI)
-            if minChipI < 0:
+                if minChipI < 0:
-                ChipI = ChipI - minChipI
+                    ChipI = ChipI - minChipI
-            if np.all(ChipI == ChipI[0,0]):
+                if np.all(ChipI == ChipI[0,0]):
-                continue
+                    continue
                 
-            minRefI = np.min(RefI)
+                minRefI = np.min(RefI)
-            if minRefI < 0:
+                if minRefI < 0:
-                RefI = RefI - minRefI
+                    RefI = RefI - minRefI
-            if np.all(RefI == RefI[0,0]):
+                if np.all(RefI == RefI[0,0]):
-                continue
+                    continue
         
             
-            if SubPixFlag:
+                if SubPixFlag:
-                # call C++
+                    # call C++
-                Dx1[ii], Dy1[ii] = np.float32(autoriftcore.arSubPixDisp_s_Py(core._autoriftcore,ChipI.shape[1],ChipI.shape[0],ChipI.ravel(),RefI.shape[1],RefI.shape[0],RefI.ravel(), oversample))
+                    Dx1[ii], Dy1[ii] = np.float32(autoriftcore.arSubPixDisp_s_Py(core._autoriftcore,ChipI.shape[1],ChipI.shape[0],ChipI.ravel(),RefI.shape[1],RefI.shape[0],RefI.ravel(),oversample))
 #                   # call Python
 #                   Dx1[ii], Dy1[ii] = arSubPixDisp(ChipI,RefI)
-            else:
+                else:
-                # call C++
+                    # call C++
-                Dx1[ii], Dy1[ii] = np.float32(autoriftcore.arPixDisp_s_Py(core._autoriftcore,ChipI.shape[1],ChipI.shape[0],ChipI.ravel(),RefI.shape[1],RefI.shape[0],RefI.ravel()))
+                    Dx1[ii], Dy1[ii] = np.float32(autoriftcore.arPixDisp_s_Py(core._autoriftcore,ChipI.shape[1],ChipI.shape[0],ChipI.ravel(),RefI.shape[1],RefI.shape[0],RefI.ravel()))
 #                   # call Python
 #                   Dx1[ii], Dy1[ii] = arPixDisp(ChipI,RefI)
+    else:
+        #   Preparation for parallel
+        in_shape = xGrid.shape
+        I_shape = I1.shape
+        shape_prod = np.asscalar(np.prod(in_shape))
         
+        #        import pdb
+        #        pdb.set_trace()
+        XI1 = mp.RawArray('f', np.asscalar(np.prod(I_shape)))
+        XI1_np = np.frombuffer(XI1,dtype=np.float32).reshape(I_shape)
+        np.copyto(XI1_np,I1)
+        del I1
+        
+        XI2 = mp.RawArray('f', np.asscalar(np.prod(I_shape)))
+        XI2_np = np.frombuffer(XI2,dtype=np.float32).reshape(I_shape)
+        np.copyto(XI2_np,I2)
+        del I2
+        
+        XxGrid = mp.RawArray('f', shape_prod)
+        XxGrid_np = np.frombuffer(XxGrid,dtype=np.float32).reshape(in_shape)
+        np.copyto(XxGrid_np,xGrid)
+        del xGrid
+        
+        XyGrid = mp.RawArray('f', shape_prod)
+        XyGrid_np = np.frombuffer(XyGrid,dtype=np.float32).reshape(in_shape)
+        np.copyto(XyGrid_np,yGrid)
+        del yGrid
+        
+        XSearchLimitX = mp.RawArray('f', shape_prod)
+        XSearchLimitX_np = np.frombuffer(XSearchLimitX,dtype=np.float32).reshape(in_shape)
+        np.copyto(XSearchLimitX_np,SearchLimitX)
+        
+        XSearchLimitY = mp.RawArray('f', shape_prod)
+        XSearchLimitY_np = np.frombuffer(XSearchLimitY,dtype=np.float32).reshape(in_shape)
+        np.copyto(XSearchLimitY_np,SearchLimitY)
+        
+        XChipSizeX = mp.RawArray('f', shape_prod)
+        XChipSizeX_np = np.frombuffer(XChipSizeX,dtype=np.float32).reshape(in_shape)
+        np.copyto(XChipSizeX_np,ChipSizeX)
+        del ChipSizeX
+        
+        XChipSizeY = mp.RawArray('f', shape_prod)
+        XChipSizeY_np = np.frombuffer(XChipSizeY,dtype=np.float32).reshape(in_shape)
+        np.copyto(XChipSizeY_np,ChipSizeY)
+        del ChipSizeY
+        
+        XDx0 = mp.RawArray('f', shape_prod)
+        XDx0_np = np.frombuffer(XDx0,dtype=np.float32).reshape(in_shape)
+        np.copyto(XDx0_np,Dx0)
+        
+        XDy0 = mp.RawArray('f', shape_prod)
+        XDy0_np = np.frombuffer(XDy0,dtype=np.float32).reshape(in_shape)
+        np.copyto(XDy0_np,Dy0)
+        #        import pdb
+        #        pdb.set_trace()
+        
+        
+#        Nchunks = mp.cpu_count() // 8 * MultiThread
+        Nchunks = MultiThread
+        chunkSize = int(np.floor(shape_prod / Nchunks))
+        chunkRem = shape_prod - chunkSize * Nchunks
+        
+        CHUNKS = []
+        
+        for k in range(Nchunks):
+            #            print(k)
+            if k == (Nchunks-1):
+                chunkInds = np.arange(k*chunkSize, (k+1)*chunkSize+chunkRem)
+            else:
+                chunkInds = np.arange(k*chunkSize, (k+1)*chunkSize)
+            CHUNKS.append(chunkInds)
+        #            print(CHUNKS)
+        
+        chunk_inputs = [(kk, CHUNKS[kk], SubPixFlag, oversample, in_shape, I_shape)
+                        for kk in range(Nchunks)]
+            
+        with mp.Pool(initializer=initializer, initargs=(XI1, XI2, XxGrid, XyGrid, XSearchLimitX, XSearchLimitY, XChipSizeX, XChipSizeY, XDx0, XDy0)) as pool:
+            Dx, Dy = zip(*pool.map(unpacking_loop_s, chunk_inputs))
+                
+        Dx = np.concatenate(Dx)
+        Dy = np.concatenate(Dy)
+        
+        Dx = np.reshape(Dx, in_shape)
+        Dy = np.reshape(Dy, in_shape)
 
     # add back 1) I1 (RefI) relative to I2 (ChipI) initial offset Dx0 and Dy0, and
     #          2) RefI relative to ChipI has a left/top boundary offset of -SearchLimitX and -SearchLimitY

contrib/geo_autoRIFT/autoRIFT/autoRIFT_ISCE.py

@@ -186,6 +186,13 @@ DATA_TYPE = Component.Parameter('DataType',
          mandatory=False,
          doc = 'Input data type: 0 -> uint8, 1 -> float32')
 
+MULTI_THREAD = Component.Parameter('MultiThread',
+         public_name = 'MULTI_THREAD',
+         default = 0,
+         type = int,
+         mandatory=False,
+         doc = 'Number of Threads; default specified by 0 uses single core and surpasses the multithreading routine')
+
 
 try:
     # Try Autorift within ISCE first
@@ -223,7 +230,8 @@ class autoRIFT_ISCE(autoRIFT, Component):
                       BUFF_DISTANCE_C,
                       COARSE_COR_CUTOFF,
                       OVER_SAMPLE_RATIO,
-                      DATA_TYPE)
+                      DATA_TYPE,
+                      MULTI_THREAD)
     
     
     

contrib/geo_autoRIFT/geogrid/Geogrid.py

@@ -77,7 +77,10 @@ class Geogrid(Component):
             raise Exception('At least the DEM parameter must be set for geogrid')
 
         from osgeo import gdal, osr
-        ds = gdal.Open(self.demname, gdal.GA_ReadOnly)
+        if self.urlflag == 1:
+            ds = gdal.Open('/vsicurl/%s' %(self.demname))
+        else:
+            ds = gdal.Open(self.demname, gdal.GA_ReadOnly)
         srs = osr.SpatialReference()
         srs.ImportFromWkt(ds.GetProjection())
         srs.AutoIdentifyEPSG()
@@ -93,7 +96,10 @@ class Geogrid(Component):
         else:
             raise Exception('Non-standard coordinate system encountered')
         if not epsgstr:  #Empty string->use shell command gdalsrsinfo for last trial
-            cmd = 'gdalsrsinfo -o epsg {0}'.format(self.demname)
+            if self.urlflag == 1:
+                cmd = 'gdalsrsinfo -o epsg /vsicurl/{0}'.format(self.demname)
+            else:
+                cmd = 'gdalsrsinfo -o epsg {0}'.format(self.demname)
             epsgstr = subprocess.check_output(cmd, shell=True)
 #            pdb.set_trace()
             epsgstr = re.findall("EPSG:(\d+)", str(epsgstr))[0]
@@ -274,6 +280,9 @@ class Geogrid(Component):
         geogrid.setRO2VYFilename_Py( self._geogrid, self.winro2vyname)
         geogrid.setLookSide_Py(self._geogrid, self.lookSide)
         geogrid.setNodataOut_Py(self._geogrid, self.nodata_out)
+        if self.urlflag is None:
+            self.urlflag = 0
+        geogrid.setUrlFlag_Py(self._geogrid, self.urlflag)
 
         self._orbit  = self.orbit.exportToC()
         geogrid.setOrbit_Py(self._geogrid, self._orbit)
@@ -321,6 +330,7 @@ class Geogrid(Component):
         self.csmaxxname = None
         self.csmaxyname = None
         self.ssmname = None
+        self.urlflag = None
         
         ##Output related parameters
         self.winlocname = None

contrib/geo_autoRIFT/geogrid/GeogridOptical.py

@@ -48,8 +48,8 @@ class GeogridOptical():
         '''
 
         ##Determine appropriate EPSG system
-        self.epsgDem = self.getProjectionSystem(self.demname)
+        self.epsgDem = self.getProjectionSystem(self.demname, self.urlflag)
-        self.epsgDat = self.getProjectionSystem(self.dat1name)
+        self.epsgDat = self.getProjectionSystem(self.dat1name, self.urlflag)
         
         ###Determine extent of data needed
         bbox = self.determineBbox()
@@ -59,7 +59,7 @@ class GeogridOptical():
         self.geogrid()
 
 
-    def getProjectionSystem(self, filename):
+    def getProjectionSystem(self, filename, urlflag):
         '''
         Testing with Greenland.
         '''
@@ -67,7 +67,10 @@ class GeogridOptical():
             raise Exception('File {0} does not exist'.format(filename))
 
         from osgeo import gdal, osr
-        ds = gdal.Open(filename, gdal.GA_ReadOnly)
+        if urlflag == 1:
+            ds = gdal.Open('/vsicurl/%s' %(filename))
+        else:
+            ds = gdal.Open(filename, gdal.GA_ReadOnly)
         srs = osr.SpatialReference()
         srs.ImportFromWkt(ds.GetProjection())
         srs.AutoIdentifyEPSG()
@@ -83,7 +86,10 @@ class GeogridOptical():
         else:
             raise Exception('Non-standard coordinate system encountered')
         if not epsgstr:  #Empty string->use shell command gdalsrsinfo for last trial
-            cmd = 'gdalsrsinfo -o epsg {0}'.format(filename)
+            if urlflag == 1:
+                cmd = 'gdalsrsinfo -o epsg /vsicurl/{0}'.format(filename)
+            else:
+                cmd = 'gdalsrsinfo -o epsg {0}'.format(filename)
             epsgstr = subprocess.check_output(cmd, shell=True)
 #            pdb.set_trace()
             epsgstr = re.findall("EPSG:(\d+)", str(epsgstr))[0]
@@ -155,6 +161,13 @@ class GeogridOptical():
         
         #   For now print inputs that were obtained
         
+        urlflag = self.urlflag
+        
+        if urlflag == 1:
+            print("\nReading input images into memory directly from URL's")
+        else:
+            print("\nReading input images locally from files")
+    
         print("\nOptical Image parameters: ")
         print("X-direction coordinate: " + str(self.startingX) + "  " + str(self.XSize))
         print("Y-direction coordinate: " + str(self.startingY) + "  " + str(self.YSize))
@@ -218,6 +231,21 @@ class GeogridOptical():
         import struct
         
 #        pdb.set_trace()
+        if urlflag == 1:
+            self.demname = '/vsicurl/%s' %(self.demname)
+            self.dhdxname = '/vsicurl/%s' %(self.dhdxname)
+            self.dhdyname = '/vsicurl/%s' %(self.dhdyname)
+            self.vxname = '/vsicurl/%s' %(self.vxname)
+            self.vyname = '/vsicurl/%s' %(self.vyname)
+            self.srxname = '/vsicurl/%s' %(self.srxname)
+            self.sryname = '/vsicurl/%s' %(self.sryname)
+            self.csminxname = '/vsicurl/%s' %(self.csminxname)
+            self.csminyname = '/vsicurl/%s' %(self.csminyname)
+            self.csmaxxname = '/vsicurl/%s' %(self.csmaxxname)
+            self.csmaxyname = '/vsicurl/%s' %(self.csmaxyname)
+            self.ssmname = '/vsicurl/%s' %(self.ssmname)
+        
+
         demDS = gdal.Open(self.demname, gdal.GA_ReadOnly)
         
         if (self.dhdxname != ""):
@@ -731,8 +759,94 @@ class GeogridOptical():
 
 
 
+    def coregister(self,in1,in2,urlflag):
+        import os
+        import numpy as np
         
+        from osgeo import gdal, osr
+        import struct
         
+        if urlflag == 1:
+            DS1 = gdal.Open('/vsicurl/%s' %(in1))
+        else:
+            DS1 = gdal.Open(in1, gdal.GA_ReadOnly)
+        trans1 = DS1.GetGeoTransform()
+        xsize1 = DS1.RasterXSize
+        ysize1 = DS1.RasterYSize
+        
+        if urlflag == 1:
+            DS2 = gdal.Open('/vsicurl/%s' %(in2))
+        else:
+            DS2 = gdal.Open(in2, gdal.GA_ReadOnly)
+        trans2 = DS2.GetGeoTransform()
+        xsize2 = DS2.RasterXSize
+        ysize2 = DS2.RasterYSize
+        
+        W = np.max([trans1[0],trans2[0]])
+        N = np.min([trans1[3],trans2[3]])
+        E = np.min([trans1[0]+xsize1*trans1[1],trans2[0]+xsize2*trans2[1]])
+        S = np.max([trans1[3]+ysize1*trans1[5],trans2[3]+ysize2*trans2[5]])
+        
+        x1a = int(np.round((W-trans1[0])/trans1[1]))
+        x1b = int(np.round((E-trans1[0])/trans1[1]))
+        y1a = int(np.round((N-trans1[3])/trans1[5]))
+        y1b = int(np.round((S-trans1[3])/trans1[5]))
+        
+        x2a = int(np.round((W-trans2[0])/trans2[1]))
+        x2b = int(np.round((E-trans2[0])/trans2[1]))
+        y2a = int(np.round((N-trans2[3])/trans2[5]))
+        y2b = int(np.round((S-trans2[3])/trans2[5]))
+        
+        x1a = np.min([x1a, xsize1-1])
+        x1b = np.min([x1b, xsize1-1])
+        y1a = np.min([y1a, ysize1-1])
+        y1b = np.min([y1b, ysize1-1])
+        x2a = np.min([x2a, xsize2-1])
+        x2b = np.min([x2b, xsize2-1])
+        y2a = np.min([y2a, ysize2-1])
+        y2b = np.min([y2b, ysize2-1])
+        
+        x1a = np.max([x1a, 0])
+        x1b = np.max([x1b, 0])
+        y1a = np.max([y1a, 0])
+        y1b = np.max([y1b, 0])
+        x2a = np.max([x2a, 0])
+        x2b = np.max([x2b, 0])
+        y2a = np.max([y2a, 0])
+        y2b = np.max([y2b, 0])
+        
+        x1a = int(x1a)
+        x1b = int(x1b)
+        y1a = int(y1a)
+        y1b = int(y1b)
+        x2a = int(x2a)
+        x2b = int(x2b)
+        y2a = int(y2a)
+        y2b = int(y2b)
+
+        trans = (W, trans1[1], 0.0, N, 0.0, trans1[5])
+
+        if urlflag == 0:
+            
+            I1 = DS1.ReadAsArray(xoff=x1a, yoff=y1a, xsize=x1b-x1a+1, ysize=y1b-y1a+1)
+            I2 = DS2.ReadAsArray(xoff=x2a, yoff=y2a, xsize=x2b-x2a+1, ysize=y2b-y2a+1)
+
+            fileformat = "GTiff"
+            driver = gdal.GetDriverByName(fileformat)
+            
+            DST1 = driver.Create(os.path.basename(in1), xsize=(x1b-x1a+1), ysize=(y1b-y1a+1), bands=1, eType=gdal.GDT_UInt16)
+            DST1.SetGeoTransform(trans)
+            DST1.SetProjection(DS1.GetProjectionRef())
+            DST1.GetRasterBand(1).WriteArray(I1)
+            DST1 = None
+            
+            DST2 = driver.Create(os.path.basename(in2), xsize=(x2b-x2a+1), ysize=(y2b-y2a+1), bands=1, eType=gdal.GDT_UInt16)
+            DST2.SetGeoTransform(trans)
+            DST2.SetProjection(DS2.GetProjectionRef())
+            DST2.GetRasterBand(1).WriteArray(I2)
+            DST2 = None
+        
+        return x1a, y1a, x1b-x1a+1, y1b-y1a+1, x2a, y2a, x2b-x2a+1, y2b-y2a+1, trans
 
 
 
@@ -753,6 +867,7 @@ class GeogridOptical():
         self.numberOfLines = None
         self.repeatTime = None
         self.chipSizeX0 = None
+        self.urlflag = None
 
         ##Input related parameters
         self.dat1name = None

contrib/geo_autoRIFT/geogrid/bindings/geogridmodule.cpp

@@ -430,6 +430,19 @@ PyObject* setNodataOut(PyObject *self, PyObject *args)
     return Py_BuildValue("i", 0);
 }
 
+PyObject* setUrlFlag(PyObject *self, PyObject *args)
+{
+    uint64_t ptr;
+    int urlflag;
+    if (!PyArg_ParseTuple(args, "Ki", &ptr, &urlflag))
+    {
+        return NULL;
+    }
+    
+    ((geoGrid*)(ptr))->urlflag = urlflag;
+    return Py_BuildValue("i", 0);
+}
+
 PyObject* setOrbit(PyObject *self, PyObject *args)
 {
     uint64_t ptr, orb;

contrib/geo_autoRIFT/geogrid/include/geogrid.h

@@ -69,6 +69,7 @@ struct geoGrid
     int nPixels;
     int lookSide;
     int nodata_out;
+    int urlflag;
     double incidenceAngle;
 
     //Output file names

contrib/geo_autoRIFT/geogrid/include/geogridmodule.h

@@ -53,6 +53,7 @@ extern "C"
         PyObject * setOrbit(PyObject *, PyObject *);
         PyObject * setLookSide(PyObject *, PyObject *);
         PyObject * setNodataOut(PyObject *, PyObject *);
+        PyObject * setUrlFlag(PyObject *, PyObject *);
 
         PyObject * setWindowLocationsFilename(PyObject *, PyObject *);
         PyObject * setWindowOffsetsFilename(PyObject *, PyObject *);
@@ -91,6 +92,7 @@ static PyMethodDef geogrid_methods[] =
         {"setOrbit_Py", setOrbit, METH_VARARGS, " "},
         {"setLookSide_Py", setLookSide, METH_VARARGS, " "},
         {"setNodataOut_Py", setNodataOut, METH_VARARGS, " "},
+        {"setUrlFlag_Py", setUrlFlag, METH_VARARGS, " "},
         {"setXLimits_Py", setXLimits, METH_VARARGS, " "},
         {"setYLimits_Py", setYLimits, METH_VARARGS, " "},
         {"setWindowLocationsFilename_Py", setWindowLocationsFilename, METH_VARARGS, " "},

contrib/geo_autoRIFT/geogrid/src/geogrid.cpp

@@ -47,6 +47,14 @@ void geoGrid::geogrid()
 
     //For now print inputs that were obtained
     
+    if (urlflag == 1){
+        std::cout << "\nReading input images into memory directly from URL's\n";
+    }
+    else
+    {
+        std::cout << "\nReading input images locally from files\n";
+    }
+    
     std::cout << "\nRadar parameters: \n";
     std::cout << "Range: " << startingRange << "  " << dr << "\n";
     std::cout << "Azimuth: " << sensingStart << "  " << prf << "\n";
@@ -145,6 +153,22 @@ void geoGrid::geogrid()
 
     double geoTrans[6];
     
+    std::string url ("/vsicurl/");
+    if (urlflag == 1)
+    {
+        demname = url + demname;
+        dhdxname = url + dhdxname;
+        dhdyname = url + dhdyname;
+        vxname = url + vxname;
+        vyname = url + vyname;
+        srxname = url + srxname;
+        sryname = url + sryname;
+        csminxname = url + csminxname;
+        csminyname = url + csminyname;
+        csmaxxname = url + csmaxxname;
+        csmaxyname = url + csmaxyname;
+        ssmname = url + ssmname;
+    }
     demDS = reinterpret_cast<GDALDataset *>(GDALOpenShared(demname.c_str(), GA_ReadOnly));
     if (dhdxname != "")
     {