From 3673da71ccf18a284497bcb005a3013835cd67fe Mon Sep 17 00:00:00 2001 From: Ryan Burns Date: Tue, 2 Feb 2021 16:29:29 -0800 Subject: [PATCH] Bump autoRIFT 1.0.7 -> 1.1.0 --- contrib/geo_autoRIFT/README.md | 7 +- contrib/geo_autoRIFT/autoRIFT/autoRIFT.py | 557 ++++++++++++++---- .../geo_autoRIFT/autoRIFT/autoRIFT_ISCE.py | 10 +- contrib/geo_autoRIFT/geogrid/Geogrid.py | 14 +- .../geo_autoRIFT/geogrid/GeogridOptical.py | 125 +++- .../geogrid/bindings/geogridmodule.cpp | 13 + .../geo_autoRIFT/geogrid/include/geogrid.h | 1 + .../geogrid/include/geogridmodule.h | 2 + contrib/geo_autoRIFT/geogrid/src/geogrid.cpp | 26 +- 9 files changed, 641 insertions(+), 114 deletions(-) diff --git a/contrib/geo_autoRIFT/README.md b/contrib/geo_autoRIFT/README.md index cf52e4b..4051bf6 100644 --- a/contrib/geo_autoRIFT/README.md +++ b/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 diff --git a/contrib/geo_autoRIFT/autoRIFT/autoRIFT.py b/contrib/geo_autoRIFT/autoRIFT/autoRIFT.py index 976a29b..78d3ff1 100755 --- a/contrib/geo_autoRIFT/autoRIFT/autoRIFT.py +++ b/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,9 +158,7 @@ class autoRIFT: - if self.zeroMask is not None: - self.zeroMask = (self.I1 == 0) - + sobelx = cv2.getDerivKernels(1,0,self.WallisFilterWidth) kernelx = np.outer(sobelx[0],sobelx[1]) @@ -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)) - validData = np.isfinite(self.I1) - S1 = np.std(self.I1[validData])*np.sqrt(self.I1[validData].size/(self.I1[validData].size-1.0)) - M1 = np.mean(self.I1[validData]) + if self.zeroMask is not None: + # validData = np.logical_not(np.isnan(self.I1)) + validData = np.isfinite(self.I1) + 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)) - validData = np.isfinite(self.I2) - S2 = np.std(self.I2[validData])*np.sqrt(self.I2[validData].size/(self.I2[validData].size-1.0)) - M2 = np.mean(self.I2[validData]) + if self.zeroMask is not None: + # validData = np.logical_not(np.isnan(self.I2)) + validData = np.isfinite(self.I2) + 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,11 +666,188 @@ 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 + + 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.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: @@ -756,50 +940,136 @@ def arImgDisp_u(I1, I2, xGrid, yGrid, ChipSizeX, ChipSizeY, SearchLimitX, Search Dx.fill(np.nan) Dy = Dx.copy() - - for jj in range(xGrid.shape[1]): - 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): + if MultiThread == 0: + for jj in range(xGrid.shape[1]): + if np.all(SearchLimitX[:,jj] == 0) & np.all(SearchLimitY[:,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] + 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) + 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 - + 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++ - 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)) + if SubPixFlag: + # 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)) # # call Python # Dx1[ii], Dy1[ii] = arSubPixDisp(ChipI,RefI) - else: - # 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())) + else: + # 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())) # # 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() - - for jj in range(xGrid.shape[1]): - 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): + if MultiThread == 0: + for jj in range(xGrid.shape[1]): + 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 + + # 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 - # 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++ - 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)) + if SubPixFlag: + # 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)) # # call Python # Dx1[ii], Dy1[ii] = arSubPixDisp(ChipI,RefI) - else: - # 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())) + else: + # 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())) # # 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 diff --git a/contrib/geo_autoRIFT/autoRIFT/autoRIFT_ISCE.py b/contrib/geo_autoRIFT/autoRIFT/autoRIFT_ISCE.py index 7345e59..39bd76b 100755 --- a/contrib/geo_autoRIFT/autoRIFT/autoRIFT_ISCE.py +++ b/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) diff --git a/contrib/geo_autoRIFT/geogrid/Geogrid.py b/contrib/geo_autoRIFT/geogrid/Geogrid.py index 77f9c6a..5e0a96a 100755 --- a/contrib/geo_autoRIFT/geogrid/Geogrid.py +++ b/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 diff --git a/contrib/geo_autoRIFT/geogrid/GeogridOptical.py b/contrib/geo_autoRIFT/geogrid/GeogridOptical.py index 1ce6f88..a120603 100755 --- a/contrib/geo_autoRIFT/geogrid/GeogridOptical.py +++ b/contrib/geo_autoRIFT/geogrid/GeogridOptical.py @@ -48,8 +48,8 @@ class GeogridOptical(): ''' ##Determine appropriate EPSG system - self.epsgDem = self.getProjectionSystem(self.demname) - self.epsgDat = self.getProjectionSystem(self.dat1name) + self.epsgDem = self.getProjectionSystem(self.demname, self.urlflag) + 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] @@ -154,6 +160,13 @@ class GeogridOptical(): def geogrid(self): # 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)) @@ -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 diff --git a/contrib/geo_autoRIFT/geogrid/bindings/geogridmodule.cpp b/contrib/geo_autoRIFT/geogrid/bindings/geogridmodule.cpp index 4a5c3f3..7dd37cc 100755 --- a/contrib/geo_autoRIFT/geogrid/bindings/geogridmodule.cpp +++ b/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; diff --git a/contrib/geo_autoRIFT/geogrid/include/geogrid.h b/contrib/geo_autoRIFT/geogrid/include/geogrid.h index 966df6c..a856c61 100755 --- a/contrib/geo_autoRIFT/geogrid/include/geogrid.h +++ b/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 diff --git a/contrib/geo_autoRIFT/geogrid/include/geogridmodule.h b/contrib/geo_autoRIFT/geogrid/include/geogridmodule.h index 22bb6d0..77c205e 100755 --- a/contrib/geo_autoRIFT/geogrid/include/geogridmodule.h +++ b/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, " "}, diff --git a/contrib/geo_autoRIFT/geogrid/src/geogrid.cpp b/contrib/geo_autoRIFT/geogrid/src/geogrid.cpp index abb2d2f..ea6e07d 100755 --- a/contrib/geo_autoRIFT/geogrid/src/geogrid.cpp +++ b/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"; @@ -144,7 +152,23 @@ void geoGrid::geogrid() GDALDataset* ssmDS = NULL; 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(GDALOpenShared(demname.c_str(), GA_ReadOnly)); if (dhdxname != "") {