From 1557a51ef569f5ce4cddf1b9cf424ab1c0a0409d Mon Sep 17 00:00:00 2001 From: Yang Lei Date: Mon, 24 May 2021 13:51:10 -0700 Subject: [PATCH] upgrade autoRIFT and Geogrid to v1.3.0 for memory and runtime improvement --- geo_autoRIFT/SConscript | 2 +- geo_autoRIFT/autoRIFT/__init__.py | 2 +- geo_autoRIFT/autoRIFT/autoRIFT.py | 167 +- geo_autoRIFT/autoRIFT/autoRIFT_ISCE.py | 8 + geo_autoRIFT/geogrid/GeogridOptical.py | 679 ++------- geo_autoRIFT/geogrid/SConscript | 15 + geo_autoRIFT/geogrid/__init__.py | 2 +- geo_autoRIFT/geogrid/bindings/SConscript1 | 28 + .../geogrid/bindings/geogridOpticalmodule.cpp | 567 +++++++ geo_autoRIFT/geogrid/include/SConscript1 | 24 + geo_autoRIFT/geogrid/include/geogridOptical.h | 93 ++ .../geogrid/include/geogridOpticalmodule.h | 124 ++ geo_autoRIFT/geogrid/src/SConscript1 | 24 + geo_autoRIFT/geogrid/src/geogrid.cpp | 35 +- geo_autoRIFT/geogrid/src/geogridOptical.cpp | 1347 +++++++++++++++++ setup.py | 23 +- testGeogridOptical.py | 2 + testGeogrid_ISCE.py | 2 +- testautoRIFT.py | 37 +- testautoRIFT_ISCE.py | 33 +- 20 files changed, 2528 insertions(+), 686 deletions(-) create mode 100755 geo_autoRIFT/geogrid/bindings/SConscript1 create mode 100755 geo_autoRIFT/geogrid/bindings/geogridOpticalmodule.cpp create mode 100755 geo_autoRIFT/geogrid/include/SConscript1 create mode 100755 geo_autoRIFT/geogrid/include/geogridOptical.h create mode 100755 geo_autoRIFT/geogrid/include/geogridOpticalmodule.h create mode 100755 geo_autoRIFT/geogrid/src/SConscript1 create mode 100755 geo_autoRIFT/geogrid/src/geogridOptical.cpp mode change 100644 => 100755 setup.py diff --git a/geo_autoRIFT/SConscript b/geo_autoRIFT/SConscript index 97d9cad..7235c9b 100755 --- a/geo_autoRIFT/SConscript +++ b/geo_autoRIFT/SConscript @@ -34,4 +34,4 @@ autorift='autoRIFT/SConscript' SConscript(autorift) geogrid='geogrid/SConscript' -SConscript(geogrid) +SConscript(geogrid) \ No newline at end of file diff --git a/geo_autoRIFT/autoRIFT/__init__.py b/geo_autoRIFT/autoRIFT/__init__.py index 7cca7dc..a73bfe6 100755 --- a/geo_autoRIFT/autoRIFT/__init__.py +++ b/geo_autoRIFT/autoRIFT/__init__.py @@ -13,4 +13,4 @@ # this means ISCE support not available. Don't raise error. Allow standalone use pass -__version__ = '1.2.0' +__version__ = '1.3.0' diff --git a/geo_autoRIFT/autoRIFT/autoRIFT.py b/geo_autoRIFT/autoRIFT/autoRIFT.py index 314e454..f67f8f4 100755 --- a/geo_autoRIFT/autoRIFT/autoRIFT.py +++ b/geo_autoRIFT/autoRIFT/autoRIFT.py @@ -105,7 +105,7 @@ def preprocess_filt_wal(self): # # self.I2 = (self.I2 - lp) - + def preprocess_filt_hps(self): ''' Do the pre processing using (orig - low-pass filter) = high-pass filter filter (3.9/5.3 min). @@ -205,7 +205,7 @@ def preprocess_filt_lap(self): - + def uniform_data_type(self): import numpy as np @@ -342,6 +342,7 @@ def autorift(self): DispFiltC.FiltWidth = (self.FiltWidth - 1) * ChipSize0_GridSpacing_oversample_ratio + 1 DispFiltC.Iter = self.Iter - 1 DispFiltC.MadScalar = self.MadScalar + DispFiltC.colfiltChunkSize = self.colfiltChunkSize DispFiltF = DISP_FILT() overlap_f = 1 - 1 / ChipSize0_GridSpacing_oversample_ratio @@ -350,6 +351,7 @@ def autorift(self): DispFiltF.FiltWidth = (self.FiltWidth - 1) * ChipSize0_GridSpacing_oversample_ratio + 1 DispFiltF.Iter = self.Iter DispFiltF.MadScalar = self.MadScalar + DispFiltF.colfiltChunkSize = self.colfiltChunkSize for i in range(ChipSizeUniX.__len__()): @@ -369,13 +371,13 @@ def autorift(self): yGrid0 = np.round(yGrid0) M0 = (ChipSizeX == 0) & (self.ChipSizeMinX <= ChipSizeUniX[i]) & (self.ChipSizeMaxX >= ChipSizeUniX[i]) - M0 = colfilt(M0.copy(), (int(1/Scale*6), int(1/Scale*6)), 0) + M0 = colfilt(M0.copy(), (int(1/Scale*6), int(1/Scale*6)), 0, self.colfiltChunkSize) M0 = cv2.resize(np.logical_not(M0).astype(np.uint8),dstShape[::-1],interpolation=cv2.INTER_NEAREST).astype(np.bool) - SearchLimitX0 = colfilt(self.SearchLimitX.copy(), (int(1/Scale), int(1/Scale)), 0) + colfilt(self.Dx0.copy(), (int(1/Scale), int(1/Scale)), 4) - SearchLimitY0 = colfilt(self.SearchLimitY.copy(), (int(1/Scale), int(1/Scale)), 0) + colfilt(self.Dy0.copy(), (int(1/Scale), int(1/Scale)), 4) - Dx00 = colfilt(self.Dx0.copy(), (int(1/Scale), int(1/Scale)), 2) - Dy00 = colfilt(self.Dy0.copy(), (int(1/Scale), int(1/Scale)), 2) + SearchLimitX0 = colfilt(self.SearchLimitX.copy(), (int(1/Scale), int(1/Scale)), 0, self.colfiltChunkSize) + colfilt(self.Dx0.copy(), (int(1/Scale), int(1/Scale)), 4, self.colfiltChunkSize) + SearchLimitY0 = colfilt(self.SearchLimitY.copy(), (int(1/Scale), int(1/Scale)), 0, self.colfiltChunkSize) + colfilt(self.Dy0.copy(), (int(1/Scale), int(1/Scale)), 4, self.colfiltChunkSize) + Dx00 = colfilt(self.Dx0.copy(), (int(1/Scale), int(1/Scale)), 2, self.colfiltChunkSize) + Dy00 = colfilt(self.Dy0.copy(), (int(1/Scale), int(1/Scale)), 2, self.colfiltChunkSize) SearchLimitX0 = np.ceil(cv2.resize(SearchLimitX0,dstShape[::-1])) SearchLimitY0 = np.ceil(cv2.resize(SearchLimitY0,dstShape[::-1])) @@ -421,8 +423,8 @@ def autorift(self): else: filtWidth = (self.sparseSearchSampleRate * ChipSize0_GridSpacing_oversample_ratio) - SearchLimitX0C = colfilt(SearchLimitX0.copy(), (int(filtWidth), int(filtWidth)), 0) - SearchLimitY0C = colfilt(SearchLimitY0.copy(), (int(filtWidth), int(filtWidth)), 0) + SearchLimitX0C = colfilt(SearchLimitX0.copy(), (int(filtWidth), int(filtWidth)), 0, self.colfiltChunkSize) + SearchLimitY0C = colfilt(SearchLimitY0.copy(), (int(filtWidth), int(filtWidth)), 0, self.colfiltChunkSize) SearchLimitX0C = SearchLimitX0C[rIdxC,cIdxC] SearchLimitY0C = SearchLimitY0C[rIdxC,cIdxC] @@ -499,8 +501,8 @@ def autorift(self): DyF[np.logical_not(M0)] = np.nan # Light interpolation with median filtered values: DxFM (filtered) and DxF (unfiltered) - DxFM = colfilt(DxF.copy(), (self.fillFiltWidth, self.fillFiltWidth), 3) - DyFM = colfilt(DyF.copy(), (self.fillFiltWidth, self.fillFiltWidth), 3) + DxFM = colfilt(DxF.copy(), (self.fillFiltWidth, self.fillFiltWidth), 3, self.colfiltChunkSize) + DyFM = colfilt(DyF.copy(), (self.fillFiltWidth, self.fillFiltWidth), 3, self.colfiltChunkSize) # M0 is mask for original valid estimates, MF is mask for filled ones, MM is mask where filtered ones exist for filling MF = np.zeros(M0.shape, dtype=np.bool) @@ -529,16 +531,16 @@ def autorift(self): # DxF0 (filtered) / Dx (unfiltered) is the result from earlier iterations, DxFM (filtered) / DxF (unfiltered) is that of the current iteration # first colfilt nans within 2-by-2 area (otherwise 1 nan will contaminate all 4 points) - DxF0 = colfilt(Dx.copy(),(int(Scale+1),int(Scale+1)),2) + DxF0 = colfilt(Dx.copy(),(int(Scale+1),int(Scale+1)),2, self.colfiltChunkSize) # then resize to half size using area (similar to averaging) to match the current iteration DxF0 = cv2.resize(DxF0,dstShape[::-1],interpolation=cv2.INTER_AREA) - DyF0 = colfilt(Dy.copy(),(int(Scale+1),int(Scale+1)),2) + DyF0 = colfilt(Dy.copy(),(int(Scale+1),int(Scale+1)),2, self.colfiltChunkSize) DyF0 = cv2.resize(DyF0,dstShape[::-1],interpolation=cv2.INTER_AREA) # Note this DxFM is almost the same as DxFM (same variable) in the light interpolation (only slightly better); however, only small portion of it will be used later at locations specified by M0 and MF that are determined in the light interpolation. So even without the following two lines, the final Dx and Dy result is still the same. # to fill out all of the missing values in DxF - DxFM = colfilt(DxF.copy(), (5,5), 3) - DyFM = colfilt(DyF.copy(), (5,5), 3) + DxFM = colfilt(DxF.copy(), (5,5), 3, self.colfiltChunkSize) + DyFM = colfilt(DyF.copy(), (5,5), 3, self.colfiltChunkSize) # fill the current-iteration result with previously determined reliable estimates that are not searched in the current iteration idx = np.isnan(DxF) & np.logical_not(np.isnan(DxF0)) @@ -578,7 +580,7 @@ def autorift(self): - + def runAutorift(self): ''' quick processing routine which calls autorift main function (user can define their own way by mimicing the workflow here). @@ -661,6 +663,7 @@ def __init__(self): self.FiltWidth = 5 self.Iter = 3 self.MadScalar = 4 + self.colfiltChunkSize = 4 self.BuffDistanceC = 8 self.CoarseCorCutoff = 0.01 self.OverSampleRatio = 16 @@ -1346,44 +1349,93 @@ def arImgDisp_s(I1, I2, xGrid, yGrid, ChipSizeX, ChipSizeY, SearchLimitX, Search -def colfilt(A, kernelSize, option): + +################## Chunked version of column filter +def colfilt(A, kernelSize, option, chunkSize=4): from skimage.util import view_as_windows as viewW import numpy as np - A = np.lib.pad(A,((int((kernelSize[0]-1)/2),int((kernelSize[0]-1)/2)),(int((kernelSize[1]-1)/2),int((kernelSize[1]-1)/2))),mode='constant',constant_values=np.nan) - - B = viewW(A, kernelSize).reshape(-1,kernelSize[0]*kernelSize[1]).T[:,::1] - - output_size = (A.shape[0]-kernelSize[0]+1,A.shape[1]-kernelSize[1]+1) - C = np.zeros(output_size,dtype=A.dtype) - if option == 0:# max - C = np.nanmax(B,axis=0).reshape(output_size) - elif option == 1:# min - C = np.nanmin(B,axis=0).reshape(output_size) - elif option == 2:# mean - C = np.nanmean(B,axis=0).reshape(output_size) - elif option == 3:# median - C = np.nanmedian(B,axis=0).reshape(output_size) - elif option == 4:# range - C = np.nanmax(B,axis=0).reshape(output_size) - np.nanmin(B,axis=0).reshape(output_size) - elif option == 6:# MAD (Median Absolute Deviation) - m = B.shape[0] - D = np.abs(B - np.dot(np.ones((m,1),dtype=A.dtype), np.array([np.nanmedian(B,axis=0)]))) - C = np.nanmedian(D,axis=0).reshape(output_size) - elif option[0] == 5:# displacement distance count with option[1] being the threshold - m = B.shape[0] - c = int(np.round((m + 1) / 2)-1) -# c = 0 - D = np.abs(B - np.dot(np.ones((m,1),dtype=A.dtype), np.array([B[c,:]]))) - C = np.sum(D= dToleranceX) & (colfilt(Dy.copy(), (self.FiltWidth, self.FiltWidth), (5,self.FracSearch)) >= dToleranceY) - + M = (colfilt(Dx.copy(), (self.FiltWidth, self.FiltWidth), (5,self.FracSearch), self.colfiltChunkSize) >= dToleranceX) & (colfilt(Dy.copy(), (self.FiltWidth, self.FiltWidth), (5,self.FracSearch), self.colfiltChunkSize) >= dToleranceY) + # if self.Iter == 3: # pdb.set_trace() @@ -1428,12 +1482,13 @@ def filtDisp(self, Dx, Dy, SearchLimitX, SearchLimitY, M, OverSampleRatio): Dx[np.logical_not(M)] = np.nan Dy[np.logical_not(M)] = np.nan - DxMad = colfilt(Dx.copy(), (self.FiltWidth, self.FiltWidth), 6) - DyMad = colfilt(Dy.copy(), (self.FiltWidth, self.FiltWidth), 6) - - DxM = colfilt(Dx.copy(), (self.FiltWidth, self.FiltWidth), 3) - DyM = colfilt(Dy.copy(), (self.FiltWidth, self.FiltWidth), 3) + DxMad = colfilt(Dx.copy(), (self.FiltWidth, self.FiltWidth), 6, self.colfiltChunkSize) + DyMad = colfilt(Dy.copy(), (self.FiltWidth, self.FiltWidth), 6, self.colfiltChunkSize) + DxM = colfilt(Dx.copy(), (self.FiltWidth, self.FiltWidth), 3, self.colfiltChunkSize) + DyM = colfilt(Dy.copy(), (self.FiltWidth, self.FiltWidth), 3, self.colfiltChunkSize) + + M = (np.abs(Dx - DxM) <= np.maximum(self.MadScalar * DxMad, DxMadmin)) & (np.abs(Dy - DyM) <= np.maximum(self.MadScalar * DyMad, DyMadmin)) & M return M @@ -1458,5 +1513,3 @@ def bwareaopen(image,size1): return image - - diff --git a/geo_autoRIFT/autoRIFT/autoRIFT_ISCE.py b/geo_autoRIFT/autoRIFT/autoRIFT_ISCE.py index bf4941d..4cb64b5 100755 --- a/geo_autoRIFT/autoRIFT/autoRIFT_ISCE.py +++ b/geo_autoRIFT/autoRIFT/autoRIFT_ISCE.py @@ -165,6 +165,13 @@ mandatory=False, doc = 'Mad Scalar') +COLFILT_CHUNK_SIZE = Component.Parameter('colfiltChunkSize', + public_name = 'COLFILT_CHUNK_SIZE', + default = 4, + type = int, + mandatory=False, + doc = 'column filter chunk size') + BUFF_DISTANCE_C = Component.Parameter('BuffDistanceC', public_name = 'BUFF_DISTANCE_C', default = 8, @@ -235,6 +242,7 @@ class autoRIFT_ISCE(autoRIFT, Component): FILT_WIDTH, ITER, MAD_SCALAR, + COLFILT_CHUNK_SIZE, BUFF_DISTANCE_C, COARSE_COR_CUTOFF, OVER_SAMPLE_RATIO, diff --git a/geo_autoRIFT/geogrid/GeogridOptical.py b/geo_autoRIFT/geogrid/GeogridOptical.py index d76edd6..4d69f5e 100755 --- a/geo_autoRIFT/geogrid/GeogridOptical.py +++ b/geo_autoRIFT/geogrid/GeogridOptical.py @@ -46,6 +46,8 @@ def runGeogrid(self): ''' Do the actual processing. ''' + + from . import geogridOptical ##Determine appropriate EPSG system self.epsgDem = self.getProjectionSystem(self.demname) @@ -54,18 +56,28 @@ def runGeogrid(self): ###Determine extent of data needed bbox = self.determineBbox() + ##Create and set parameters + self.setState() + ##check parameters self.checkState() ##Run - self.geogrid() + geogridOptical.geogridOptical_Py(self._geogridOptical) self.get_center_latlon() + ##Get parameters + self.getState() + + ##Clean up + self.finalize() + def get_center_latlon(self): ''' Get center lat/lon of the image. ''' from osgeo import gdal + gdal.AllRegister() self.epsgDem = 4326 self.epsgDat = self.getProjectionSystem(self.dat1name) self.determineBbox() @@ -166,604 +178,96 @@ def determineBbox(self, zrange=[-200,4000]): self._ylim = [np.min(xyzs[:,1]), np.max(xyzs[:,1])] - def checkState(self): + def getState(self): + + from . import geogridOptical + + self.pOff = geogridOptical.getXOff_Py(self._geogridOptical) + self.lOff = geogridOptical.getYOff_Py(self._geogridOptical) + self.pCount = geogridOptical.getXCount_Py(self._geogridOptical) + self.lCount = geogridOptical.getYCount_Py(self._geogridOptical) + self.X_res = geogridOptical.getXPixelSize_Py(self._geogridOptical) + self.Y_res = geogridOptical.getYPixelSize_Py(self._geogridOptical) + + def setState(self): ''' Create C object and populate. ''' - if self.repeatTime < 0: - raise Exception('Input image 1 must be older than input image 2') - - - - def geogrid(self): - - # For now print inputs that were obtained - - print("\nOptical Image parameters: ") - print("X-direction coordinate: " + str(self.startingX) + " " + str(self.XSize)) - print("Y-direction coordinate: " + str(self.startingY) + " " + str(self.YSize)) - print("Dimensions: " + str(self.numberOfSamples) + " " + str(self.numberOfLines) + "\n") - - print("Map inputs: ") - print("EPSG: " + str(self.epsgDem)) - print("Smallest Allowable Chip Size in m: " + str(self.chipSizeX0)) - print("Grid spacing in m: " + str(self.gridSpacingX)) - print("Repeat Time: " + str(self.repeatTime)) - print("XLimits: " + str(self._xlim[0]) + " " + str(self._xlim[1])) - print("YLimits: " + str(self._ylim[0]) + " " + str(self._ylim[1])) - print("Extent in km: " + str((self._xlim[1]-self._xlim[0])/1000.0) + " " + str((self._ylim[1]-self._ylim[0])/1000.0)) - if (self.demname != ""): - print("DEM: " + str(self.demname)) - if (self.dhdxname != ""): - print("Slopes: " + str(self.dhdxname) + " " + str(self.dhdyname)) - if (self.vxname != ""): - print("Velocities: " + str(self.vxname) + " " + str(self.vyname)) - if (self.srxname != ""): - print("Search Range: " + str(self.srxname) + " " + str(self.sryname)) - if (self.csminxname != ""): - print("Chip Size Min: " + str(self.csminxname) + " " + str(self.csminyname)) - if (self.csmaxxname != ""): - print("Chip Size Max: " + str(self.csmaxxname) + " " + str(self.csmaxyname)) - if (self.ssmname != ""): - print("Stable Surface Mask: " + str(self.ssmname)) - - - print("\nOutputs: ") - - print("Window locations: " + str(self.winlocname)) - - if (self.dhdxname != ""): - if (self.vxname != ""): - print("Window offsets: " + str(self.winoffname)) - - print("Window rdr_off2vel_x vector: " + str(self.winro2vxname)) - print("Window rdr_off2vel_y vector: " + str(self.winro2vyname)) - - if (self.srxname != ""): - print("Window search range: " + str(self.winsrname)) - - if (self.csminxname != ""): - print("Window chip size min: " + str(self.wincsminname)) - if (self.csmaxxname != ""): - print("Window chip size max: " + str(self.wincsmaxname)) - if (self.ssmname != ""): - print("Window stable surface mask: " + str(self.winssmname)) - - print("Output Nodata Value: " + str(self.nodata_out) + "\n") - - - - print("Starting processing .... ") - - dt_unity = self.srs_dt_unity - max_factor = self.srs_max_scale - upper_thld = self.srs_max_search - lower_thld = self.srs_min_search - - - from osgeo import gdal, osr - import numpy as np - import struct - -# pdb.set_trace() - - demDS = gdal.Open(self.demname, gdal.GA_ReadOnly) - - if (self.dhdxname != ""): - sxDS = gdal.Open(self.dhdxname, gdal.GA_ReadOnly) - syDS = gdal.Open(self.dhdyname, gdal.GA_ReadOnly) - - if (self.vxname != ""): - vxDS = gdal.Open(self.vxname, gdal.GA_ReadOnly) - vyDS = gdal.Open(self.vyname, gdal.GA_ReadOnly) - - if (self.srxname != ""): - srxDS = gdal.Open(self.srxname, gdal.GA_ReadOnly) - sryDS = gdal.Open(self.sryname, gdal.GA_ReadOnly) - - if (self.csminxname != ""): - csminxDS = gdal.Open(self.csminxname, gdal.GA_ReadOnly) - csminyDS = gdal.Open(self.csminyname, gdal.GA_ReadOnly) - - if (self.csmaxxname != ""): - csmaxxDS = gdal.Open(self.csmaxxname, gdal.GA_ReadOnly) - csmaxyDS = gdal.Open(self.csmaxyname, gdal.GA_ReadOnly) - - if (self.ssmname != ""): - ssmDS = gdal.Open(self.ssmname, gdal.GA_ReadOnly) - - if demDS is None: - raise Exception('Error opening DEM file {0}'.format(self.demname)) - - if (self.dhdxname != ""): - if (sxDS is None): - raise Exception('Error opening x-direction slope file {0}'.format(self.dhdxname)) - if (syDS is None): - raise Exception('Error opening y-direction slope file {0}'.format(self.dhdyname)) - - if (self.vxname != ""): - if (vxDS is None): - raise Exception('Error opening x-direction velocity file {0}'.format(self.vxname)) - if (vyDS is None): - raise Exception('Error opening y-direction velocity file {0}'.format(self.vyname)) - - if (self.srxname != ""): - if (srxDS is None): - raise Exception('Error opening x-direction search range file {0}'.format(self.srxname)) - if (sryDS is None): - raise Exception('Error opening y-direction search range file {0}'.format(self.sryname)) - - if (self.csminxname != ""): - if (csminxDS is None): - raise Exception('Error opening x-direction chip size min file {0}'.format(self.csminxname)) - if (csminyDS is None): - raise Exception('Error opening y-direction chip size min file {0}'.format(self.csminyname)) - - if (self.csmaxxname != ""): - if (csmaxxDS is None): - raise Exception('Error opening x-direction chip size max file {0}'.format(self.csmaxxname)) - if (csmaxyDS is None): - raise Exception('Error opening y-direction chip size max file {0}'.format(self.csmaxyname)) - - if (self.ssmname != ""): - if (ssmDS is None): - raise Exception('Error opening stable surface mask file {0}'.format(self.ssmname)) - - geoTrans = demDS.GetGeoTransform() - demXSize = demDS.RasterXSize - demYSize = demDS.RasterYSize - - - # Get offsets and size to read from DEM - lOff = int(np.max( [np.floor((self._ylim[1] - geoTrans[3])/geoTrans[5]), 0.])) -# pdb.set_trace() - lCount = int(np.min([ np.ceil((self._ylim[0] - geoTrans[3])/geoTrans[5]), demYSize-1.]) - lOff) - - pOff = int(np.max([ np.floor((self._xlim[0] - geoTrans[0])/geoTrans[1]), 0.])) - pCount = int(np.min([ np.ceil((self._xlim[1] - geoTrans[0])/geoTrans[1]), demXSize-1.]) - pOff) - - print("Xlimits : " + str(geoTrans[0] + pOff * geoTrans[1]) + " " + str(geoTrans[0] + (pOff + pCount) * geoTrans[1])) - - print("Ylimits : " + str(geoTrans[3] + (lOff + lCount) * geoTrans[5]) + " " + str(geoTrans[3] + lOff * geoTrans[5])) - - print("Origin index (in DEM) of geogrid: " + str(pOff) + " " + str(lOff)) - self.pOff = pOff - self.lOff = lOff - - print("Dimensions of geogrid: " + str(pCount) + " x " + str(lCount)) - self.pCount = pCount - self.lCount = lCount - - projDem = osr.SpatialReference() - if self.epsgDem: - projDem.ImportFromEPSG(self.epsgDem) - else: - raise Exception('EPSG code does not exist for DEM') - - projDat = osr.SpatialReference() - if self.epsgDat: - projDat.ImportFromEPSG(self.epsgDat) - else: - raise Exception('EPSG code does not exist for image data') - fwdTrans = osr.CoordinateTransformation(projDem, projDat) - invTrans = osr.CoordinateTransformation(projDat, projDem) + from . import geogridOptical - if (self.vxname != ""): - nodata = vxDS.GetRasterBand(1).GetNoDataValue() - else: - nodata = 0 - - nodata_out = self.nodata_out - - - pszFormat = "GTiff" - adfGeoTransform = ( geoTrans[0] + pOff * geoTrans[1], geoTrans[1], 0, geoTrans[3] + lOff * geoTrans[5], 0, geoTrans[5] ) - oSRS = osr.SpatialReference() - pszSRS_WKT = projDem.ExportToWkt() - - - - poDriver = gdal.GetDriverByName(pszFormat) - if( poDriver is None ): - raise Exception('Cannot create gdal driver for output') - - pszDstFilename = self.winlocname - poDstDS = poDriver.Create(pszDstFilename, xsize=pCount, ysize=lCount, bands=2, eType=gdal.GDT_Int32) - poDstDS.SetGeoTransform( adfGeoTransform ) - poDstDS.SetProjection( pszSRS_WKT ) - - poBand1 = poDstDS.GetRasterBand(1) - poBand2 = poDstDS.GetRasterBand(2) - poBand1.SetNoDataValue(nodata_out) - poBand2.SetNoDataValue(nodata_out) - - - - if ((self.dhdxname != "")&(self.vxname != "")): - poDriverOff = gdal.GetDriverByName(pszFormat) - if( poDriverOff is None ): - raise Exception('Cannot create gdal driver for output') - - pszDstFilenameOff = self.winoffname - poDstDSOff = poDriverOff.Create(pszDstFilenameOff, xsize=pCount, ysize=lCount, bands=2, eType=gdal.GDT_Int32) - poDstDSOff.SetGeoTransform( adfGeoTransform ) - poDstDSOff.SetProjection( pszSRS_WKT ) - - poBand1Off = poDstDSOff.GetRasterBand(1) - poBand2Off = poDstDSOff.GetRasterBand(2) - poBand1Off.SetNoDataValue(nodata_out) - poBand2Off.SetNoDataValue(nodata_out) - - - if ((self.dhdxname != "")&(self.srxname != "")): - poDriverSch = gdal.GetDriverByName(pszFormat) - if( poDriverSch is None ): - raise Exception('Cannot create gdal driver for output') - - pszDstFilenameSch = self.winsrname - poDstDSSch = poDriverSch.Create(pszDstFilenameSch, xsize=pCount, ysize=lCount, bands=2, eType=gdal.GDT_Int32) - poDstDSSch.SetGeoTransform( adfGeoTransform ) - poDstDSSch.SetProjection( pszSRS_WKT ) - - poBand1Sch = poDstDSSch.GetRasterBand(1) - poBand2Sch = poDstDSSch.GetRasterBand(2) - poBand1Sch.SetNoDataValue(nodata_out) - poBand2Sch.SetNoDataValue(nodata_out) - - if (self.csminxname != ""): - poDriverMin = gdal.GetDriverByName(pszFormat) - if( poDriverMin is None ): - raise Exception('Cannot create gdal driver for output') - - pszDstFilenameMin = self.wincsminname - poDstDSMin = poDriverMin.Create(pszDstFilenameMin, xsize=pCount, ysize=lCount, bands=2, eType=gdal.GDT_Int32) - poDstDSMin.SetGeoTransform( adfGeoTransform ) - poDstDSMin.SetProjection( pszSRS_WKT ) - - poBand1Min = poDstDSMin.GetRasterBand(1) - poBand2Min = poDstDSMin.GetRasterBand(2) - poBand1Min.SetNoDataValue(nodata_out) - poBand2Min.SetNoDataValue(nodata_out) - - if (self.csmaxxname != ""): - poDriverMax = gdal.GetDriverByName(pszFormat) - if( poDriverMax is None ): - raise Exception('Cannot create gdal driver for output') - - pszDstFilenameMax = self.wincsmaxname - poDstDSMax = poDriverMax.Create(pszDstFilenameMax, xsize=pCount, ysize=lCount, bands=2, eType=gdal.GDT_Int32) - poDstDSMax.SetGeoTransform( adfGeoTransform ) - poDstDSMax.SetProjection( pszSRS_WKT ) - - poBand1Max = poDstDSMax.GetRasterBand(1) - poBand2Max = poDstDSMax.GetRasterBand(2) - poBand1Max.SetNoDataValue(nodata_out) - poBand2Max.SetNoDataValue(nodata_out) - - - if (self.ssmname != ""): - poDriverMsk = gdal.GetDriverByName(pszFormat) - if( poDriverMsk is None ): - raise Exception('Cannot create gdal driver for output') - - pszDstFilenameMsk = self.winssmname - poDstDSMsk = poDriverMsk.Create(pszDstFilenameMsk, xsize=pCount, ysize=lCount, bands=1, eType=gdal.GDT_Int32) - poDstDSMsk.SetGeoTransform( adfGeoTransform ) - poDstDSMsk.SetProjection( pszSRS_WKT ) - - poBand1Msk = poDstDSMsk.GetRasterBand(1) - poBand1Msk.SetNoDataValue(nodata_out) - - - - - if (self.dhdxname != ""): - poDriverRO2VX = gdal.GetDriverByName(pszFormat) - if( poDriverRO2VX is None ): - raise Exception('Cannot create gdal driver for output') - - pszDstFilenameRO2VX = self.winro2vxname - poDstDSRO2VX = poDriverRO2VX.Create(pszDstFilenameRO2VX, xsize=pCount, ysize=lCount, bands=2, eType=gdal.GDT_Float64) - poDstDSRO2VX.SetGeoTransform( adfGeoTransform ) - poDstDSRO2VX.SetProjection( pszSRS_WKT ) - - poBand1RO2VX = poDstDSRO2VX.GetRasterBand(1) - poBand2RO2VX = poDstDSRO2VX.GetRasterBand(2) - poBand1RO2VX.SetNoDataValue(nodata_out) - poBand2RO2VX.SetNoDataValue(nodata_out) - - - poDriverRO2VY = gdal.GetDriverByName(pszFormat) - if( poDriverRO2VY is None ): - raise Exception('Cannot create gdal driver for output') - - pszDstFilenameRO2VY = self.winro2vyname - poDstDSRO2VY = poDriverRO2VY.Create(pszDstFilenameRO2VY, xsize=pCount, ysize=lCount, bands=2, eType=gdal.GDT_Float64) - poDstDSRO2VY.SetGeoTransform( adfGeoTransform ) - poDstDSRO2VY.SetProjection( pszSRS_WKT ) - - poBand1RO2VY = poDstDSRO2VY.GetRasterBand(1) - poBand2RO2VY = poDstDSRO2VY.GetRasterBand(2) - poBand1RO2VY.SetNoDataValue(nodata_out) - poBand2RO2VY.SetNoDataValue(nodata_out) - - - - raster1 = np.zeros(pCount,dtype=np.int32) - raster2 = np.zeros(pCount,dtype=np.int32) - raster11 = np.zeros(pCount,dtype=np.int32) - raster22 = np.zeros(pCount,dtype=np.int32) - sr_raster11 = np.zeros(pCount,dtype=np.int32) - sr_raster22 = np.zeros(pCount,dtype=np.int32) - csmin_raster11 = np.zeros(pCount,dtype=np.int32) - csmin_raster22 = np.zeros(pCount,dtype=np.int32) - csmax_raster11 = np.zeros(pCount,dtype=np.int32) - csmax_raster22 = np.zeros(pCount,dtype=np.int32) - ssm_raster = np.zeros(pCount,dtype=np.int32) - raster1a = np.zeros(pCount,dtype=np.float64) - raster1b = np.zeros(pCount,dtype=np.float64) - raster2a = np.zeros(pCount,dtype=np.float64) - raster2b = np.zeros(pCount,dtype=np.float64) - - - - # X- and Y-direction pixel size - X_res = np.abs(self.XSize) - Y_res = np.abs(self.YSize) - print("X-direction pixel size: " + str(X_res)) - print("Y-direction pixel size: " + str(Y_res)) - self.X_res = X_res - self.Y_res = Y_res - - ChipSizeX0_PIX_X = np.ceil(self.chipSizeX0 / X_res / 4) * 4 - ChipSizeX0_PIX_Y = np.ceil(self.chipSizeX0 / Y_res / 4) * 4 - - - - - - for ii in range(lCount): - y = geoTrans[3] + (lOff+ii+0.5) * geoTrans[5] - demLine = demDS.GetRasterBand(1).ReadRaster(xoff=pOff, yoff=lOff+ii, xsize=pCount, ysize=1, buf_xsize=pCount, buf_ysize=1, buf_type=gdal.GDT_Float64) - demLine = struct.unpack('d' * pCount, demLine) - - if (self.dhdxname != ""): - sxLine = sxDS.GetRasterBand(1).ReadRaster(xoff=pOff, yoff=lOff+ii, xsize=pCount, ysize=1, buf_xsize=pCount, buf_ysize=1, buf_type=gdal.GDT_Float64) - sxLine = struct.unpack('d' * pCount, sxLine) - syLine = syDS.GetRasterBand(1).ReadRaster(xoff=pOff, yoff=lOff+ii, xsize=pCount, ysize=1, buf_xsize=pCount, buf_ysize=1, buf_type=gdal.GDT_Float64) - syLine = struct.unpack('d' * pCount, syLine) - - if (self.vxname != ""): - vxLine = vxDS.GetRasterBand(1).ReadRaster(xoff=pOff, yoff=lOff+ii, xsize=pCount, ysize=1, buf_xsize=pCount, buf_ysize=1, buf_type=gdal.GDT_Float64) - vxLine = struct.unpack('d' * pCount, vxLine) - vyLine = vyDS.GetRasterBand(1).ReadRaster(xoff=pOff, yoff=lOff+ii, xsize=pCount, ysize=1, buf_xsize=pCount, buf_ysize=1, buf_type=gdal.GDT_Float64) - vyLine = struct.unpack('d' * pCount, vyLine) - - if (self.srxname != ""): - srxLine = srxDS.GetRasterBand(1).ReadRaster(xoff=pOff, yoff=lOff+ii, xsize=pCount, ysize=1, buf_xsize=pCount, buf_ysize=1, buf_type=gdal.GDT_Float64) - srxLine = struct.unpack('d' * pCount, srxLine) - sryLine = sryDS.GetRasterBand(1).ReadRaster(xoff=pOff, yoff=lOff+ii, xsize=pCount, ysize=1, buf_xsize=pCount, buf_ysize=1, buf_type=gdal.GDT_Float64) - sryLine = struct.unpack('d' * pCount, sryLine) - - if (self.csminxname != ""): - csminxLine = csminxDS.GetRasterBand(1).ReadRaster(xoff=pOff, yoff=lOff+ii, xsize=pCount, ysize=1, buf_xsize=pCount, buf_ysize=1, buf_type=gdal.GDT_Float64) - csminxLine = struct.unpack('d' * pCount, csminxLine) - csminyLine = csminyDS.GetRasterBand(1).ReadRaster(xoff=pOff, yoff=lOff+ii, xsize=pCount, ysize=1, buf_xsize=pCount, buf_ysize=1, buf_type=gdal.GDT_Float64) - csminyLine = struct.unpack('d' * pCount, csminyLine) - - if (self.csmaxxname != ""): - csmaxxLine = csmaxxDS.GetRasterBand(1).ReadRaster(xoff=pOff, yoff=lOff+ii, xsize=pCount, ysize=1, buf_xsize=pCount, buf_ysize=1, buf_type=gdal.GDT_Float64) - csmaxxLine = struct.unpack('d' * pCount, csmaxxLine) - csmaxyLine = csmaxyDS.GetRasterBand(1).ReadRaster(xoff=pOff, yoff=lOff+ii, xsize=pCount, ysize=1, buf_xsize=pCount, buf_ysize=1, buf_type=gdal.GDT_Float64) - csmaxyLine = struct.unpack('d' * pCount, csmaxyLine) - - if (self.ssmname != ""): - ssmLine = ssmDS.GetRasterBand(1).ReadRaster(xoff=pOff, yoff=lOff+ii, xsize=pCount, ysize=1, buf_xsize=pCount, buf_ysize=1, buf_type=gdal.GDT_Float64) - ssmLine = struct.unpack('d' * pCount, ssmLine) - - for jj in range(pCount): - xyzs = np.array([geoTrans[0] + (jj+pOff+0.5)*geoTrans[1], y, demLine[jj]]) - targxyz0 = xyzs.copy() - if (self.dhdxname != ""): - slp = np.array([sxLine[jj], syLine[jj], -1.0]) - if (self.vxname != ""): - vel = np.array([vxLine[jj], vyLine[jj], 0.0]) - else: - vel = np.array([0., 0., 0.]) - if (self.srxname != ""): - schrng1 = np.array([srxLine[jj], sryLine[jj], 0.0]) - schrng1[0] *= np.max((max_factor*((dt_unity-1)*max_factor+(max_factor-1)-(max_factor-1)*self.repeatTime/24.0/3600.0)/((dt_unity-1)*max_factor),1)) - schrng1[0] = np.min((np.max((schrng1[0],lower_thld)),upper_thld)) - schrng1[1] *= np.max((max_factor*((dt_unity-1)*max_factor+(max_factor-1)-(max_factor-1)*self.repeatTime/24.0/3600.0)/((dt_unity-1)*max_factor),1)) - schrng1[1] = np.min((np.max((schrng1[1],lower_thld)),upper_thld)) - schrng2 = np.array([-schrng1[0], schrng1[1], 0.0]) - targutm0 = np.array(fwdTrans.TransformPoint(targxyz0[0],targxyz0[1],targxyz0[2])) - xind = np.round((targutm0[0] - self.startingX) / self.XSize) + 1. - yind = np.round((targutm0[1] - self.startingY) / self.YSize) + 1. - - # x-direction vector - targutm = targutm0.copy() - targutm[0] = targutm0[0] + self.XSize - targxyz = np.array(invTrans.TransformPoint(targutm[0],targutm[1],targutm[2])) - xunit = (targxyz-targxyz0) / np.linalg.norm(targxyz-targxyz0) - - # y-direction vector - targutm = targutm0.copy() - targutm[1] = targutm0[1] + self.YSize - targxyz = np.array(invTrans.TransformPoint(targutm[0],targutm[1],targutm[2])) - yunit = (targxyz-targxyz0) / np.linalg.norm(targxyz-targxyz0) - - # local normal vector - if (self.dhdxname != ""): - normal = -slp / np.linalg.norm(slp) - else: - normal = np.array([0., 0., 0.]) - - if (self.vxname != ""): - vel[2] = -(vel[0]*normal[0]+vel[1]*normal[1])/normal[2] - - if (self.srxname != ""): - schrng1[2] = -(schrng1[0]*normal[0]+schrng1[1]*normal[1])/normal[2] - schrng2[2] = -(schrng2[0]*normal[0]+schrng2[1]*normal[1])/normal[2] - - - - if ((xind > self.numberOfSamples)|(xind < 1)|(yind > self.numberOfLines)|(yind < 1)): -# pdb.set_trace() - raster1[jj] = nodata_out - raster2[jj] = nodata_out - raster11[jj] = nodata_out - raster22[jj] = nodata_out - - sr_raster11[jj] = nodata_out - sr_raster22[jj] = nodata_out - csmin_raster11[jj] = nodata_out - csmin_raster22[jj] = nodata_out - csmax_raster11[jj] = nodata_out - csmax_raster22[jj] = nodata_out - ssm_raster[jj] = nodata_out - - raster1a[jj] = nodata_out - raster1b[jj] = nodata_out - raster2a[jj] = nodata_out - raster2b[jj] = nodata_out - else: - raster1[jj] = xind; - raster2[jj] = yind; -# pdb.set_trace() -# if ((self.vxname != "")&(vel[0] != nodata)): -## pdb.set_trace() -# raster11[jj] = np.round(np.dot(vel,xunit)*self.repeatTime/self.XSize/365.0/24.0/3600.0*1) -# raster22[jj] = np.round(np.dot(vel,yunit)*self.repeatTime/self.YSize/365.0/24.0/3600.0*1) -# else: -# raster11[jj] = 0. -# raster22[jj] = 0. - if (self.dhdxname != ""): - - if (self.vxname != ""): - if (vel[0] == nodata): - raster11[jj] = 0. - raster22[jj] = 0. - else: - raster11[jj] = np.round(np.dot(vel,xunit)*self.repeatTime/self.XSize/365.0/24.0/3600.0*1) - raster22[jj] = np.round(np.dot(vel,yunit)*self.repeatTime/self.YSize/365.0/24.0/3600.0*1) - - cross = np.cross(xunit,yunit) - cross = cross / np.linalg.norm(cross) - cross_check = np.abs(np.arccos(np.dot(normal,cross))/np.pi*180.0-90.0) - - if (cross_check > 1.0): - raster1a[jj] = normal[2]/(self.repeatTime/self.XSize/365.0/24.0/3600.0)*(normal[2]*yunit[1]-normal[1]*yunit[2])/((normal[2]*xunit[0]-normal[0]*xunit[2])*(normal[2]*yunit[1]-normal[1]*yunit[2])-(normal[2]*yunit[0]-normal[0]*yunit[2])*(normal[2]*xunit[1]-normal[1]*xunit[2])); - raster1b[jj] = -normal[2]/(self.repeatTime/self.YSize/365.0/24.0/3600.0)*(normal[2]*xunit[1]-normal[1]*xunit[2])/((normal[2]*xunit[0]-normal[0]*xunit[2])*(normal[2]*yunit[1]-normal[1]*yunit[2])-(normal[2]*yunit[0]-normal[0]*yunit[2])*(normal[2]*xunit[1]-normal[1]*xunit[2])); - raster2a[jj] = -normal[2]/(self.repeatTime/self.XSize/365.0/24.0/3600.0)*(normal[2]*yunit[0]-normal[0]*yunit[2])/((normal[2]*xunit[0]-normal[0]*xunit[2])*(normal[2]*yunit[1]-normal[1]*yunit[2])-(normal[2]*yunit[0]-normal[0]*yunit[2])*(normal[2]*xunit[1]-normal[1]*xunit[2])); - raster2b[jj] = normal[2]/(self.repeatTime/self.YSize/365.0/24.0/3600.0)*(normal[2]*xunit[0]-normal[0]*xunit[2])/((normal[2]*xunit[0]-normal[0]*xunit[2])*(normal[2]*yunit[1]-normal[1]*yunit[2])-(normal[2]*yunit[0]-normal[0]*yunit[2])*(normal[2]*xunit[1]-normal[1]*xunit[2])); - else: - raster1a[jj] = nodata_out - raster1b[jj] = nodata_out - raster2a[jj] = nodata_out - raster2b[jj] = nodata_out - - if (self.srxname != ""): - if ((self.vxname != "")&(vel[0] == nodata)): - sr_raster11[jj] = 0 - sr_raster22[jj] = 0 - else: - sr_raster11[jj] = np.abs(np.round(np.dot(schrng1,xunit)*self.repeatTime/self.XSize/365.0/24.0/3600.0*1)) - sr_raster22[jj] = np.abs(np.round(np.dot(schrng1,yunit)*self.repeatTime/self.YSize/365.0/24.0/3600.0*1)) - if (np.abs(np.round(np.dot(schrng2,xunit)*self.repeatTime/self.XSize/365.0/24.0/3600.0*1)) > sr_raster11[jj]): - sr_raster11[jj] = np.abs(np.round(np.dot(schrng2,xunit)*self.repeatTime/self.XSize/365.0/24.0/3600.0*1)) - if (np.abs(np.round(np.dot(schrng2,yunit)*self.repeatTime/self.YSize/365.0/24.0/3600.0*1)) > sr_raster22[jj]): - sr_raster22[jj] = np.abs(np.round(np.dot(schrng2,yunit)*self.repeatTime/self.YSize/365.0/24.0/3600.0*1)) - if (sr_raster11[jj] == 0): - sr_raster11[jj] = 1 - if (sr_raster22[jj] == 0): - sr_raster22[jj] = 1 - - if (self.csminxname != ""): - csmin_raster11[jj] = csminxLine[jj] / self.chipSizeX0 * ChipSizeX0_PIX_X - csmin_raster22[jj] = csminyLine[jj] / self.chipSizeX0 * ChipSizeX0_PIX_Y - - - if (self.csmaxxname != ""): - csmax_raster11[jj] = csmaxxLine[jj] / self.chipSizeX0 * ChipSizeX0_PIX_X - csmax_raster22[jj] = csmaxyLine[jj] / self.chipSizeX0 * ChipSizeX0_PIX_Y - - - - if (self.ssmname != ""): - ssm_raster[jj] = ssmLine[jj] - - - - - - -# pdb.set_trace() - - poBand1.WriteRaster(xoff=0, yoff=ii, xsize=pCount, ysize=1, buf_len=raster1.tostring(), buf_xsize=pCount, buf_ysize=1, buf_type=gdal.GDT_Int32) - poBand2.WriteRaster(xoff=0, yoff=ii, xsize=pCount, ysize=1, buf_len=raster2.tostring(), buf_xsize=pCount, buf_ysize=1, buf_type=gdal.GDT_Int32) - if ((self.dhdxname != "")&(self.vxname != "")): - poBand1Off.WriteRaster(xoff=0, yoff=ii, xsize=pCount, ysize=1, buf_len=raster11.tostring(), buf_xsize=pCount, buf_ysize=1, buf_type=gdal.GDT_Int32) - poBand2Off.WriteRaster(xoff=0, yoff=ii, xsize=pCount, ysize=1, buf_len=raster22.tostring(), buf_xsize=pCount, buf_ysize=1, buf_type=gdal.GDT_Int32) - if ((self.dhdxname != "")&(self.srxname != "")): - poBand1Sch.WriteRaster(xoff=0, yoff=ii, xsize=pCount, ysize=1, buf_len=sr_raster11.tostring(), buf_xsize=pCount, buf_ysize=1, buf_type=gdal.GDT_Int32) - poBand2Sch.WriteRaster(xoff=0, yoff=ii, xsize=pCount, ysize=1, buf_len=sr_raster22.tostring(), buf_xsize=pCount, buf_ysize=1, buf_type=gdal.GDT_Int32) - if (self.csminxname != ""): - poBand1Min.WriteRaster(xoff=0, yoff=ii, xsize=pCount, ysize=1, buf_len=csmin_raster11.tostring(), buf_xsize=pCount, buf_ysize=1, buf_type=gdal.GDT_Int32) - poBand2Min.WriteRaster(xoff=0, yoff=ii, xsize=pCount, ysize=1, buf_len=csmin_raster22.tostring(), buf_xsize=pCount, buf_ysize=1, buf_type=gdal.GDT_Int32) - if (self.csmaxxname != ""): - poBand1Max.WriteRaster(xoff=0, yoff=ii, xsize=pCount, ysize=1, buf_len=csmax_raster11.tostring(), buf_xsize=pCount, buf_ysize=1, buf_type=gdal.GDT_Int32) - poBand2Max.WriteRaster(xoff=0, yoff=ii, xsize=pCount, ysize=1, buf_len=csmax_raster22.tostring(), buf_xsize=pCount, buf_ysize=1, buf_type=gdal.GDT_Int32) - if (self.ssmname != ""): - poBand1Msk.WriteRaster(xoff=0, yoff=ii, xsize=pCount, ysize=1, buf_len=ssm_raster.tostring(), buf_xsize=pCount, buf_ysize=1, buf_type=gdal.GDT_Int32) - if (self.dhdxname != ""): - poBand1RO2VX.WriteRaster(xoff=0, yoff=ii, xsize=pCount, ysize=1, buf_len=raster1a.tostring(), buf_xsize=pCount, buf_ysize=1, buf_type=gdal.GDT_Float64) - poBand2RO2VX.WriteRaster(xoff=0, yoff=ii, xsize=pCount, ysize=1, buf_len=raster1b.tostring(), buf_xsize=pCount, buf_ysize=1, buf_type=gdal.GDT_Float64) - poBand1RO2VY.WriteRaster(xoff=0, yoff=ii, xsize=pCount, ysize=1, buf_len=raster2a.tostring(), buf_xsize=pCount, buf_ysize=1, buf_type=gdal.GDT_Float64) - poBand2RO2VY.WriteRaster(xoff=0, yoff=ii, xsize=pCount, ysize=1, buf_len=raster2b.tostring(), buf_xsize=pCount, buf_ysize=1, buf_type=gdal.GDT_Float64) - - - poDstDS = None - if ((self.dhdxname != "")&(self.vxname != "")): - poDstDSOff = None - if ((self.dhdxname != "")&(self.srxname != "")): - poDstDSSch = None - if (self.csminxname != ""): - poDstDSMin = None - if (self.csmaxxname != ""): - poDstDSMax = None - if (self.ssmname != ""): - poDstDSMsk = None - if (self.dhdxname != ""): - poDstDSRO2VX = None + if self._geogridOptical is not None: + geogridOptical.destroyGeoGridOptical_Py(self._geogridOptical) + + self._geogridOptical = geogridOptical.createGeoGridOptical_Py() + geogridOptical.setOpticalImageDimensions_Py( self._geogridOptical, self.numberOfSamples, self.numberOfLines) + geogridOptical.setXParameters_Py( self._geogridOptical, self.startingX, self.XSize) + geogridOptical.setYParameters_Py( self._geogridOptical, self.startingY, self.YSize) + geogridOptical.setRepeatTime_Py(self._geogridOptical, self.repeatTime) + + geogridOptical.setDtUnity_Py( self._geogridOptical, self.srs_dt_unity) + geogridOptical.setMaxFactor_Py( self._geogridOptical, self.srs_max_scale) + geogridOptical.setUpperThreshold_Py( self._geogridOptical, self.srs_max_search) + geogridOptical.setLowerThreshold_Py(self._geogridOptical, self.srs_min_search) + + geogridOptical.setEPSG_Py(self._geogridOptical, self.epsgDem, self.epsgDat) + geogridOptical.setChipSizeX0_Py(self._geogridOptical, self.chipSizeX0) + geogridOptical.setGridSpacingX_Py(self._geogridOptical, self.gridSpacingX) + + geogridOptical.setXLimits_Py(self._geogridOptical, self._xlim[0], self._xlim[1]) + geogridOptical.setYLimits_Py(self._geogridOptical, self._ylim[0], self._ylim[1]) + if self.demname: + geogridOptical.setDEM_Py(self._geogridOptical, self.demname) - poDstDSRO2VY = None + if (self.dhdxname is not None) and (self.dhdyname is not None): + geogridOptical.setSlopes_Py(self._geogridOptical, self.dhdxname, self.dhdyname) + + if (self.vxname is not None) and (self.vyname is not None): + geogridOptical.setVelocities_Py(self._geogridOptical, self.vxname, self.vyname) + + if (self.srxname is not None) and (self.sryname is not None): + geogridOptical.setSearchRange_Py(self._geogridOptical, self.srxname, self.sryname) - demDS = None + if (self.csminxname is not None) and (self.csminyname is not None): + geogridOptical.setChipSizeMin_Py(self._geogridOptical, self.csminxname, self.csminyname) + + if (self.csmaxxname is not None) and (self.csmaxyname is not None): + geogridOptical.setChipSizeMax_Py(self._geogridOptical, self.csmaxxname, self.csmaxyname) - if (self.dhdxname != ""): - sxDS = None - syDS = None + if (self.ssmname is not None): + geogridOptical.setStableSurfaceMask_Py(self._geogridOptical, self.ssmname) + + geogridOptical.setWindowLocationsFilename_Py( self._geogridOptical, self.winlocname) + geogridOptical.setWindowOffsetsFilename_Py( self._geogridOptical, self.winoffname) + geogridOptical.setWindowSearchRangeFilename_Py( self._geogridOptical, self.winsrname) + geogridOptical.setWindowChipSizeMinFilename_Py( self._geogridOptical, self.wincsminname) + geogridOptical.setWindowChipSizeMaxFilename_Py( self._geogridOptical, self.wincsmaxname) + geogridOptical.setWindowStableSurfaceMaskFilename_Py( self._geogridOptical, self.winssmname) + geogridOptical.setRO2VXFilename_Py( self._geogridOptical, self.winro2vxname) + geogridOptical.setRO2VYFilename_Py( self._geogridOptical, self.winro2vyname) + geogridOptical.setNodataOut_Py(self._geogridOptical, self.nodata_out) + + + def checkState(self): + ''' + Create C object and populate. + ''' + if self.repeatTime < 0: + raise Exception('Input image 1 must be older than input image 2') - if (self.vxname != ""): - vxDS = None - vyDS = None - if (self.srxname != ""): - srxDS = None - sryDS = None + def finalize(self): + ''' + Clean up all the C pointers. + ''' - if (self.csminxname != ""): - csminxDS = None - csminyDS = None + from . import geogridOptical + + geogridOptical.destroyGeoGridOptical_Py(self._geogridOptical) + self._geogridOptical = None - if (self.csmaxxname != ""): - csmaxxDS = None - csmaxyDS = None - if (self.ssmname != ""): - ssmDS = None + @@ -882,6 +386,9 @@ def __init__(self): self._xlim = None self._ylim = None self.nodata_out = None + + ##Pointer to C + self._geogridOptical = None ##parameters for autoRIFT self.pOff = None diff --git a/geo_autoRIFT/geogrid/SConscript b/geo_autoRIFT/geogrid/SConscript index c7128ad..a11c390 100755 --- a/geo_autoRIFT/geogrid/SConscript +++ b/geo_autoRIFT/geogrid/SConscript @@ -26,13 +26,28 @@ if not os.path.exists(initFile): fout.write("#!/usr/bin/env python3") fout.close() + listFiles = ['Geogrid.py','GeogridOptical.py',initFile] envgeogrid.Install(install,listFiles) envgeogrid.Alias('install',install) Export('envgeogrid') + + + + bindingsScons = 'bindings/SConscript' SConscript(bindingsScons,variant_dir = envgeogrid['PRJ_SCONS_BUILD'] + '/' + package + '/' + project + '/bindings') includeScons = 'include/SConscript' SConscript(includeScons) srcScons = 'src/SConscript' SConscript(srcScons,variant_dir = envgeogrid['PRJ_SCONS_BUILD'] + '/' + package + '/' + project + '/src') + + + +bindingsScons1 = 'bindings/SConscript1' +SConscript(bindingsScons1,variant_dir = envgeogrid['PRJ_SCONS_BUILD'] + '/' + package + '/' + project + '/bindings') +includeScons1 = 'include/SConscript1' +SConscript(includeScons1) +srcScons1 = 'src/SConscript1' +SConscript(srcScons1,variant_dir = envgeogrid['PRJ_SCONS_BUILD'] + '/' + package + '/' + project + '/src') + diff --git a/geo_autoRIFT/geogrid/__init__.py b/geo_autoRIFT/geogrid/__init__.py index 58dec40..b2cff13 100755 --- a/geo_autoRIFT/geogrid/__init__.py +++ b/geo_autoRIFT/geogrid/__init__.py @@ -4,7 +4,7 @@ # from .splitSpectrum import PySplitRangeSpectrum # return PySplitRangeSpectrum() - +# should always work - standalone or with ISCE from .GeogridOptical import GeogridOptical try: diff --git a/geo_autoRIFT/geogrid/bindings/SConscript1 b/geo_autoRIFT/geogrid/bindings/SConscript1 new file mode 100755 index 0000000..db21694 --- /dev/null +++ b/geo_autoRIFT/geogrid/bindings/SConscript1 @@ -0,0 +1,28 @@ +#!/usr/bin/env python +#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +# +# Author: Piyush Agram +# Copyright 2019, by the California Institute of Technology. ALL RIGHTS RESERVED. United States Government Sponsorship acknowledged. +# Any commercial use must be negotiated with the Office of Technology Transfer at the California Institute of Technology. +# +# This software may be subject to U.S. export control laws. By accepting this software, the user agrees to comply with all applicable U.S. +# export laws and regulations. User has the responsibility to obtain export licenses, or other export authority as may be required before +# exporting such information to foreign countries or providing access to foreign persons. +# +#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + + +import os + +Import('envgeogrid') +package = envgeogrid['PACKAGE'] +project = envgeogrid['PROJECT'] +install = envgeogrid['PRJ_SCONS_INSTALL'] + '/' + package + '/' + project +build = envgeogrid['PRJ_SCONS_BUILD'] + '/' + package + '/' + project +libList = ['gomp','geogridOptical','combinedLib','gdal'] +envgeogrid.PrependUnique(LIBS = libList) +module = envgeogrid.LoadableModule(target = 'geogridOptical.abi3.so', source = 'geogridOpticalmodule.cpp') +envgeogrid.Install(install,module) +envgeogrid.Alias('install',install) +envgeogrid.Install(build,module) +envgeogrid.Alias('build',build) diff --git a/geo_autoRIFT/geogrid/bindings/geogridOpticalmodule.cpp b/geo_autoRIFT/geogrid/bindings/geogridOpticalmodule.cpp new file mode 100755 index 0000000..da87fe7 --- /dev/null +++ b/geo_autoRIFT/geogrid/bindings/geogridOpticalmodule.cpp @@ -0,0 +1,567 @@ +/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + * Copyright 2019 California Institute of Technology. ALL RIGHTS RESERVED. + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + * United States Government Sponsorship acknowledged. This software is subject to + * U.S. export control laws and regulations and has been classified as 'EAR99 NLR' + * (No [Export] License Required except when exporting to an embargoed country, + * end user, or in support of a prohibited end use). By downloading this software, + * the user agrees to comply with all applicable U.S. export laws and regulations. + * The user has the responsibility to obtain export licenses, or other export + * authority as may be required before exporting this software to any 'EAR99' + * embargoed foreign country or citizen of those countries. + * + * Authors: Piyush Agram, Yang Lei + *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + */ + + + +#include +#include +#include "geogridOptical.h" +#include "geogridOpticalmodule.h" + +static const char * const __doc__ = "Python extension for geogrid"; + +PyModuleDef moduledef = { + //header + PyModuleDef_HEAD_INIT, + //name of the module + "geogridOptical", + //module documentation string + __doc__, + //size of the per-interpreter state of the module; + -1, + geogrid_methods, +}; + +//Initialization function for the module +PyMODINIT_FUNC +PyInit_geogridOptical() +{ + PyObject* module = PyModule_Create(&moduledef); + if (!module) + { + return module; + } + return module; +} + +PyObject* createGeoGridOptical(PyObject* self, PyObject *args) +{ + geoGridOptical* ptr = new geoGridOptical; + return Py_BuildValue("K", (uint64_t) ptr); +} + +PyObject* destroyGeoGridOptical(PyObject *self, PyObject *args) +{ + uint64_t ptr; + if (!PyArg_ParseTuple(args, "K", &ptr)) + { + return NULL; + } + + if (((geoGridOptical*)(ptr))!=NULL) + { + delete ((geoGridOptical*)(ptr)); + } + return Py_BuildValue("i", 0); +} + +PyObject* setEPSG(PyObject *self, PyObject *args) +{ + uint64_t ptr; + int code1, code2; + if (!PyArg_ParseTuple(args, "Kii", &ptr, &code1, &code2)) + { + return NULL; + } + ((geoGridOptical*)(ptr))->epsgDem = code1; + ((geoGridOptical*)(ptr))->epsgDat = code2; + return Py_BuildValue("i", 0); +} + + +PyObject* setChipSizeX0(PyObject *self, PyObject *args) +{ + uint64_t ptr; + double chipSizeX0; + if (!PyArg_ParseTuple(args, "Kd", &ptr, &chipSizeX0)) + { + return NULL; + } + ((geoGridOptical*)(ptr))->chipSizeX0 = chipSizeX0; + return Py_BuildValue("i", 0); +} + +PyObject* setGridSpacingX(PyObject *self, PyObject *args) +{ + uint64_t ptr; + double gridSpacingX; + if (!PyArg_ParseTuple(args, "Kd", &ptr, &gridSpacingX)) + { + return NULL; + } + ((geoGridOptical*)(ptr))->gridSpacingX = gridSpacingX; + return Py_BuildValue("i", 0); +} + +PyObject* setRepeatTime(PyObject *self, PyObject *args) +{ + uint64_t ptr; + double repeatTime; + if (!PyArg_ParseTuple(args, "Kd", &ptr, &repeatTime)) + { + return NULL; + } + ((geoGridOptical*)(ptr))->dt = repeatTime; + return Py_BuildValue("i", 0); +} + +PyObject* setOpticalImageDimensions(PyObject *self, PyObject *args) +{ + uint64_t ptr; + int wid, len; + if (!PyArg_ParseTuple(args, "Kii", &ptr, &wid, &len)) + { + return NULL; + } + + ((geoGridOptical*)(ptr))->nPixels = wid; + ((geoGridOptical*)(ptr))->nLines = len; + return Py_BuildValue("i", 0); +} + +PyObject* setXParameters(PyObject *self, PyObject *args) +{ + uint64_t ptr; + double r0, rspace; + if (!PyArg_ParseTuple(args, "Kdd", &ptr, &r0, &rspace)) + { + return NULL; + } + + ((geoGridOptical*)(ptr))->startingX = r0; + ((geoGridOptical*)(ptr))->XSize = rspace; + return Py_BuildValue("i", 0); +} + +PyObject* setYParameters(PyObject *self, PyObject *args) +{ + uint64_t ptr; + double t0, prf; + if (!PyArg_ParseTuple(args, "Kdd", &ptr, &t0, &prf)) + { + return NULL; + } + + ((geoGridOptical*)(ptr))->startingY = t0; + ((geoGridOptical*)(ptr))->YSize = prf; + return Py_BuildValue("i", 0); +} + +PyObject* setXLimits(PyObject *self, PyObject *args) +{ + uint64_t ptr; + double x0, x1; + if (!PyArg_ParseTuple(args, "Kdd", &ptr, &x0, &x1)) + { + return NULL; + } + + ((geoGridOptical*)(ptr))->xmin = x0; + ((geoGridOptical*)(ptr))->xmax = x1; + return Py_BuildValue("i", 0); +} + +PyObject* setYLimits(PyObject *self, PyObject *args) +{ + uint64_t ptr; + double x0, x1; + if (!PyArg_ParseTuple(args, "Kdd", &ptr, &x0, &x1)) + { + return NULL; + } + + ((geoGridOptical*)(ptr))->ymin = x0; + ((geoGridOptical*)(ptr))->ymax = x1; + return Py_BuildValue("i", 0); +} + +PyObject* setDEM(PyObject *self, PyObject *args) +{ + uint64_t ptr; + char* name; + if (!PyArg_ParseTuple(args, "Ks", &ptr, &name)) + { + return NULL; + } + + ((geoGridOptical*)(ptr))->demname = std::string(name); + return Py_BuildValue("i", 0); +} + +PyObject* setVelocities(PyObject *self, PyObject* args) +{ + uint64_t ptr; + char *vx; + char *vy; + if (!PyArg_ParseTuple(args, "Kss", &ptr, &vx, &vy)) + { + return NULL; + } + + ((geoGridOptical*)(ptr))->vxname = std::string(vx); + ((geoGridOptical*)(ptr))->vyname = std::string(vy); + return Py_BuildValue("i", 0); +} + +PyObject* setSearchRange(PyObject *self, PyObject* args) +{ + uint64_t ptr; + char *srx; + char *sry; + if (!PyArg_ParseTuple(args, "Kss", &ptr, &srx, &sry)) + { + return NULL; + } + + ((geoGridOptical*)(ptr))->srxname = std::string(srx); + ((geoGridOptical*)(ptr))->sryname = std::string(sry); + return Py_BuildValue("i", 0); +} + +PyObject* setChipSizeMin(PyObject *self, PyObject* args) +{ + uint64_t ptr; + char *csminx; + char *csminy; + if (!PyArg_ParseTuple(args, "Kss", &ptr, &csminx, &csminy)) + { + return NULL; + } + + ((geoGridOptical*)(ptr))->csminxname = std::string(csminx); + ((geoGridOptical*)(ptr))->csminyname = std::string(csminy); + return Py_BuildValue("i", 0); +} + +PyObject* setChipSizeMax(PyObject *self, PyObject* args) +{ + uint64_t ptr; + char *csmaxx; + char *csmaxy; + if (!PyArg_ParseTuple(args, "Kss", &ptr, &csmaxx, &csmaxy)) + { + return NULL; + } + + ((geoGridOptical*)(ptr))->csmaxxname = std::string(csmaxx); + ((geoGridOptical*)(ptr))->csmaxyname = std::string(csmaxy); + return Py_BuildValue("i", 0); +} + +PyObject* setStableSurfaceMask(PyObject *self, PyObject* args) +{ + uint64_t ptr; + char *ssm; + if (!PyArg_ParseTuple(args, "Ks", &ptr, &ssm)) + { + return NULL; + } + + ((geoGridOptical*)(ptr))->ssmname = std::string(ssm); + return Py_BuildValue("i", 0); +} + +PyObject* setSlopes(PyObject *self, PyObject* args) +{ + uint64_t ptr; + char *sx; + char *sy; + if (!PyArg_ParseTuple(args, "Kss", &ptr, &sx, &sy)) + { + return NULL; + } + + ((geoGridOptical*)(ptr))->dhdxname = std::string(sx); + ((geoGridOptical*)(ptr))->dhdyname = std::string(sy); + return Py_BuildValue("i", 0); +} + +PyObject* setWindowLocationsFilename(PyObject *self, PyObject *args) +{ + uint64_t ptr; + char* name; + if (!PyArg_ParseTuple(args, "Ks", &ptr, &name)) + { + return NULL; + } + + ((geoGridOptical*)(ptr))->pixlinename = std::string(name); + return Py_BuildValue("i", 0); +} + +PyObject* setWindowOffsetsFilename(PyObject *self, PyObject *args) +{ + uint64_t ptr; + char* name; + if (!PyArg_ParseTuple(args, "Ks", &ptr, &name)) + { + return NULL; + } + + ((geoGridOptical*)(ptr))->offsetname = std::string(name); + return Py_BuildValue("i", 0); +} + +PyObject* setWindowSearchRangeFilename(PyObject *self, PyObject *args) +{ + uint64_t ptr; + char* name; + if (!PyArg_ParseTuple(args, "Ks", &ptr, &name)) + { + return NULL; + } + + ((geoGridOptical*)(ptr))->searchrangename = std::string(name); + return Py_BuildValue("i", 0); +} + +PyObject* setWindowChipSizeMinFilename(PyObject *self, PyObject *args) +{ + uint64_t ptr; + char* name; + if (!PyArg_ParseTuple(args, "Ks", &ptr, &name)) + { + return NULL; + } + + ((geoGridOptical*)(ptr))->chipsizeminname = std::string(name); + return Py_BuildValue("i", 0); +} + +PyObject* setWindowChipSizeMaxFilename(PyObject *self, PyObject *args) +{ + uint64_t ptr; + char* name; + if (!PyArg_ParseTuple(args, "Ks", &ptr, &name)) + { + return NULL; + } + + ((geoGridOptical*)(ptr))->chipsizemaxname = std::string(name); + return Py_BuildValue("i", 0); +} + +PyObject* setWindowStableSurfaceMaskFilename(PyObject *self, PyObject *args) +{ + uint64_t ptr; + char* name; + if (!PyArg_ParseTuple(args, "Ks", &ptr, &name)) + { + return NULL; + } + + ((geoGridOptical*)(ptr))->stablesurfacemaskname = std::string(name); + return Py_BuildValue("i", 0); +} + +PyObject* setRO2VXFilename(PyObject *self, PyObject *args) +{ + uint64_t ptr; + char* name; + if (!PyArg_ParseTuple(args, "Ks", &ptr, &name)) + { + return NULL; + } + + ((geoGridOptical*)(ptr))->ro2vx_name = std::string(name); + return Py_BuildValue("i", 0); +} + +PyObject* setRO2VYFilename(PyObject *self, PyObject *args) +{ + uint64_t ptr; + char* name; + if (!PyArg_ParseTuple(args, "Ks", &ptr, &name)) + { + return NULL; + } + + ((geoGridOptical*)(ptr))->ro2vy_name = std::string(name); + return Py_BuildValue("i", 0); +} + + +PyObject* setNodataOut(PyObject *self, PyObject *args) +{ + uint64_t ptr; + int nodata; + if (!PyArg_ParseTuple(args, "Ki", &ptr, &nodata)) + { + return NULL; + } + + ((geoGridOptical*)(ptr))->nodata_out = nodata; + return Py_BuildValue("i", 0); +} + + + +PyObject* getXOff(PyObject *self, PyObject *args) +{ + int var; + uint64_t ptr; + + if (!PyArg_ParseTuple(args, "K", &ptr)) + { + return NULL; + } + + var = ((geoGridOptical*)(ptr))->pOff; + return Py_BuildValue("i",var); +} + +PyObject* getYOff(PyObject *self, PyObject *args) +{ + int var; + uint64_t ptr; + + if (!PyArg_ParseTuple(args, "K", &ptr)) + { + return NULL; + } + + var = ((geoGridOptical*)(ptr))->lOff; + return Py_BuildValue("i",var); +} + +PyObject* getXCount(PyObject *self, PyObject *args) +{ + int var; + uint64_t ptr; + + if (!PyArg_ParseTuple(args, "K", &ptr)) + { + return NULL; + } + + var = ((geoGridOptical*)(ptr))->pCount; + return Py_BuildValue("i",var); +} + +PyObject* getYCount(PyObject *self, PyObject *args) +{ + int var; + uint64_t ptr; + + if (!PyArg_ParseTuple(args, "K", &ptr)) + { + return NULL; + } + + var = ((geoGridOptical*)(ptr))->lCount; + return Py_BuildValue("i",var); +} + +PyObject* getXPixelSize(PyObject *self, PyObject *args) +{ + double var; + uint64_t ptr; + + if (!PyArg_ParseTuple(args, "K", &ptr)) + { + return NULL; + } + + var = ((geoGridOptical*)(ptr))->X_res; + return Py_BuildValue("d",var); +} + +PyObject* getYPixelSize(PyObject *self, PyObject *args) +{ + double var; + uint64_t ptr; + + if (!PyArg_ParseTuple(args, "K", &ptr)) + { + return NULL; + } + + var = ((geoGridOptical*)(ptr))->Y_res; + return Py_BuildValue("d",var); +} + +PyObject* setDtUnity(PyObject *self, PyObject *args) +{ + uint64_t ptr; + double dt_unity; + if (!PyArg_ParseTuple(args, "Kd", &ptr, &dt_unity)) + { + return NULL; + } + ((geoGridOptical*)(ptr))->dt_unity = dt_unity; + return Py_BuildValue("i", 0); +} + +PyObject* setMaxFactor(PyObject *self, PyObject *args) +{ + uint64_t ptr; + double max_factor; + if (!PyArg_ParseTuple(args, "Kd", &ptr, &max_factor)) + { + return NULL; + } + ((geoGridOptical*)(ptr))->max_factor = max_factor; + return Py_BuildValue("i", 0); +} + +PyObject* setUpperThreshold(PyObject *self, PyObject *args) +{ + uint64_t ptr; + double upper_thld; + if (!PyArg_ParseTuple(args, "Kd", &ptr, &upper_thld)) + { + return NULL; + } + ((geoGridOptical*)(ptr))->upper_thld = upper_thld; + return Py_BuildValue("i", 0); +} + +PyObject* setLowerThreshold(PyObject *self, PyObject *args) +{ + uint64_t ptr; + double lower_thld; + if (!PyArg_ParseTuple(args, "Kd", &ptr, &lower_thld)) + { + return NULL; + } + ((geoGridOptical*)(ptr))->lower_thld = lower_thld; + return Py_BuildValue("i", 0); +} + + +PyObject* geogridOptical(PyObject* self, PyObject* args) +{ + uint64_t ptr; + if (!PyArg_ParseTuple(args, "K", &ptr)) + { + return NULL; + } + + ((geoGridOptical*)(ptr))->geogridOptical(); + return Py_BuildValue("i", 0); +} diff --git a/geo_autoRIFT/geogrid/include/SConscript1 b/geo_autoRIFT/geogrid/include/SConscript1 new file mode 100755 index 0000000..afa3d57 --- /dev/null +++ b/geo_autoRIFT/geogrid/include/SConscript1 @@ -0,0 +1,24 @@ +#!/usr/bin/env python +#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +# +# Author: Piyush Agram +# Copyright 2019, by the California Institute of Technology. ALL RIGHTS RESERVED. United States Government Sponsorship acknowledged. +# Any commercial use must be negotiated with the Office of Technology Transfer at the California Institute of Technology. +# +# This software may be subject to U.S. export control laws. By accepting this software, the user agrees to comply with all applicable U.S. +# export laws and regulations. User has the responsibility to obtain export licenses, or other export authority as may be required before +# exporting such information to foreign countries or providing access to foreign persons. +# +#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + + +import os + +Import('envgeogrid') +package = envgeogrid['PACKAGE'] +project = envgeogrid['PROJECT'] +build = envgeogrid['PRJ_SCONS_BUILD'] + '/' + package + '/' + project + '/include' +envgeogrid.AppendUnique(CPPPATH = [build]) +listFiles = ['geogridOptical.h', 'geogridOpticalmodule.h'] +envgeogrid.Install(build,listFiles) +envgeogrid.Alias('build',build) diff --git a/geo_autoRIFT/geogrid/include/geogridOptical.h b/geo_autoRIFT/geogrid/include/geogridOptical.h new file mode 100755 index 0000000..5dc1957 --- /dev/null +++ b/geo_autoRIFT/geogrid/include/geogridOptical.h @@ -0,0 +1,93 @@ +/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + * Copyright 2019 California Institute of Technology. ALL RIGHTS RESERVED. + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + * United States Government Sponsorship acknowledged. This software is subject to + * U.S. export control laws and regulations and has been classified as 'EAR99 NLR' + * (No [Export] License Required except when exporting to an embargoed country, + * end user, or in support of a prohibited end use). By downloading this software, + * the user agrees to comply with all applicable U.S. export laws and regulations. + * The user has the responsibility to obtain export licenses, or other export + * authority as may be required before exporting this software to any 'EAR99' + * embargoed foreign country or citizen of those countries. + * + * Authors: Piyush Agram, Yang Lei + *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + */ + +#ifndef GEOGRIDOPTICAL_H +#define GEOGRIDOPTICAL_H + +#include + + + +struct geoGridOptical +{ + //DEM related inputs + std::string demname; //DEM + std::string dhdxname; //Slope in X + std::string dhdyname; //Slope in Y + std::string vxname; //Velocity in X + std::string vyname; //Velocity in Y + std::string srxname; //Search range in X + std::string sryname; //Search range in Y + std::string csminxname; //Chip size minimum in x + std::string csminyname; //Chip size minimum in y + std::string csmaxxname; //Chip size maximum in x + std::string csmaxyname; //Chip size maximum in y + std::string ssmname; //Stable surface mask + int epsgDem, epsgDat; + double chipSizeX0; + double gridSpacingX; + + //Bounding box related + double xmin, xmax; + double ymin, ymax; + + //Radar image related inputs + double startingX, startingY; + double XSize, YSize; + int nLines, nPixels; + double dt; + int nodata_out; + int pOff, lOff, pCount, lCount; + double X_res, Y_res; + + //dt-varying search range rountine parameters + double dt_unity; + double max_factor; + double upper_thld, lower_thld; + + //Output file names + std::string pixlinename; + std::string offsetname; + std::string searchrangename; + std::string chipsizeminname; + std::string chipsizemaxname; + std::string stablesurfacemaskname; + std::string ro2vx_name; + std::string ro2vy_name; + + //Functions + void computeBbox(double *); + void geogridOptical(); + void cross_C(double r_u[3], double r_v[3], double r_w[3]); + double dot_C(double r_v[3], double r_w[3]); + double norm_C(double a[3]); + void unitvec_C(double v[3], double u[3]); +}; + + +#endif diff --git a/geo_autoRIFT/geogrid/include/geogridOpticalmodule.h b/geo_autoRIFT/geogrid/include/geogridOpticalmodule.h new file mode 100755 index 0000000..5ceb7ef --- /dev/null +++ b/geo_autoRIFT/geogrid/include/geogridOpticalmodule.h @@ -0,0 +1,124 @@ +/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + * Copyright 2019 California Institute of Technology. ALL RIGHTS RESERVED. + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + * United States Government Sponsorship acknowledged. This software is subject to + * U.S. export control laws and regulations and has been classified as 'EAR99 NLR' + * (No [Export] License Required except when exporting to an embargoed country, + * end user, or in support of a prohibited end use). By downloading this software, + * the user agrees to comply with all applicable U.S. export laws and regulations. + * The user has the responsibility to obtain export licenses, or other export + * authority as may be required before exporting this software to any 'EAR99' + * embargoed foreign country or citizen of those countries. + * + * Authors: Piyush Agram, Yang Lei + *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + */ + + +#ifndef geogridOpticalmodule_h +#define geogridOpticalmodule_h + +#include +#include + +extern "C" +{ + PyObject * createGeoGridOptical(PyObject*, PyObject*); + PyObject * destroyGeoGridOptical(PyObject*, PyObject*); + PyObject * geogridOptical(PyObject *, PyObject *); + PyObject * setOpticalImageDimensions(PyObject *, PyObject *); + PyObject * setXParameters(PyObject *, PyObject *); + PyObject * setYParameters(PyObject*, PyObject *); + PyObject * setRepeatTime(PyObject *, PyObject *); + + PyObject * setDEM(PyObject *, PyObject *); + PyObject * setVelocities(PyObject*, PyObject*); + PyObject * setSearchRange(PyObject*, PyObject*); + PyObject * setChipSizeMin(PyObject*, PyObject*); + PyObject * setChipSizeMax(PyObject*, PyObject*); + PyObject * setStableSurfaceMask(PyObject*, PyObject*); + PyObject * setSlopes(PyObject*, PyObject*); + PyObject * setNodataOut(PyObject *, PyObject *); + + PyObject * setDtUnity(PyObject *, PyObject *); + PyObject * setMaxFactor(PyObject *, PyObject *); + PyObject * setUpperThreshold(PyObject*, PyObject *); + PyObject * setLowerThreshold(PyObject *, PyObject *); + + PyObject * setWindowLocationsFilename(PyObject *, PyObject *); + PyObject * setWindowOffsetsFilename(PyObject *, PyObject *); + PyObject * setWindowSearchRangeFilename(PyObject *, PyObject *); + PyObject * setWindowChipSizeMinFilename(PyObject *, PyObject *); + PyObject * setWindowChipSizeMaxFilename(PyObject *, PyObject *); + PyObject * setWindowStableSurfaceMaskFilename(PyObject *, PyObject *); + PyObject * setRO2VXFilename(PyObject *, PyObject *); + PyObject * setRO2VYFilename(PyObject *, PyObject *); + PyObject * setEPSG(PyObject *, PyObject *); + PyObject * setChipSizeX0(PyObject *, PyObject *); + PyObject * setGridSpacingX(PyObject *, PyObject *); + PyObject * setXLimits(PyObject *, PyObject *); + PyObject * setYLimits(PyObject *, PyObject *); + PyObject * getXPixelSize(PyObject *, PyObject *); + PyObject * getYPixelSize(PyObject *, PyObject *); + PyObject * getXOff(PyObject *, PyObject *); + PyObject * getYOff(PyObject *, PyObject *); + PyObject * getXCount(PyObject *, PyObject *); + PyObject * getYCount(PyObject *, PyObject *); +} + +static PyMethodDef geogrid_methods[] = +{ + {"createGeoGridOptical_Py", createGeoGridOptical, METH_VARARGS, " "}, + {"destroyGeoGridOptical_Py", destroyGeoGridOptical, METH_VARARGS, " "}, + {"geogridOptical_Py", geogridOptical, METH_VARARGS, " "}, + {"setOpticalImageDimensions_Py", setOpticalImageDimensions, METH_VARARGS, " "}, + {"setXParameters_Py", setXParameters, METH_VARARGS, " "}, + {"setYParameters_Py", setYParameters, METH_VARARGS, " "}, + {"setRepeatTime_Py", setRepeatTime, METH_VARARGS, " "}, + {"setDEM_Py", setDEM, METH_VARARGS, " "}, + {"setEPSG_Py", setEPSG, METH_VARARGS, " "}, + {"setChipSizeX0_Py", setChipSizeX0, METH_VARARGS, " "}, + {"setGridSpacingX_Py", setGridSpacingX, METH_VARARGS, " "}, + {"setVelocities_Py", setVelocities, METH_VARARGS, " "}, + {"setSearchRange_Py", setSearchRange, METH_VARARGS, " "}, + {"setChipSizeMin_Py", setChipSizeMin, METH_VARARGS, " "}, + {"setChipSizeMax_Py", setChipSizeMax, METH_VARARGS, " "}, + {"setStableSurfaceMask_Py", setStableSurfaceMask, METH_VARARGS, " "}, + {"setSlopes_Py", setSlopes, METH_VARARGS, " "}, + {"setNodataOut_Py", setNodataOut, METH_VARARGS, " "}, + {"setDtUnity_Py", setDtUnity, METH_VARARGS, " "}, + {"setMaxFactor_Py", setMaxFactor, METH_VARARGS, " "}, + {"setUpperThreshold_Py", setUpperThreshold, METH_VARARGS, " "}, + {"setLowerThreshold_Py", setLowerThreshold, METH_VARARGS, " "}, + {"setXLimits_Py", setXLimits, METH_VARARGS, " "}, + {"setYLimits_Py", setYLimits, METH_VARARGS, " "}, + {"getXPixelSize_Py", getXPixelSize, METH_VARARGS, " "}, + {"getYPixelSize_Py", getYPixelSize, METH_VARARGS, " "}, + {"getXOff_Py", getXOff, METH_VARARGS, " "}, + {"getYOff_Py", getYOff, METH_VARARGS, " "}, + {"getXCount_Py", getXCount, METH_VARARGS, " "}, + {"getYCount_Py", getYCount, METH_VARARGS, " "}, + {"setWindowLocationsFilename_Py", setWindowLocationsFilename, METH_VARARGS, " "}, + {"setWindowOffsetsFilename_Py", setWindowOffsetsFilename, METH_VARARGS, " "}, + {"setWindowSearchRangeFilename_Py", setWindowSearchRangeFilename, METH_VARARGS, " "}, + {"setWindowChipSizeMinFilename_Py", setWindowChipSizeMinFilename, METH_VARARGS, " "}, + {"setWindowChipSizeMaxFilename_Py", setWindowChipSizeMaxFilename, METH_VARARGS, " "}, + {"setWindowStableSurfaceMaskFilename_Py", setWindowStableSurfaceMaskFilename, METH_VARARGS, " "}, + {"setRO2VXFilename_Py", setRO2VXFilename, METH_VARARGS, " "}, + {"setRO2VYFilename_Py", setRO2VYFilename, METH_VARARGS, " "}, + {NULL, NULL, 0, NULL} +}; +#endif //geoGridOpticalmodule_h + diff --git a/geo_autoRIFT/geogrid/src/SConscript1 b/geo_autoRIFT/geogrid/src/SConscript1 new file mode 100755 index 0000000..104fc03 --- /dev/null +++ b/geo_autoRIFT/geogrid/src/SConscript1 @@ -0,0 +1,24 @@ +#!/usr/bin/env python +#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +# +# Author: Piyush Agram +# Copyright 2019, by the California Institute of Technology. ALL RIGHTS RESERVED. United States Government Sponsorship acknowledged. +# Any commercial use must be negotiated with the Office of Technology Transfer at the California Institute of Technology. +# +# This software may be subject to U.S. export control laws. By accepting this software, the user agrees to comply with all applicable U.S. +# export laws and regulations. User has the responsibility to obtain export licenses, or other export authority as may be required before +# exporting such information to foreign countries or providing access to foreign persons. +# +#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + + +import os + +Import('envgeogrid') +build = envgeogrid['PRJ_LIB_DIR'] +envgeogrid.AppendUnique(CXXFLAGS = '-fopenmp') +envgeogrid.AppendUnique(CXXFLAGS = '-std=c++11') +listFiles = ['geogridOptical.cpp'] +lib = envgeogrid.Library(target = 'geogridOptical', source = listFiles) +envgeogrid.Install(build,lib) +envgeogrid.Alias('build',build) diff --git a/geo_autoRIFT/geogrid/src/geogrid.cpp b/geo_autoRIFT/geogrid/src/geogrid.cpp index 1d7f06d..a859816 100755 --- a/geo_autoRIFT/geogrid/src/geogrid.cpp +++ b/geo_autoRIFT/geogrid/src/geogrid.cpp @@ -1265,7 +1265,7 @@ void geoGrid::geogrid() if (srxname != "") { - if ((vxname != "")&(vel[0] == nodata)) + if ((schrng1[0] == nodata)|(schrng1[0] == 0)) { sr_raster11[jj] = 0; sr_raster22[jj] = 0; @@ -1299,21 +1299,44 @@ void geoGrid::geogrid() if (csminxname != "") { - csmin_raster11[jj] = csminxLine[jj] / ChipSizeX0 * ChipSizeX0_PIX_grd; - csmin_raster22[jj] = csminyLine[jj] / ChipSizeX0 * ChipSizeX0_PIX_azm; + if (csminxLine[jj] == nodata) + { + csmin_raster11[jj] = nodata_out; + csmin_raster22[jj] = nodata_out; + } + else + { + csmin_raster11[jj] = csminxLine[jj] / ChipSizeX0 * ChipSizeX0_PIX_grd; + csmin_raster22[jj] = csminyLine[jj] / ChipSizeX0 * ChipSizeX0_PIX_azm; + } } if (csmaxxname != "") { - csmax_raster11[jj] = csmaxxLine[jj] / ChipSizeX0 * ChipSizeX0_PIX_grd; - csmax_raster22[jj] = csmaxyLine[jj] / ChipSizeX0 * ChipSizeX0_PIX_azm; + if (csmaxxLine[jj] == nodata) + { + csmax_raster11[jj] = nodata_out; + csmax_raster22[jj] = nodata_out; + } + else + { + csmax_raster11[jj] = csmaxxLine[jj] / ChipSizeX0 * ChipSizeX0_PIX_grd; + csmax_raster22[jj] = csmaxyLine[jj] / ChipSizeX0 * ChipSizeX0_PIX_azm; + } } if (ssmname != "") { - ssm_raster[jj] = ssmLine[jj]; + if (ssmLine[jj] == nodata) + { + ssm_raster[jj] = nodata_out; + } + else + { + ssm_raster[jj] = ssmLine[jj]; + } } diff --git a/geo_autoRIFT/geogrid/src/geogridOptical.cpp b/geo_autoRIFT/geogrid/src/geogridOptical.cpp new file mode 100755 index 0000000..0a2f9db --- /dev/null +++ b/geo_autoRIFT/geogrid/src/geogridOptical.cpp @@ -0,0 +1,1347 @@ +/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + * Copyright 2019 California Institute of Technology. ALL RIGHTS RESERVED. + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + * United States Government Sponsorship acknowledged. This software is subject to + * U.S. export control laws and regulations and has been classified as 'EAR99 NLR' + * (No [Export] License Required except when exporting to an embargoed country, + * end user, or in support of a prohibited end use). By downloading this software, + * the user agrees to comply with all applicable U.S. export laws and regulations. + * The user has the responsibility to obtain export licenses, or other export + * authority as may be required before exporting this software to any 'EAR99' + * embargoed foreign country or citizen of those countries. + * + * Authors: Piyush Agram, Yang Lei + *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + */ + +#include "geogridOptical.h" +#include +#include +#include +#include +#include + + + + +void geoGridOptical::geogridOptical() +{ + //Some constants + double deg2rad = M_PI/180.0; + + //For now print inputs that were obtained + + std::cout << "\nOptical Image parameters: \n"; + std::cout << "X-direction coordinate: " << startingX << " " << XSize << "\n"; + std::cout << "Y-direction coordinate: " << startingY << " " << YSize << "\n"; + std::cout << "Dimensions: " << nPixels << " " << nLines << "\n"; + + std::cout << "\nMap inputs: \n"; + std::cout << "EPSG: " << epsgDem << "\n"; + std::cout << "Smallest Allowable Chip Size in m: " << chipSizeX0 << "\n"; + std::cout << "Grid spacing in m: " << gridSpacingX << "\n"; + std::cout << "Repeat Time: " << dt << "\n"; + std::cout << "XLimits: " << xmin << " " << xmax << "\n"; + std::cout << "YLimits: " << ymin << " " << ymax << "\n"; + std::cout << "Extent in km: " << (xmax - xmin)/1000. << " " << (ymax - ymin)/1000. << "\n"; + if (demname != "") + { + std::cout << "DEM: " << demname << "\n"; + } + if (dhdxname != "") + { + std::cout << "Slopes: " << dhdxname << " " << dhdyname << "\n"; + } + if (vxname != "") + { + std::cout << "Velocities: " << vxname << " " << vyname << "\n"; + } + if (srxname != "") + { + std::cout << "Search Range: " << srxname << " " << sryname << "\n"; + } + if (csminxname != "") + { + std::cout << "Chip Size Min: " << csminxname << " " << csminyname << "\n"; + } + if (csmaxxname != "") + { + std::cout << "Chip Size Max: " << csmaxxname << " " << csmaxyname << "\n"; + } + if (ssmname != "") + { + std::cout << "Stable Surface Mask: " << ssmname << "\n"; + } + + + std::cout << "\nOutputs: \n"; + std::cout << "Window locations: " << pixlinename << "\n"; + if (dhdxname != "") + { + if (vxname != "") + { + std::cout << "Window offsets: " << offsetname << "\n"; + } + + std::cout << "Window rdr_off2vel_x vector: " << ro2vx_name << "\n"; + std::cout << "Window rdr_off2vel_y vector: " << ro2vy_name << "\n"; + + if (srxname != "") + { + std::cout << "Window search range: " << searchrangename << "\n"; + } + } + + if (csminxname != "") + { + std::cout << "Window chip size min: " << chipsizeminname << "\n"; + } + if (csmaxxname != "") + { + std::cout << "Window chip size max: " << chipsizemaxname << "\n"; + } + if (ssmname != "") + { + std::cout << "Window stable surface mask: " << stablesurfacemaskname << "\n"; + } + + std::cout << "Output Nodata Value: " << nodata_out << "\n"; + + + std::cout << "\nStarting processing .... \n"; + + //Startup GDAL + GDALAllRegister(); + + //DEM related information + GDALDataset* demDS = NULL; + GDALDataset* sxDS = NULL; + GDALDataset* syDS = NULL; + GDALDataset* vxDS = NULL; + GDALDataset* vyDS = NULL; + GDALDataset* srxDS = NULL; + GDALDataset* sryDS = NULL; + GDALDataset* csminxDS = NULL; + GDALDataset* csminyDS = NULL; + GDALDataset* csmaxxDS = NULL; + GDALDataset* csmaxyDS = NULL; + GDALDataset* ssmDS = NULL; + + double geoTrans[6]; + + demDS = reinterpret_cast(GDALOpenShared(demname.c_str(), GA_ReadOnly)); + if (dhdxname != "") + { + sxDS = reinterpret_cast(GDALOpenShared(dhdxname.c_str(), GA_ReadOnly)); + syDS = reinterpret_cast(GDALOpenShared(dhdyname.c_str(), GA_ReadOnly)); + } + if (vxname != "") + { + vxDS = reinterpret_cast(GDALOpenShared(vxname.c_str(), GA_ReadOnly)); + vyDS = reinterpret_cast(GDALOpenShared(vyname.c_str(), GA_ReadOnly)); + } + if (srxname != "") + { + srxDS = reinterpret_cast(GDALOpenShared(srxname.c_str(), GA_ReadOnly)); + sryDS = reinterpret_cast(GDALOpenShared(sryname.c_str(), GA_ReadOnly)); + } + if (csminxname != "") + { + csminxDS = reinterpret_cast(GDALOpenShared(csminxname.c_str(), GA_ReadOnly)); + csminyDS = reinterpret_cast(GDALOpenShared(csminyname.c_str(), GA_ReadOnly)); + } + if (csmaxxname != "") + { + csmaxxDS = reinterpret_cast(GDALOpenShared(csmaxxname.c_str(), GA_ReadOnly)); + csmaxyDS = reinterpret_cast(GDALOpenShared(csmaxyname.c_str(), GA_ReadOnly)); + } + if (ssmname != "") + { + ssmDS = reinterpret_cast(GDALOpenShared(ssmname.c_str(), GA_ReadOnly)); + } + if (demDS == NULL) + { + std::cout << "Error opening DEM file { " << demname << " }\n"; + std::cout << "Exiting with error code .... (101) \n"; + GDALDestroyDriverManager(); + exit(101); + } + if (dhdxname != "") + { + if (sxDS == NULL) + { + std::cout << "Error opening x-direction slope file { " << dhdxname << " }\n"; + std::cout << "Exiting with error code .... (101) \n"; + GDALDestroyDriverManager(); + exit(101); + } + if (syDS == NULL) + { + std::cout << "Error opening y-direction slope file { " << dhdyname << " }\n"; + std::cout << "Exiting with error code .... (101) \n"; + GDALDestroyDriverManager(); + exit(101); + } + } + if (vxname != "") + { + if (vxDS == NULL) + { + std::cout << "Error opening x-direction velocity file { " << vxname << " }\n"; + std::cout << "Exiting with error code .... (101) \n"; + GDALDestroyDriverManager(); + exit(101); + } + if (vyDS == NULL) + { + std::cout << "Error opening y-direction velocity file { " << vyname << " }\n"; + std::cout << "Exiting with error code .... (101) \n"; + GDALDestroyDriverManager(); + exit(101); + } + } + if (srxname != "") + { + if (srxDS == NULL) + { + std::cout << "Error opening x-direction search range file { " << srxname << " }\n"; + std::cout << "Exiting with error code .... (101) \n"; + GDALDestroyDriverManager(); + exit(101); + } + if (sryDS == NULL) + { + std::cout << "Error opening y-direction search range file { " << sryname << " }\n"; + std::cout << "Exiting with error code .... (101) \n"; + GDALDestroyDriverManager(); + exit(101); + } + } + if (csminxname != "") + { + if (csminxDS == NULL) + { + std::cout << "Error opening x-direction chip size min file { " << csminxname << " }\n"; + std::cout << "Exiting with error code .... (101) \n"; + GDALDestroyDriverManager(); + exit(101); + } + if (csminyDS == NULL) + { + std::cout << "Error opening y-direction chip size min file { " << csminyname << " }\n"; + std::cout << "Exiting with error code .... (101) \n"; + GDALDestroyDriverManager(); + exit(101); + } + } + if (csmaxxname != "") + { + if (csmaxxDS == NULL) + { + std::cout << "Error opening x-direction chip size max file { " << csmaxxname << " }\n"; + std::cout << "Exiting with error code .... (101) \n"; + GDALDestroyDriverManager(); + exit(101); + } + if (csmaxyDS == NULL) + { + std::cout << "Error opening y-direction chip size max file { " << csmaxyname << " }\n"; + std::cout << "Exiting with error code .... (101) \n"; + GDALDestroyDriverManager(); + exit(101); + } + } + if (ssmname != "") + { + if (ssmDS == NULL) + { + std::cout << "Error opening stable surface mask file { " << ssmname << " }\n"; + std::cout << "Exiting with error code .... (101) \n"; + GDALDestroyDriverManager(); + exit(101); + } + } + + demDS->GetGeoTransform(geoTrans); + int demXSize = demDS->GetRasterXSize(); + int demYSize = demDS->GetRasterYSize(); + + + //Get offsets and size to read from DEM +// int lOff = std::max( std::floor((ymax - geoTrans[3])/geoTrans[5]), 0.); +// int lCount = std::min( std::ceil((ymin - geoTrans[3])/geoTrans[5]), demYSize-1.) - lOff; +// +// int pOff = std::max( std::floor((xmin - geoTrans[0])/geoTrans[1]), 0.); +// int pCount = std::min( std::ceil((xmax - geoTrans[0])/geoTrans[1]), demXSize-1.) - pOff; + lOff = std::max( std::floor((ymax - geoTrans[3])/geoTrans[5]), 0.); + lCount = std::min( std::ceil((ymin - geoTrans[3])/geoTrans[5]), demYSize-1.) - lOff; + + pOff = std::max( std::floor((xmin - geoTrans[0])/geoTrans[1]), 0.); + pCount = std::min( std::ceil((xmax - geoTrans[0])/geoTrans[1]), demXSize-1.) - pOff; + + + std::cout << "Xlimits : " << geoTrans[0] + pOff * geoTrans[1] << " " + << geoTrans[0] + (pOff + pCount) * geoTrans[1] << "\n"; + + + std::cout << "Ylimits : " << geoTrans[3] + (lOff + lCount) * geoTrans[5] << " " + << geoTrans[3] + lOff * geoTrans[5] << "\n"; + + std::cout << "Origin index (in DEM) of geogrid: " << pOff << " " << lOff << "\n"; + + std::cout << "Dimensions of geogrid: " << pCount << " x " << lCount << "\n"; + + + //Create GDAL Transformers + OGRSpatialReference demSRS(nullptr); + if (demSRS.importFromEPSG(epsgDem) != 0) + { + std::cout << "Could not create OGR spatial reference for DEM EPSG code: " << epsgDem << "\n"; + GDALClose(demDS); + GDALDestroyDriverManager(); + exit(102); + } + + OGRSpatialReference datSRS(nullptr); + if (datSRS.importFromEPSG(epsgDat) != 0) + { + std::cout << "Could not create OGR spatil reference for Data EPSG code: " << epsgDat << "\n"; + exit(103); + } + + OGRCoordinateTransformation *fwdTrans = OGRCreateCoordinateTransformation( &demSRS, &datSRS); + OGRCoordinateTransformation *invTrans = OGRCreateCoordinateTransformation( &datSRS, &demSRS); + + + + std::vector demLine(pCount); + std::vector sxLine(pCount); + std::vector syLine(pCount); + std::vector vxLine(pCount); + std::vector vyLine(pCount); + std::vector srxLine(pCount); + std::vector sryLine(pCount); + std::vector csminxLine(pCount); + std::vector csminyLine(pCount); + std::vector csmaxxLine(pCount); + std::vector csmaxyLine(pCount); + std::vector ssmLine(pCount); + + GInt32 raster1[pCount]; + GInt32 raster2[pCount]; + GInt32 raster11[pCount]; + GInt32 raster22[pCount]; + + GInt32 sr_raster11[pCount]; + GInt32 sr_raster22[pCount]; + GInt32 csmin_raster11[pCount]; + GInt32 csmin_raster22[pCount]; + GInt32 csmax_raster11[pCount]; + GInt32 csmax_raster22[pCount]; + GInt32 ssm_raster[pCount]; + + double raster1a[pCount]; + double raster1b[pCount]; +// double raster1c[pCount]; + + double raster2a[pCount]; + double raster2b[pCount]; +// double raster2c[pCount]; + + GDALRasterBand *poBand1 = NULL; + GDALRasterBand *poBand2 = NULL; + GDALRasterBand *poBand1Off = NULL; + GDALRasterBand *poBand2Off = NULL; + GDALRasterBand *poBand1Sch = NULL; + GDALRasterBand *poBand2Sch = NULL; + GDALRasterBand *poBand1Min = NULL; + GDALRasterBand *poBand2Min = NULL; + GDALRasterBand *poBand1Max = NULL; + GDALRasterBand *poBand2Max = NULL; + GDALRasterBand *poBand1Msk = NULL; + GDALRasterBand *poBand1RO2VX = NULL; + GDALRasterBand *poBand1RO2VY = NULL; + GDALRasterBand *poBand2RO2VX = NULL; + GDALRasterBand *poBand2RO2VY = NULL; + + GDALDataset *poDstDS = NULL; + GDALDataset *poDstDSOff = NULL; + GDALDataset *poDstDSSch = NULL; + GDALDataset *poDstDSMin = NULL; + GDALDataset *poDstDSMax = NULL; + GDALDataset *poDstDSMsk = NULL; + GDALDataset *poDstDSRO2VX = NULL; + GDALDataset *poDstDSRO2VY = NULL; + + + double nodata; + int* pbSuccess = NULL; + nodata = demDS->GetRasterBand(1)->GetNoDataValue(pbSuccess); + + + const char *pszFormat = "GTiff"; + char **papszOptions = NULL; + std::string str = ""; + double adfGeoTransform[6] = { geoTrans[0] + pOff * geoTrans[1], geoTrans[1], 0, geoTrans[3] + lOff * geoTrans[5], 0, geoTrans[5]}; + OGRSpatialReference oSRS; + char *pszSRS_WKT = NULL; + demSRS.exportToWkt( &pszSRS_WKT ); + + + + GDALDriver *poDriver; + poDriver = GetGDALDriverManager()->GetDriverByName(pszFormat); + if( poDriver == NULL ) + exit(107); +// GDALDataset *poDstDS; + + str = pixlinename; + const char * pszDstFilename = str.c_str(); + poDstDS = poDriver->Create( pszDstFilename, pCount, lCount, 2, GDT_Int32, + papszOptions ); + + + poDstDS->SetGeoTransform( adfGeoTransform ); + poDstDS->SetProjection( pszSRS_WKT ); +// CPLFree( pszSRS_WKT ); + + +// GDALRasterBand *poBand1; +// GDALRasterBand *poBand2; + poBand1 = poDstDS->GetRasterBand(1); + poBand2 = poDstDS->GetRasterBand(2); + poBand1->SetNoDataValue(nodata_out); + poBand2->SetNoDataValue(nodata_out); + + + if ((dhdxname != "")&(vxname != "")) + { + + GDALDriver *poDriverOff; + poDriverOff = GetGDALDriverManager()->GetDriverByName(pszFormat); + if( poDriverOff == NULL ) + exit(107); +// GDALDataset *poDstDSOff; + + str = offsetname; + const char * pszDstFilenameOff = str.c_str(); + poDstDSOff = poDriverOff->Create( pszDstFilenameOff, pCount, lCount, 2, GDT_Int32, + papszOptions ); + + poDstDSOff->SetGeoTransform( adfGeoTransform ); + poDstDSOff->SetProjection( pszSRS_WKT ); + // CPLFree( pszSRS_WKT ); + +// GDALRasterBand *poBand1Off; +// GDALRasterBand *poBand2Off; + poBand1Off = poDstDSOff->GetRasterBand(1); + poBand2Off = poDstDSOff->GetRasterBand(2); + poBand1Off->SetNoDataValue(nodata_out); + poBand2Off->SetNoDataValue(nodata_out); + + } + + if ((dhdxname != "")&(srxname != "")) + { + + GDALDriver *poDriverSch; + poDriverSch = GetGDALDriverManager()->GetDriverByName(pszFormat); + if( poDriverSch == NULL ) + exit(107); +// GDALDataset *poDstDSSch; + + str = searchrangename; + const char * pszDstFilenameSch = str.c_str(); + poDstDSSch = poDriverSch->Create( pszDstFilenameSch, pCount, lCount, 2, GDT_Int32, + papszOptions ); + + poDstDSSch->SetGeoTransform( adfGeoTransform ); + poDstDSSch->SetProjection( pszSRS_WKT ); + // CPLFree( pszSRS_WKT ); + +// GDALRasterBand *poBand1Sch; +// GDALRasterBand *poBand2Sch; + poBand1Sch = poDstDSSch->GetRasterBand(1); + poBand2Sch = poDstDSSch->GetRasterBand(2); + poBand1Sch->SetNoDataValue(nodata_out); + poBand2Sch->SetNoDataValue(nodata_out); + + } + + if (csminxname != "") + { + + GDALDriver *poDriverMin; + poDriverMin = GetGDALDriverManager()->GetDriverByName(pszFormat); + if( poDriverMin == NULL ) + exit(107); +// GDALDataset *poDstDSMin; + + str = chipsizeminname; + const char * pszDstFilenameMin = str.c_str(); + poDstDSMin = poDriverMin->Create( pszDstFilenameMin, pCount, lCount, 2, GDT_Int32, + papszOptions ); + + poDstDSMin->SetGeoTransform( adfGeoTransform ); + poDstDSMin->SetProjection( pszSRS_WKT ); + // CPLFree( pszSRS_WKT ); + +// GDALRasterBand *poBand1Min; +// GDALRasterBand *poBand2Min; + poBand1Min = poDstDSMin->GetRasterBand(1); + poBand2Min = poDstDSMin->GetRasterBand(2); + poBand1Min->SetNoDataValue(nodata_out); + poBand2Min->SetNoDataValue(nodata_out); + + } + + + if (csmaxxname != "") + { + + GDALDriver *poDriverMax; + poDriverMax = GetGDALDriverManager()->GetDriverByName(pszFormat); + if( poDriverMax == NULL ) + exit(107); +// GDALDataset *poDstDSMax; + + str = chipsizemaxname; + const char * pszDstFilenameMax = str.c_str(); + poDstDSMax = poDriverMax->Create( pszDstFilenameMax, pCount, lCount, 2, GDT_Int32, + papszOptions ); + + poDstDSMax->SetGeoTransform( adfGeoTransform ); + poDstDSMax->SetProjection( pszSRS_WKT ); + // CPLFree( pszSRS_WKT ); + +// GDALRasterBand *poBand1Max; +// GDALRasterBand *poBand2Max; + poBand1Max = poDstDSMax->GetRasterBand(1); + poBand2Max = poDstDSMax->GetRasterBand(2); + poBand1Max->SetNoDataValue(nodata_out); + poBand2Max->SetNoDataValue(nodata_out); + + } + + + + if (ssmname != "") + { + + GDALDriver *poDriverMsk; + poDriverMsk = GetGDALDriverManager()->GetDriverByName(pszFormat); + if( poDriverMsk == NULL ) + exit(107); +// GDALDataset *poDstDSMsk; + + str = stablesurfacemaskname; + const char * pszDstFilenameMsk = str.c_str(); + poDstDSMsk = poDriverMsk->Create( pszDstFilenameMsk, pCount, lCount, 1, GDT_Int32, + papszOptions ); + + poDstDSMsk->SetGeoTransform( adfGeoTransform ); + poDstDSMsk->SetProjection( pszSRS_WKT ); + // CPLFree( pszSRS_WKT ); + +// GDALRasterBand *poBand1Msk; + poBand1Msk = poDstDSMsk->GetRasterBand(1); + poBand1Msk->SetNoDataValue(nodata_out); + + } + + + if (dhdxname != "") + { + + GDALDriver *poDriverRO2VX; + poDriverRO2VX = GetGDALDriverManager()->GetDriverByName(pszFormat); + if( poDriverRO2VX == NULL ) + exit(107); +// GDALDataset *poDstDSRO2VX; + + str = ro2vx_name; + const char * pszDstFilenameRO2VX = str.c_str(); + poDstDSRO2VX = poDriverRO2VX->Create( pszDstFilenameRO2VX, pCount, lCount, 2, GDT_Float64, + papszOptions ); + + poDstDSRO2VX->SetGeoTransform( adfGeoTransform ); + poDstDSRO2VX->SetProjection( pszSRS_WKT ); + // CPLFree( pszSRS_WKT ); + +// GDALRasterBand *poBand1RO2VX; +// GDALRasterBand *poBand2RO2VX; + // GDALRasterBand *poBand3Los; + poBand1RO2VX = poDstDSRO2VX->GetRasterBand(1); + poBand2RO2VX = poDstDSRO2VX->GetRasterBand(2); + // poBand3Los = poDstDSLos->GetRasterBand(3); + poBand1RO2VX->SetNoDataValue(nodata_out); + poBand2RO2VX->SetNoDataValue(nodata_out); + // poBand3Los->SetNoDataValue(nodata_out); + + + GDALDriver *poDriverRO2VY; + poDriverRO2VY = GetGDALDriverManager()->GetDriverByName(pszFormat); + if( poDriverRO2VY == NULL ) + exit(107); +// GDALDataset *poDstDSRO2VY; + + str = ro2vy_name; + const char * pszDstFilenameRO2VY = str.c_str(); + poDstDSRO2VY = poDriverRO2VY->Create( pszDstFilenameRO2VY, pCount, lCount, 2, GDT_Float64, + papszOptions ); + + poDstDSRO2VY->SetGeoTransform( adfGeoTransform ); + poDstDSRO2VY->SetProjection( pszSRS_WKT ); +// CPLFree( pszSRS_WKT ); + +// GDALRasterBand *poBand1RO2VY; +// GDALRasterBand *poBand2RO2VY; + // GDALRasterBand *poBand3Alt; + poBand1RO2VY = poDstDSRO2VY->GetRasterBand(1); + poBand2RO2VY = poDstDSRO2VY->GetRasterBand(2); + // poBand3Alt = poDstDSAlt->GetRasterBand(3); + poBand1RO2VY->SetNoDataValue(nodata_out); + poBand2RO2VY->SetNoDataValue(nodata_out); + // poBand3Alt->SetNoDataValue(nodata_out); + + } + + CPLFree( pszSRS_WKT ); + + + + + + // ground range and azimuth pixel size + + X_res = std::abs(XSize); + Y_res = std::abs(YSize); + std::cout << "X-direction pixel size: " << X_res << "\n"; + std::cout << "Y-direction pixel size: " << Y_res << "\n"; +// int ChipSizeX0 = 240; + double ChipSizeX0 = chipSizeX0; + int ChipSizeX0_PIX_X = std::ceil(ChipSizeX0 / X_res / 4) * 4; + int ChipSizeX0_PIX_Y = std::ceil(ChipSizeX0 / Y_res / 4) * 4; + + double xind, yind; + + for (int ii=0; iiGetRasterBand(1)->RasterIO(GF_Read, + pOff, lOff+ii, + pCount, 1, + (void*) (demLine.data()), + pCount, 1, GDT_Float64, + sizeof(double), sizeof(double)*pCount, NULL); + + if (status != 0) + { + std::cout << "Error read line " << lOff + ii << " from DEM file: " << demname << "\n"; + GDALClose(demDS); + GDALDestroyDriverManager(); + exit(105); + } + + if (dhdxname != "") + { + status = sxDS->GetRasterBand(1)->RasterIO(GF_Read, + pOff, lOff+ii, + pCount, 1, + (void*) (sxLine.data()), + pCount, 1, GDT_Float64, + sizeof(double), sizeof(double)*pCount, NULL); + + if (status != 0) + { + std::cout << "Error read line " << lOff + ii << " from x-direction slope file: " << dhdxname << "\n"; + GDALClose(sxDS); + GDALDestroyDriverManager(); + exit(105); + } + + + status = syDS->GetRasterBand(1)->RasterIO(GF_Read, + pOff, lOff+ii, + pCount, 1, + (void*) (syLine.data()), + pCount, 1, GDT_Float64, + sizeof(double), sizeof(double)*pCount, NULL); + + if (status != 0) + { + std::cout << "Error read line " << lOff + ii << " from y-direction slope file: " << dhdyname << "\n"; + GDALClose(syDS); + GDALDestroyDriverManager(); + exit(105); + } + + } + + if (vxname != "") + { + status = vxDS->GetRasterBand(1)->RasterIO(GF_Read, + pOff, lOff+ii, + pCount, 1, + (void*) (vxLine.data()), + pCount, 1, GDT_Float64, + sizeof(double), sizeof(double)*pCount, NULL); + + if (status != 0) + { + std::cout << "Error read line " << lOff + ii << " from x-direction velocity file: " << vxname << "\n"; + GDALClose(vxDS); + GDALDestroyDriverManager(); + exit(105); + } + + + status = vyDS->GetRasterBand(1)->RasterIO(GF_Read, + pOff, lOff+ii, + pCount, 1, + (void*) (vyLine.data()), + pCount, 1, GDT_Float64, + sizeof(double), sizeof(double)*pCount, NULL); + + if (status != 0) + { + std::cout << "Error read line " << lOff + ii << " from y-direction velocity file: " << vyname << "\n"; + GDALClose(vyDS); + GDALDestroyDriverManager(); + exit(105); + } + } + + if (srxname != "") + { + status = srxDS->GetRasterBand(1)->RasterIO(GF_Read, + pOff, lOff+ii, + pCount, 1, + (void*) (srxLine.data()), + pCount, 1, GDT_Float64, + sizeof(double), sizeof(double)*pCount, NULL); + + if (status != 0) + { + std::cout << "Error read line " << lOff + ii << " from x-direction search range file: " << srxname << "\n"; + GDALClose(srxDS); + GDALDestroyDriverManager(); + exit(105); + } + + + status = sryDS->GetRasterBand(1)->RasterIO(GF_Read, + pOff, lOff+ii, + pCount, 1, + (void*) (sryLine.data()), + pCount, 1, GDT_Float64, + sizeof(double), sizeof(double)*pCount, NULL); + + if (status != 0) + { + std::cout << "Error read line " << lOff + ii << " from y-direction search range file: " << sryname << "\n"; + GDALClose(sryDS); + GDALDestroyDriverManager(); + exit(105); + } + } + + + if (csminxname != "") + { + status = csminxDS->GetRasterBand(1)->RasterIO(GF_Read, + pOff, lOff+ii, + pCount, 1, + (void*) (csminxLine.data()), + pCount, 1, GDT_Float64, + sizeof(double), sizeof(double)*pCount, NULL); + + if (status != 0) + { + std::cout << "Error read line " << lOff + ii << " from x-direction chip size min file: " << csminxname << "\n"; + GDALClose(csminxDS); + GDALDestroyDriverManager(); + exit(105); + } + + + status = csminyDS->GetRasterBand(1)->RasterIO(GF_Read, + pOff, lOff+ii, + pCount, 1, + (void*) (csminyLine.data()), + pCount, 1, GDT_Float64, + sizeof(double), sizeof(double)*pCount, NULL); + + if (status != 0) + { + std::cout << "Error read line " << lOff + ii << " from y-direction chip size min file: " << csminyname << "\n"; + GDALClose(csminyDS); + GDALDestroyDriverManager(); + exit(105); + } + } + + + if (csmaxxname != "") + { + status = csmaxxDS->GetRasterBand(1)->RasterIO(GF_Read, + pOff, lOff+ii, + pCount, 1, + (void*) (csmaxxLine.data()), + pCount, 1, GDT_Float64, + sizeof(double), sizeof(double)*pCount, NULL); + + if (status != 0) + { + std::cout << "Error read line " << lOff + ii << " from x-direction chip size max file: " << csmaxxname << "\n"; + GDALClose(csmaxxDS); + GDALDestroyDriverManager(); + exit(105); + } + + + status = csmaxyDS->GetRasterBand(1)->RasterIO(GF_Read, + pOff, lOff+ii, + pCount, 1, + (void*) (csmaxyLine.data()), + pCount, 1, GDT_Float64, + sizeof(double), sizeof(double)*pCount, NULL); + + if (status != 0) + { + std::cout << "Error read line " << lOff + ii << " from y-direction chip size max file: " << csmaxyname << "\n"; + GDALClose(csmaxyDS); + GDALDestroyDriverManager(); + exit(105); + } + } + + + + if (ssmname != "") + { + status = ssmDS->GetRasterBand(1)->RasterIO(GF_Read, + pOff, lOff+ii, + pCount, 1, + (void*) (ssmLine.data()), + pCount, 1, GDT_Float64, + sizeof(double), sizeof(double)*pCount, NULL); + + if (status != 0) + { + std::cout << "Error read line " << lOff + ii << " from stable surface mask file: " << ssmname << "\n"; + GDALClose(ssmDS); + GDALDestroyDriverManager(); + exit(105); + } + + } + + + + + + for (int jj=0; jjTransform(1, xyzs, xyzs+1, xyzs+2); + + for(int pp=0; pp<3; pp++) + { + targutm0[pp] = xyzs[pp]; + } + + xind = std::round((targutm0[0] - startingX) / XSize) + 1.; + yind = std::round((targutm0[1] - startingY) / YSize) + 1.; + + + + + + + //*********************Slant-range vector + + for(int pp=0; pp<3; pp++) + { + targutm[pp] = targutm0[pp]; + } + targutm[0] += XSize; + + + //Convert from LLH inplace to DEM coordinates + invTrans->Transform(1, targutm, targutm+1, targutm+2); + + for(int pp=0; pp<3; pp++) + { + xdiff[pp] = targutm[pp] - targxyz0[pp]; + } + unitvec_C(xdiff, xunit); + + + + + //*********************Along-track vector + + for(int pp=0; pp<3; pp++) + { + targutm[pp] = targutm0[pp]; + } + targutm[1] += YSize; + + + //Convert from LLH inplace to DEM coordinates + invTrans->Transform(1, targutm, targutm+1, targutm+2); + + for(int pp=0; pp<3; pp++) + { + ydiff[pp] = targutm[pp] - targxyz0[pp]; + } + unitvec_C(ydiff, yunit); + + + + + //*********************Local normal vector + if (dhdxname != "") + { + unitvec_C(slp, normal); + for(int pp=0; pp<3; pp++) + { + normal[pp] = -normal[pp]; + } + } + else + { + for(int pp=0; pp<3; pp++) + { + normal[pp] = 0.0; + } + } + + if (vxname != "") + { + vel[2] = -(vel[0]*normal[0]+vel[1]*normal[1])/normal[2]; + } + + if (srxname != "") + { + schrng1[2] = -(schrng1[0]*normal[0]+schrng1[1]*normal[1])/normal[2]; + schrng2[2] = -(schrng2[0]*normal[0]+schrng2[1]*normal[1])/normal[2]; + } + + + if ((xind > nPixels)|(xind < 1)|(yind > nLines)|(yind < 1)) + { + raster1[jj] = nodata_out; + raster2[jj] = nodata_out; + raster11[jj] = nodata_out; + raster22[jj] = nodata_out; + + sr_raster11[jj] = nodata_out; + sr_raster22[jj] = nodata_out; + csmin_raster11[jj] = nodata_out; + csmin_raster22[jj] = nodata_out; + csmax_raster11[jj] = nodata_out; + csmax_raster22[jj] = nodata_out; + ssm_raster[jj] = nodata_out; + + raster1a[jj] = nodata_out; + raster1b[jj] = nodata_out; +// raster1c[jj] = nodata_out; + raster2a[jj] = nodata_out; + raster2b[jj] = nodata_out; +// raster2c[jj] = nodata_out; + + } + else + { + raster1[jj] = xind; + raster2[jj] = yind; + + if (dhdxname != "") + { + + if (vxname != "") + { + if (vel[0] == nodata) + { + raster11[jj] = 0.; + raster22[jj] = 0.; + } + else + { + raster11[jj] = std::round(dot_C(vel,xunit)*dt/XSize/365.0/24.0/3600.0*1); + raster22[jj] = std::round(dot_C(vel,yunit)*dt/YSize/365.0/24.0/3600.0*1); + } + + } + + cross_C(xunit,yunit,cross); + unitvec_C(cross, cross); + cross_check = std::abs(std::acos(dot_C(normal,cross))/deg2rad-90.0); + + if (cross_check > 1.0) + { + raster1a[jj] = normal[2]/(dt/XSize/365.0/24.0/3600.0)*(normal[2]*yunit[1]-normal[1]*yunit[2])/((normal[2]*xunit[0]-normal[0]*xunit[2])*(normal[2]*yunit[1]-normal[1]*yunit[2])-(normal[2]*yunit[0]-normal[0]*yunit[2])*(normal[2]*xunit[1]-normal[1]*xunit[2])); + raster1b[jj] = -normal[2]/(dt/YSize/365.0/24.0/3600.0)*(normal[2]*xunit[1]-normal[1]*xunit[2])/((normal[2]*xunit[0]-normal[0]*xunit[2])*(normal[2]*yunit[1]-normal[1]*yunit[2])-(normal[2]*yunit[0]-normal[0]*yunit[2])*(normal[2]*xunit[1]-normal[1]*xunit[2])); + raster2a[jj] = -normal[2]/(dt/XSize/365.0/24.0/3600.0)*(normal[2]*yunit[0]-normal[0]*yunit[2])/((normal[2]*xunit[0]-normal[0]*xunit[2])*(normal[2]*yunit[1]-normal[1]*yunit[2])-(normal[2]*yunit[0]-normal[0]*yunit[2])*(normal[2]*xunit[1]-normal[1]*xunit[2])); + raster2b[jj] = normal[2]/(dt/YSize/365.0/24.0/3600.0)*(normal[2]*xunit[0]-normal[0]*xunit[2])/((normal[2]*xunit[0]-normal[0]*xunit[2])*(normal[2]*yunit[1]-normal[1]*yunit[2])-(normal[2]*yunit[0]-normal[0]*yunit[2])*(normal[2]*xunit[1]-normal[1]*xunit[2])); + } + else + { + raster1a[jj] = nodata_out; + raster1b[jj] = nodata_out; + raster2a[jj] = nodata_out; + raster2b[jj] = nodata_out; + } + + if (srxname != "") + { + if ((schrng1[0] == nodata)|(schrng1[0] == 0)) + { + sr_raster11[jj] = 0; + sr_raster22[jj] = 0; + } + else + { + sr_raster11[jj] = std::abs(std::round(dot_C(schrng1,xunit)*dt/XSize/365.0/24.0/3600.0*1)); + sr_raster22[jj] = std::abs(std::round(dot_C(schrng1,yunit)*dt/YSize/365.0/24.0/3600.0*1)); + if (std::abs(std::round(dot_C(schrng2,xunit)*dt/XSize/365.0/24.0/3600.0*1)) > sr_raster11[jj]) + { + sr_raster11[jj] = std::abs(std::round(dot_C(schrng2,xunit)*dt/XSize/365.0/24.0/3600.0*1)); + } + if (std::abs(std::round(dot_C(schrng2,yunit)*dt/YSize/365.0/24.0/3600.0*1)) > sr_raster22[jj]) + { + sr_raster22[jj] = std::abs(std::round(dot_C(schrng2,yunit)*dt/YSize/365.0/24.0/3600.0*1)); + } + if (sr_raster11[jj] == 0) + { + sr_raster11[jj] = 1; + } + if (sr_raster22[jj] == 0) + { + sr_raster22[jj] = 1; + } + } + } + + } + + + + if (csminxname != "") + { + if (csminxLine[jj] == nodata) + { + csmin_raster11[jj] = nodata_out; + csmin_raster22[jj] = nodata_out; + } + else + { + csmin_raster11[jj] = csminxLine[jj] / ChipSizeX0 * ChipSizeX0_PIX_X; + csmin_raster22[jj] = csminyLine[jj] / ChipSizeX0 * ChipSizeX0_PIX_Y; + } + } + + + if (csmaxxname != "") + { + if (csmaxxLine[jj] == nodata) + { + csmax_raster11[jj] = nodata_out; + csmax_raster22[jj] = nodata_out; + } + else + { + csmax_raster11[jj] = csmaxxLine[jj] / ChipSizeX0 * ChipSizeX0_PIX_X; + csmax_raster22[jj] = csmaxyLine[jj] / ChipSizeX0 * ChipSizeX0_PIX_Y; + } + } + + + if (ssmname != "") + { + if (ssmLine[jj] == nodata) + { + ssm_raster[jj] = nodata_out; + } + else + { + ssm_raster[jj] = ssmLine[jj]; + } + } + + + +// raster1a[jj] = los[0]*dt/dr/365.0/24.0/3600.0; +// raster1b[jj] = los[1]*dt/dr/365.0/24.0/3600.0; +// raster1c[jj] = los[2]*dt/dr/365.0/24.0/3600.0; +// raster2a[jj] = temp[0]*dt/norm_C(alt)/365.0/24.0/3600.0; +// raster2b[jj] = temp[1]*dt/norm_C(alt)/365.0/24.0/3600.0; +// raster2c[jj] = temp[2]*dt/norm_C(alt)/365.0/24.0/3600.0; + } + + +// std::cout << ii << " " << jj << "\n"; +// std::cout << rgind << " " << azind << "\n"; +// std::cout << raster1[jj][ii] << " " << raster2[jj][ii] << "\n"; +// std::cout << raster1[ii][jj] << "\n"; + } + + + + poBand1->RasterIO( GF_Write, 0, ii, pCount, 1, + raster1, pCount, 1, GDT_Int32, 0, 0 ); + poBand2->RasterIO( GF_Write, 0, ii, pCount, 1, + raster2, pCount, 1, GDT_Int32, 0, 0 ); + + if ((dhdxname != "")&(vxname != "")) + { + poBand1Off->RasterIO( GF_Write, 0, ii, pCount, 1, + raster11, pCount, 1, GDT_Int32, 0, 0 ); + poBand2Off->RasterIO( GF_Write, 0, ii, pCount, 1, + raster22, pCount, 1, GDT_Int32, 0, 0 ); + } + + if ((dhdxname != "")&(srxname != "")) + { + poBand1Sch->RasterIO( GF_Write, 0, ii, pCount, 1, + sr_raster11, pCount, 1, GDT_Int32, 0, 0 ); + poBand2Sch->RasterIO( GF_Write, 0, ii, pCount, 1, + sr_raster22, pCount, 1, GDT_Int32, 0, 0 ); + } + + if (csminxname != "") + { + poBand1Min->RasterIO( GF_Write, 0, ii, pCount, 1, + csmin_raster11, pCount, 1, GDT_Int32, 0, 0 ); + poBand2Min->RasterIO( GF_Write, 0, ii, pCount, 1, + csmin_raster22, pCount, 1, GDT_Int32, 0, 0 ); + } + + if (csmaxxname != "") + { + poBand1Max->RasterIO( GF_Write, 0, ii, pCount, 1, + csmax_raster11, pCount, 1, GDT_Int32, 0, 0 ); + poBand2Max->RasterIO( GF_Write, 0, ii, pCount, 1, + csmax_raster22, pCount, 1, GDT_Int32, 0, 0 ); + } + + if (ssmname != "") + { + poBand1Msk->RasterIO( GF_Write, 0, ii, pCount, 1, + ssm_raster, pCount, 1, GDT_Int32, 0, 0 ); + } + + if (dhdxname != "") + { + poBand1RO2VX->RasterIO( GF_Write, 0, ii, pCount, 1, + raster1a, pCount, 1, GDT_Float64, 0, 0 ); + poBand2RO2VX->RasterIO( GF_Write, 0, ii, pCount, 1, + raster1b, pCount, 1, GDT_Float64, 0, 0 ); + // poBand3Los->RasterIO( GF_Write, 0, ii, pCount, 1, + // raster1c, pCount, 1, GDT_Float64, 0, 0 ); + poBand1RO2VY->RasterIO( GF_Write, 0, ii, pCount, 1, + raster2a, pCount, 1, GDT_Float64, 0, 0 ); + poBand2RO2VY->RasterIO( GF_Write, 0, ii, pCount, 1, + raster2b, pCount, 1, GDT_Float64, 0, 0 ); + // poBand3Alt->RasterIO( GF_Write, 0, ii, pCount, 1, + // raster2c, pCount, 1, GDT_Float64, 0, 0 ); + } + + + } + + /* Once we're done, close properly the dataset */ + GDALClose( (GDALDatasetH) poDstDS ); + + if ((dhdxname != "")&(vxname != "")) + { + /* Once we're done, close properly the dataset */ + GDALClose( (GDALDatasetH) poDstDSOff ); + } + + if ((dhdxname != "")&(srxname != "")) + { + /* Once we're done, close properly the dataset */ + GDALClose( (GDALDatasetH) poDstDSSch ); + } + + if (csminxname != "") + { + /* Once we're done, close properly the dataset */ + GDALClose( (GDALDatasetH) poDstDSMin ); + } + + if (csmaxxname != "") + { + /* Once we're done, close properly the dataset */ + GDALClose( (GDALDatasetH) poDstDSMax ); + } + + if (ssmname != "") + { + /* Once we're done, close properly the dataset */ + GDALClose( (GDALDatasetH) poDstDSMsk ); + } + + if (dhdxname != "") + { + /* Once we're done, close properly the dataset */ + GDALClose( (GDALDatasetH) poDstDSRO2VX ); + + /* Once we're done, close properly the dataset */ + GDALClose( (GDALDatasetH) poDstDSRO2VY ); + } + + + GDALClose(demDS); + + if (dhdxname != "") + { + GDALClose(sxDS); + GDALClose(syDS); + } + + if (vxname != "") + { + GDALClose(vxDS); + GDALClose(vyDS); + } + + if (srxname != "") + { + GDALClose(srxDS); + GDALClose(sryDS); + } + + if (csminxname != "") + { + GDALClose(csminxDS); + GDALClose(csminyDS); + } + + if (csmaxxname != "") + { + GDALClose(csmaxxDS); + GDALClose(csmaxyDS); + } + + if (ssmname != "") + { + GDALClose(ssmDS); + } + + GDALDestroyDriverManager(); + +} + +void geoGridOptical::computeBbox(double *wesn) +{ + std::cout << "\nEstimated bounding box: \n" + << "West: " << wesn[0] << "\n" + << "East: " << wesn[1] << "\n" + << "South: " << wesn[2] << "\n" + << "North: " << wesn[3] << "\n"; +} + +double geoGridOptical::dot_C(double r_v[3], double r_w[3]) +{ + double dot; + dot = r_v[0]*r_w[0] + r_v[1]*r_w[1] + r_v[2]*r_w[2]; + return dot; +} + +void geoGridOptical::cross_C(double r_u[3], double r_v[3], double r_w[3]) +{ + r_w[0] = r_u[1]*r_v[2] - r_u[2]*r_v[1]; + r_w[1] = r_u[2]*r_v[0] - r_u[0]*r_v[2]; + r_w[2] = r_u[0]*r_v[1] - r_u[1]*r_v[0]; +} + +double geoGridOptical::norm_C(double r_v[3]) +{ + double norm; + norm = std::sqrt(r_v[0]*r_v[0] + r_v[1]*r_v[1] + r_v[2]*r_v[2]); + return norm; +} + + +void geoGridOptical::unitvec_C(double r_v[3], double r_w[3]) +{ + double norm; + norm = std::sqrt(r_v[0]*r_v[0] + r_v[1]*r_v[1] + r_v[2]*r_v[2]); + r_w[0] = r_v[0] / norm; + r_w[1] = r_v[1] / norm; + r_w[2] = r_v[2] / norm; +} diff --git a/setup.py b/setup.py old mode 100644 new mode 100755 index 6a57221..4a1a02c --- a/setup.py +++ b/setup.py @@ -30,18 +30,29 @@ sources= ['geo_autoRIFT/autoRIFT/bindings/autoriftcoremodule.cpp'], include_dirs=[np.get_include()] + ['geo_autoRIFT/autoRIFT/include', - os.path.join(path, 'include')], + os.path.join(path, 'include/opencv4/')], library_dirs = [os.path.join(path, 'lib')], libraries=['opencv_core', 'opencv_highgui', 'opencv_imgproc'], extra_compile_args=['-std=c++11'], language="c++" + ), + Extension( + name="geogrid/geogridOptical", + sources= ['geo_autoRIFT/geogrid/bindings/geogridOpticalmodule.cpp','geo_autoRIFT/geogrid/src/geogridOptical.cpp'], + include_dirs=[np.get_include()] + + ['geo_autoRIFT/geogrid/include', + os.path.join(path, 'include')], + library_dirs = [os.path.join(path, 'lib')], + libraries=['gomp','gdal'], + extra_compile_args=['-std=c++11'], + language="c++" ) ] -setup (name = 'autoRIFT', - version = '1.0.0', +setup (name = 'geo_autoRIFT', + version = '1.3.0', description = 'This is the autoRIFT python package', - package_dir={'autoRIFT': 'geo_autoRIFT/autoRIFT', 'geogrid': 'geo_autoRIFT/geogrid'}, - packages=['autoRIFT', 'geogrid'], -# scripts=['geo_autoRIFT/geogrid/GeogridOptical.py'], + package_dir={'autoRIFT': 'geo_autoRIFT/autoRIFT','geogrid': 'geo_autoRIFT/geogrid'}, + packages=['autoRIFT','geogrid'], +# scripts=['geo_autoRIFT/geogrid/GeogridOptical.py'], ext_modules = extensions) diff --git a/testGeogridOptical.py b/testGeogridOptical.py index 2cc7531..f54cfb4 100755 --- a/testGeogridOptical.py +++ b/testGeogridOptical.py @@ -84,6 +84,7 @@ def coregisterLoadMetadata(indir_m, indir_s): import re from geogrid import GeogridOptical +# import isce # from components.contrib.geo_autoRIFT.geogrid import GeogridOptical obj = GeogridOptical() @@ -132,6 +133,7 @@ def runGeogrid(info, info1, dem, dhdx, dhdy, vx, vy, srx, sry, csminx, csminy, c ''' from geogrid import GeogridOptical +# import isce # from components.contrib.geo_autoRIFT.geogrid import GeogridOptical from osgeo import gdal diff --git a/testGeogrid_ISCE.py b/testGeogrid_ISCE.py index e1b7eec..ac7b338 100755 --- a/testGeogrid_ISCE.py +++ b/testGeogrid_ISCE.py @@ -204,7 +204,7 @@ def runGeogrid(info, info1, dem, dhdx, dhdy, vx, vy, srx, sry, csminx, csminy, c import isce from components.contrib.geo_autoRIFT.geogrid import Geogrid -# from geogrid import Geogrid +# from geogrid import Geogrid from osgeo import gdal dem_info = gdal.Info(dem, format='json') diff --git a/testautoRIFT.py b/testautoRIFT.py index b7e940b..11e6a80 100755 --- a/testautoRIFT.py +++ b/testautoRIFT.py @@ -107,6 +107,7 @@ def loadProductOptical(file_m, file_s): Load the product using Product Manager. ''' from geogrid import GeogridOptical +# import isce # from components.contrib.geo_autoRIFT.geogrid import GeogridOptical obj = GeogridOptical() @@ -225,12 +226,12 @@ def runAutorift(I1, I2, xGrid, yGrid, Dx0, Dy0, SRx0, SRy0, CSMINx0, CSMINy0, CS obj.ChipSizeMinX = CSMINx0 if geogrid_run_info is None: - gridspacingx = float(str.split(subprocess.getoutput('fgrep "Grid spacing in m:" testGeogrid.txt'))[-1]) - chipsizex0 = float(str.split(subprocess.getoutput('fgrep "Smallest Allowable Chip Size in m:" testGeogrid.txt'))[-1]) + gridspacingx = float(str.split(runCmd('fgrep "Grid spacing in m:" testGeogrid.txt'))[-1]) + chipsizex0 = float(str.split(runCmd('fgrep "Smallest Allowable Chip Size in m:" testGeogrid.txt'))[-1]) try: - pixsizex = float(str.split(subprocess.getoutput('fgrep "Ground range pixel size:" testGeogrid.txt'))[-1]) + pixsizex = float(str.split(runCmd('fgrep "Ground range pixel size:" testGeogrid.txt'))[-1]) except: - pixsizex = float(str.split(subprocess.getoutput('fgrep "X-direction pixel size:" testGeogrid.txt'))[-1]) + pixsizex = float(str.split(runCmd('fgrep "X-direction pixel size:" testGeogrid.txt'))[-1]) else: gridspacingx = geogrid_run_info['gridspacingx'] chipsizex0 = geogrid_run_info['chipsizex0'] @@ -298,6 +299,7 @@ def runAutorift(I1, I2, xGrid, yGrid, Dx0, Dy0, SRx0, SRy0, CSMINx0, CSMINy0, CS # obj.sparseSearchSampleRate = 16 obj.OverSampleRatio = 64 +# obj.colfiltChunkSize = 4 # OverSampleRatio can be assigned as a scalar (such as the above line) or as a Python dictionary below for intellgient use (ChipSize-dependent). # Here, four chip sizes are used: ChipSize0X*[1,2,4,8] and four OverSampleRatio are considered [16,32,64,128]. The intelligent selection of OverSampleRatio (as a function of chip size) was determined by analyzing various combinations of (OverSampleRatio and chip size) and comparing the resulting image quality and statistics with the reference scenario (where the largest OverSampleRatio of 128 and chip size of ChipSize0X*8 are considered). @@ -398,7 +400,8 @@ def generateAutoriftProduct(indir_m, indir_s, grid_location, init_offset, search import numpy as np import time - # from components.contrib.geo_autoRIFT.autoRIFT import __version__ as version +# import isce +# from components.contrib.geo_autoRIFT.autoRIFT import __version__ as version from autoRIFT import __version__ as version if optical_flag == 1: @@ -719,6 +722,7 @@ def generateAutoriftProduct(indir_m, indir_s, grid_location, init_offset, search if nc_sensor == "S": if geogrid_run_info is None: chipsizex0 = float(str.split(runCmd('fgrep "Smallest Allowable Chip Size in m:" testGeogrid.txt'))[-1]) + gridspacingx = float(str.split(runCmd('fgrep "Grid spacing in m:" testGeogrid.txt'))[-1]) rangePixelSize = float(str.split(runCmd('fgrep "Ground range pixel size:" testGeogrid.txt'))[4]) azimuthPixelSize = float(str.split(runCmd('fgrep "Azimuth pixel size:" testGeogrid.txt'))[3]) dt = float(str.split(runCmd('fgrep "Repeat Time:" testGeogrid.txt'))[2]) @@ -726,6 +730,7 @@ def generateAutoriftProduct(indir_m, indir_s, grid_location, init_offset, search # print (str(rangePixelSize)+" "+str(azimuthPixelSize)) else: chipsizex0 = geogrid_run_info['chipsizex0'] + gridspacingx = geogrid_run_info['gridspacingx'] rangePixelSize = geogrid_run_info['XPixelSize'] azimuthPixelSize = geogrid_run_info['YPixelSize'] dt = geogrid_run_info['dt'] @@ -750,9 +755,9 @@ def generateAutoriftProduct(indir_m, indir_s, grid_location, init_offset, search PPP = roi_valid_percentage * 100 if ncname is None: out_nc_filename = f"./{master_filename[0:-4]}_X_{slave_filename[0:-4]}" \ - f"_G{chipsizex0:04.0f}V02_P{np.floor(PPP):03.0f}.nc" + f"_G{gridspacingx:04.0f}V02_P{np.floor(PPP):03.0f}.nc" else: - out_nc_filename = f"{ncname}_G{chipsizex0:04.0f}V02_P{np.floor(PPP):03.0f}.nc" + out_nc_filename = f"{ncname}_G{gridspacingx:04.0f}V02_P{np.floor(PPP):03.0f}.nc" CHIPSIZEY = np.round(CHIPSIZEX * ScaleChipSizeY / 2) * 2 @@ -771,7 +776,7 @@ def generateAutoriftProduct(indir_m, indir_s, grid_location, init_offset, search date_ct = d0 + (d1 - d0)/2 date_center = date_ct.strftime("%Y%m%d") - IMG_INFO_DICT = {'mission_img1':master_split[0][0],'sensor_img1':'C','satellite_img1':master_split[0][1:3],'acquisition_img1':master_dt,'absolute_orbit_number_img1':master_split[7],'mission_data_take_ID_img1':master_split[8],'product_unique_ID_img1':master_split[9][0:4],'mission_img2':slave_split[0][0],'sensor_img2':'C','satellite_img2':slave_split[0][1:3],'acquisition_img2':slave_dt,'absolute_orbit_number_img2':slave_split[7],'mission_data_take_ID_img2':slave_split[8],'product_unique_ID_img2':slave_split[9][0:4],'date_dt':date_dt,'date_center':date_center,'latitude':cen_lat,'longitude':cen_lon,'roi_valid_percentage':roi_valid_percentage,'autoRIFT_software_version':version} + IMG_INFO_DICT = {'mission_img1':master_split[0][0],'sensor_img1':'C','satellite_img1':master_split[0][1:3],'acquisition_img1':master_dt,'time_standard_img1':'UTC','absolute_orbit_number_img1':master_split[7],'mission_data_take_ID_img1':master_split[8],'product_unique_ID_img1':master_split[9][0:4],'mission_img2':slave_split[0][0],'sensor_img2':'C','satellite_img2':slave_split[0][1:3],'acquisition_img2':slave_dt,'time_standard_img2':'UTC','absolute_orbit_number_img2':slave_split[7],'mission_data_take_ID_img2':slave_split[8],'product_unique_ID_img2':slave_split[9][0:4],'date_dt':date_dt,'date_center':date_center,'latitude':cen_lat,'longitude':cen_lon,'roi_valid_percentage':PPP,'autoRIFT_software_version':version} error_vector = np.array([[0.0356, 0.0501, 0.0266, 0.0622, 0.0357, 0.0501], [0.5194, 1.1638, 0.3319, 1.3701, 0.5191, 1.1628]]) @@ -786,11 +791,13 @@ def generateAutoriftProduct(indir_m, indir_s, grid_location, init_offset, search elif nc_sensor == "L": if geogrid_run_info is None: chipsizex0 = float(str.split(runCmd('fgrep "Smallest Allowable Chip Size in m:" testGeogrid.txt'))[-1]) + gridspacingx = float(str.split(runCmd('fgrep "Grid spacing in m:" testGeogrid.txt'))[-1]) XPixelSize = float(str.split(runCmd('fgrep "X-direction pixel size:" testGeogrid.txt'))[3]) YPixelSize = float(str.split(runCmd('fgrep "Y-direction pixel size:" testGeogrid.txt'))[3]) epsg = float(str.split(runCmd('fgrep "EPSG:" testGeogrid.txt'))[1]) else: chipsizex0 = geogrid_run_info['chipsizex0'] + gridspacingx = geogrid_run_info['gridspacingx'] XPixelSize = geogrid_run_info['XPixelSize'] YPixelSize = geogrid_run_info['YPixelSize'] epsg = geogrid_run_info['epsg'] @@ -826,9 +833,9 @@ def generateAutoriftProduct(indir_m, indir_s, grid_location, init_offset, search PPP = roi_valid_percentage * 100 if ncname is None: out_nc_filename = f"./{master_filename[0:-8]}_X_{slave_filename[0:-8]}" \ - f"_G{chipsizex0:04.0f}V02_P{np.floor(PPP):03.0f}.nc" + f"_G{gridspacingx:04.0f}V02_P{np.floor(PPP):03.0f}.nc" else: - out_nc_filename = f"{ncname}_G{chipsizex0:04.0f}V02_P{np.floor(PPP):03.0f}.nc" + out_nc_filename = f"{ncname}_G{gridspacingx:04.0f}V02_P{np.floor(PPP):03.0f}.nc" CHIPSIZEY = np.round(CHIPSIZEX * ScaleChipSizeY / 2) * 2 @@ -849,7 +856,7 @@ def generateAutoriftProduct(indir_m, indir_s, grid_location, init_offset, search master_dt = master_split[3][0:8] + master_time.strftime("T%H:%M:%S") slave_dt = slave_split[3][0:8] + slave_time.strftime("T%H:%M:%S") - IMG_INFO_DICT = {'mission_img1':master_split[0][0],'sensor_img1':master_split[0][1],'satellite_img1':np.float64(master_split[0][2:4]),'correction_level_img1':master_split[1],'path_img1':np.float64(master_split[2][0:3]),'row_img1':np.float64(master_split[2][3:6]),'acquisition_date_img1':master_dt,'processing_date_img1':master_split[4][0:8],'collection_number_img1':np.float64(master_split[5]),'collection_category_img1':master_split[6],'mission_img2':slave_split[0][0],'sensor_img2':slave_split[0][1],'satellite_img2':np.float64(slave_split[0][2:4]),'correction_level_img2':slave_split[1],'path_img2':np.float64(slave_split[2][0:3]),'row_img2':np.float64(slave_split[2][3:6]),'acquisition_date_img2':slave_dt,'processing_date_img2':slave_split[4][0:8],'collection_number_img2':np.float64(slave_split[5]),'collection_category_img2':slave_split[6],'date_dt':date_dt,'date_center':date_center,'latitude':cen_lat,'longitude':cen_lon,'roi_valid_percentage':roi_valid_percentage,'autoRIFT_software_version':version} + IMG_INFO_DICT = {'mission_img1':master_split[0][0],'sensor_img1':master_split[0][1],'satellite_img1':np.float64(master_split[0][2:4]),'correction_level_img1':master_split[1],'path_img1':np.float64(master_split[2][0:3]),'row_img1':np.float64(master_split[2][3:6]),'acquisition_date_img1':master_dt,'time_standard_img1':'UTC','processing_date_img1':master_split[4][0:8],'collection_number_img1':np.float64(master_split[5]),'collection_category_img1':master_split[6],'mission_img2':slave_split[0][0],'sensor_img2':slave_split[0][1],'satellite_img2':np.float64(slave_split[0][2:4]),'correction_level_img2':slave_split[1],'path_img2':np.float64(slave_split[2][0:3]),'row_img2':np.float64(slave_split[2][3:6]),'acquisition_date_img2':slave_dt,'time_standard_img2':'UTC','processing_date_img2':slave_split[4][0:8],'collection_number_img2':np.float64(slave_split[5]),'collection_category_img2':slave_split[6],'date_dt':date_dt,'date_center':date_center,'latitude':cen_lat,'longitude':cen_lon,'roi_valid_percentage':PPP,'autoRIFT_software_version':version} error_vector = np.array([25.5,25.5]) @@ -864,11 +871,13 @@ def generateAutoriftProduct(indir_m, indir_s, grid_location, init_offset, search elif nc_sensor == "S2": if geogrid_run_info is None: chipsizex0 = float(str.split(runCmd('fgrep "Smallest Allowable Chip Size in m:" testGeogrid.txt'))[-1]) + gridspacingx = float(str.split(runCmd('fgrep "Grid spacing in m:" testGeogrid.txt'))[-1]) XPixelSize = float(str.split(runCmd('fgrep "X-direction pixel size:" testGeogrid.txt'))[3]) YPixelSize = float(str.split(runCmd('fgrep "Y-direction pixel size:" testGeogrid.txt'))[3]) epsg = float(str.split(runCmd('fgrep "EPSG:" testGeogrid.txt'))[1]) else: chipsizex0 = geogrid_run_info['chipsizex0'] + gridspacingx = geogrid_run_info['gridspacingx'] XPixelSize = geogrid_run_info['XPixelSize'] YPixelSize = geogrid_run_info['YPixelSize'] epsg = geogrid_run_info['epsg'] @@ -899,9 +908,9 @@ def generateAutoriftProduct(indir_m, indir_s, grid_location, init_offset, search PPP = roi_valid_percentage * 100 if ncname is None: out_nc_filename = f"./{master_filename[0:-8]}_X_{slave_filename[0:-8]}" \ - f"_G{chipsizex0:04.0f}V02_P{np.floor(PPP):03.0f}.nc" + f"_G{gridspacingx:04.0f}V02_P{np.floor(PPP):03.0f}.nc" else: - out_nc_filename = f"{ncname}_G{chipsizex0:04.0f}V02_P{np.floor(PPP):03.0f}.nc" + out_nc_filename = f"{ncname}_G{gridspacingx:04.0f}V02_P{np.floor(PPP):03.0f}.nc" CHIPSIZEY = np.round(CHIPSIZEX * ScaleChipSizeY / 2) * 2 from datetime import date @@ -921,7 +930,7 @@ def generateAutoriftProduct(indir_m, indir_s, grid_location, init_offset, search master_dt = master_split[2] + master_time.strftime("T%H:%M:%S") slave_dt = slave_split[2] + slave_time.strftime("T%H:%M:%S") - IMG_INFO_DICT = {'mission_img1':master_split[0][-3],'satellite_img1':master_split[0][-2:],'correction_level_img1':master_split[4][:3],'acquisition_date_img1':master_dt,'mission_img2':slave_split[0][-3],'satellite_img2':slave_split[0][-2:],'correction_level_img2':slave_split[4][:3],'acquisition_date_img2':slave_dt,'date_dt':date_dt,'date_center':date_center,'latitude':cen_lat,'longitude':cen_lon,'roi_valid_percentage':roi_valid_percentage,'autoRIFT_software_version':version} + IMG_INFO_DICT = {'mission_img1':master_split[0][-3],'satellite_img1':master_split[0][-2:],'correction_level_img1':master_split[4][:3],'acquisition_date_img1':master_dt,'time_standard_img1':'UTC','mission_img2':slave_split[0][-3],'satellite_img2':slave_split[0][-2:],'correction_level_img2':slave_split[4][:3],'acquisition_date_img2':slave_dt,'time_standard_img2':'UTC','date_dt':date_dt,'date_center':date_center,'latitude':cen_lat,'longitude':cen_lon,'roi_valid_percentage':PPP,'autoRIFT_software_version':version} error_vector = np.array([25.5,25.5]) diff --git a/testautoRIFT_ISCE.py b/testautoRIFT_ISCE.py index e02c7bc..bd6d353 100755 --- a/testautoRIFT_ISCE.py +++ b/testautoRIFT_ISCE.py @@ -225,12 +225,12 @@ def runAutorift(I1, I2, xGrid, yGrid, Dx0, Dy0, SRx0, SRy0, CSMINx0, CSMINy0, CS obj.ChipSizeMinX = CSMINx0 if geogrid_run_info is None: - gridspacingx = float(str.split(subprocess.getoutput('fgrep "Grid spacing in m:" testGeogrid.txt'))[-1]) - chipsizex0 = float(str.split(subprocess.getoutput('fgrep "Smallest Allowable Chip Size in m:" testGeogrid.txt'))[-1]) + gridspacingx = float(str.split(runCmd('fgrep "Grid spacing in m:" testGeogrid.txt'))[-1]) + chipsizex0 = float(str.split(runCmd('fgrep "Smallest Allowable Chip Size in m:" testGeogrid.txt'))[-1]) try: - pixsizex = float(str.split(subprocess.getoutput('fgrep "Ground range pixel size:" testGeogrid.txt'))[-1]) + pixsizex = float(str.split(runCmd('fgrep "Ground range pixel size:" testGeogrid.txt'))[-1]) except: - pixsizex = float(str.split(subprocess.getoutput('fgrep "X-direction pixel size:" testGeogrid.txt'))[-1]) + pixsizex = float(str.split(runCmd('fgrep "X-direction pixel size:" testGeogrid.txt'))[-1]) else: gridspacingx = geogrid_run_info['gridspacingx'] chipsizex0 = geogrid_run_info['chipsizex0'] @@ -298,6 +298,7 @@ def runAutorift(I1, I2, xGrid, yGrid, Dx0, Dy0, SRx0, SRy0, CSMINx0, CSMINy0, CS # obj.sparseSearchSampleRate = 16 obj.OverSampleRatio = 64 +# obj.colfiltChunkSize = 4 # OverSampleRatio can be assigned as a scalar (such as the above line) or as a Python dictionary below for intellgient use (ChipSize-dependent). # Here, four chip sizes are used: ChipSize0X*[1,2,4,8] and four OverSampleRatio are considered [16,32,64,128]. The intelligent selection of OverSampleRatio (as a function of chip size) was determined by analyzing various combinations of (OverSampleRatio and chip size) and comparing the resulting image quality and statistics with the reference scenario (where the largest OverSampleRatio of 128 and chip size of ChipSize0X*8 are considered). @@ -721,6 +722,7 @@ def generateAutoriftProduct(indir_m, indir_s, grid_location, init_offset, search if nc_sensor == "S": if geogrid_run_info is None: chipsizex0 = float(str.split(runCmd('fgrep "Smallest Allowable Chip Size in m:" testGeogrid.txt'))[-1]) + gridspacingx = float(str.split(runCmd('fgrep "Grid spacing in m:" testGeogrid.txt'))[-1]) rangePixelSize = float(str.split(runCmd('fgrep "Ground range pixel size:" testGeogrid.txt'))[4]) azimuthPixelSize = float(str.split(runCmd('fgrep "Azimuth pixel size:" testGeogrid.txt'))[3]) dt = float(str.split(runCmd('fgrep "Repeat Time:" testGeogrid.txt'))[2]) @@ -728,6 +730,7 @@ def generateAutoriftProduct(indir_m, indir_s, grid_location, init_offset, search # print (str(rangePixelSize)+" "+str(azimuthPixelSize)) else: chipsizex0 = geogrid_run_info['chipsizex0'] + gridspacingx = geogrid_run_info['gridspacingx'] rangePixelSize = geogrid_run_info['XPixelSize'] azimuthPixelSize = geogrid_run_info['YPixelSize'] dt = geogrid_run_info['dt'] @@ -752,9 +755,9 @@ def generateAutoriftProduct(indir_m, indir_s, grid_location, init_offset, search PPP = roi_valid_percentage * 100 if ncname is None: out_nc_filename = f"./{master_filename[0:-4]}_X_{slave_filename[0:-4]}" \ - f"_G{chipsizex0:04.0f}V02_P{np.floor(PPP):03.0f}.nc" + f"_G{gridspacingx:04.0f}V02_P{np.floor(PPP):03.0f}.nc" else: - out_nc_filename = f"{ncname}_G{chipsizex0:04.0f}V02_P{np.floor(PPP):03.0f}.nc" + out_nc_filename = f"{ncname}_G{gridspacingx:04.0f}V02_P{np.floor(PPP):03.0f}.nc" CHIPSIZEY = np.round(CHIPSIZEX * ScaleChipSizeY / 2) * 2 @@ -773,7 +776,7 @@ def generateAutoriftProduct(indir_m, indir_s, grid_location, init_offset, search date_ct = d0 + (d1 - d0)/2 date_center = date_ct.strftime("%Y%m%d") - IMG_INFO_DICT = {'mission_img1':master_split[0][0],'sensor_img1':'C','satellite_img1':master_split[0][1:3],'acquisition_img1':master_dt,'absolute_orbit_number_img1':master_split[7],'mission_data_take_ID_img1':master_split[8],'product_unique_ID_img1':master_split[9][0:4],'mission_img2':slave_split[0][0],'sensor_img2':'C','satellite_img2':slave_split[0][1:3],'acquisition_img2':slave_dt,'absolute_orbit_number_img2':slave_split[7],'mission_data_take_ID_img2':slave_split[8],'product_unique_ID_img2':slave_split[9][0:4],'date_dt':date_dt,'date_center':date_center,'latitude':cen_lat,'longitude':cen_lon,'roi_valid_percentage':roi_valid_percentage,'autoRIFT_software_version':version} + IMG_INFO_DICT = {'mission_img1':master_split[0][0],'sensor_img1':'C','satellite_img1':master_split[0][1:3],'acquisition_img1':master_dt,'time_standard_img1':'UTC','absolute_orbit_number_img1':master_split[7],'mission_data_take_ID_img1':master_split[8],'product_unique_ID_img1':master_split[9][0:4],'mission_img2':slave_split[0][0],'sensor_img2':'C','satellite_img2':slave_split[0][1:3],'acquisition_img2':slave_dt,'time_standard_img2':'UTC','absolute_orbit_number_img2':slave_split[7],'mission_data_take_ID_img2':slave_split[8],'product_unique_ID_img2':slave_split[9][0:4],'date_dt':date_dt,'date_center':date_center,'latitude':cen_lat,'longitude':cen_lon,'roi_valid_percentage':PPP,'autoRIFT_software_version':version} error_vector = np.array([[0.0356, 0.0501, 0.0266, 0.0622, 0.0357, 0.0501], [0.5194, 1.1638, 0.3319, 1.3701, 0.5191, 1.1628]]) @@ -788,11 +791,13 @@ def generateAutoriftProduct(indir_m, indir_s, grid_location, init_offset, search elif nc_sensor == "L": if geogrid_run_info is None: chipsizex0 = float(str.split(runCmd('fgrep "Smallest Allowable Chip Size in m:" testGeogrid.txt'))[-1]) + gridspacingx = float(str.split(runCmd('fgrep "Grid spacing in m:" testGeogrid.txt'))[-1]) XPixelSize = float(str.split(runCmd('fgrep "X-direction pixel size:" testGeogrid.txt'))[3]) YPixelSize = float(str.split(runCmd('fgrep "Y-direction pixel size:" testGeogrid.txt'))[3]) epsg = float(str.split(runCmd('fgrep "EPSG:" testGeogrid.txt'))[1]) else: chipsizex0 = geogrid_run_info['chipsizex0'] + gridspacingx = geogrid_run_info['gridspacingx'] XPixelSize = geogrid_run_info['XPixelSize'] YPixelSize = geogrid_run_info['YPixelSize'] epsg = geogrid_run_info['epsg'] @@ -828,9 +833,9 @@ def generateAutoriftProduct(indir_m, indir_s, grid_location, init_offset, search PPP = roi_valid_percentage * 100 if ncname is None: out_nc_filename = f"./{master_filename[0:-8]}_X_{slave_filename[0:-8]}" \ - f"_G{chipsizex0:04.0f}V02_P{np.floor(PPP):03.0f}.nc" + f"_G{gridspacingx:04.0f}V02_P{np.floor(PPP):03.0f}.nc" else: - out_nc_filename = f"{ncname}_G{chipsizex0:04.0f}V02_P{np.floor(PPP):03.0f}.nc" + out_nc_filename = f"{ncname}_G{gridspacingx:04.0f}V02_P{np.floor(PPP):03.0f}.nc" CHIPSIZEY = np.round(CHIPSIZEX * ScaleChipSizeY / 2) * 2 from datetime import date @@ -850,7 +855,7 @@ def generateAutoriftProduct(indir_m, indir_s, grid_location, init_offset, search master_dt = master_split[3][0:8] + master_time.strftime("T%H:%M:%S") slave_dt = slave_split[3][0:8] + slave_time.strftime("T%H:%M:%S") - IMG_INFO_DICT = {'mission_img1':master_split[0][0],'sensor_img1':master_split[0][1],'satellite_img1':np.float64(master_split[0][2:4]),'correction_level_img1':master_split[1],'path_img1':np.float64(master_split[2][0:3]),'row_img1':np.float64(master_split[2][3:6]),'acquisition_date_img1':master_dt,'processing_date_img1':master_split[4][0:8],'collection_number_img1':np.float64(master_split[5]),'collection_category_img1':master_split[6],'mission_img2':slave_split[0][0],'sensor_img2':slave_split[0][1],'satellite_img2':np.float64(slave_split[0][2:4]),'correction_level_img2':slave_split[1],'path_img2':np.float64(slave_split[2][0:3]),'row_img2':np.float64(slave_split[2][3:6]),'acquisition_date_img2':slave_dt,'processing_date_img2':slave_split[4][0:8],'collection_number_img2':np.float64(slave_split[5]),'collection_category_img2':slave_split[6],'date_dt':date_dt,'date_center':date_center,'latitude':cen_lat,'longitude':cen_lon,'roi_valid_percentage':roi_valid_percentage,'autoRIFT_software_version':version} + IMG_INFO_DICT = {'mission_img1':master_split[0][0],'sensor_img1':master_split[0][1],'satellite_img1':np.float64(master_split[0][2:4]),'correction_level_img1':master_split[1],'path_img1':np.float64(master_split[2][0:3]),'row_img1':np.float64(master_split[2][3:6]),'acquisition_date_img1':master_dt,'time_standard_img1':'UTC','processing_date_img1':master_split[4][0:8],'collection_number_img1':np.float64(master_split[5]),'collection_category_img1':master_split[6],'mission_img2':slave_split[0][0],'sensor_img2':slave_split[0][1],'satellite_img2':np.float64(slave_split[0][2:4]),'correction_level_img2':slave_split[1],'path_img2':np.float64(slave_split[2][0:3]),'row_img2':np.float64(slave_split[2][3:6]),'acquisition_date_img2':slave_dt,'time_standard_img2':'UTC','processing_date_img2':slave_split[4][0:8],'collection_number_img2':np.float64(slave_split[5]),'collection_category_img2':slave_split[6],'date_dt':date_dt,'date_center':date_center,'latitude':cen_lat,'longitude':cen_lon,'roi_valid_percentage':PPP,'autoRIFT_software_version':version} error_vector = np.array([25.5,25.5]) @@ -865,11 +870,13 @@ def generateAutoriftProduct(indir_m, indir_s, grid_location, init_offset, search elif nc_sensor == "S2": if geogrid_run_info is None: chipsizex0 = float(str.split(runCmd('fgrep "Smallest Allowable Chip Size in m:" testGeogrid.txt'))[-1]) + gridspacingx = float(str.split(runCmd('fgrep "Grid spacing in m:" testGeogrid.txt'))[-1]) XPixelSize = float(str.split(runCmd('fgrep "X-direction pixel size:" testGeogrid.txt'))[3]) YPixelSize = float(str.split(runCmd('fgrep "Y-direction pixel size:" testGeogrid.txt'))[3]) epsg = float(str.split(runCmd('fgrep "EPSG:" testGeogrid.txt'))[1]) else: chipsizex0 = geogrid_run_info['chipsizex0'] + gridspacingx = geogrid_run_info['gridspacingx'] XPixelSize = geogrid_run_info['XPixelSize'] YPixelSize = geogrid_run_info['YPixelSize'] epsg = geogrid_run_info['epsg'] @@ -900,9 +907,9 @@ def generateAutoriftProduct(indir_m, indir_s, grid_location, init_offset, search PPP = roi_valid_percentage * 100 if ncname is None: out_nc_filename = f"./{master_filename[0:-8]}_X_{slave_filename[0:-8]}" \ - f"_G{chipsizex0:04.0f}V02_P{np.floor(PPP):03.0f}.nc" + f"_G{gridspacingx:04.0f}V02_P{np.floor(PPP):03.0f}.nc" else: - out_nc_filename = f"{ncname}_G{chipsizex0:04.0f}V02_P{np.floor(PPP):03.0f}.nc" + out_nc_filename = f"{ncname}_G{gridspacingx:04.0f}V02_P{np.floor(PPP):03.0f}.nc" CHIPSIZEY = np.round(CHIPSIZEX * ScaleChipSizeY / 2) * 2 from datetime import date @@ -922,7 +929,7 @@ def generateAutoriftProduct(indir_m, indir_s, grid_location, init_offset, search master_dt = master_split[2] + master_time.strftime("T%H:%M:%S") slave_dt = slave_split[2] + slave_time.strftime("T%H:%M:%S") - IMG_INFO_DICT = {'mission_img1':master_split[0][-3],'satellite_img1':master_split[0][-2:],'correction_level_img1':master_split[4][:3],'acquisition_date_img1':master_dt,'mission_img2':slave_split[0][-3],'satellite_img2':slave_split[0][-2:],'correction_level_img2':slave_split[4][:3],'acquisition_date_img2':slave_dt,'date_dt':date_dt,'date_center':date_center,'latitude':cen_lat,'longitude':cen_lon,'roi_valid_percentage':roi_valid_percentage,'autoRIFT_software_version':version} + IMG_INFO_DICT = {'mission_img1':master_split[0][-3],'satellite_img1':master_split[0][-2:],'correction_level_img1':master_split[4][:3],'acquisition_date_img1':master_dt,'time_standard_img1':'UTC','mission_img2':slave_split[0][-3],'satellite_img2':slave_split[0][-2:],'correction_level_img2':slave_split[4][:3],'acquisition_date_img2':slave_dt,'time_standard_img2':'UTC','date_dt':date_dt,'date_center':date_center,'latitude':cen_lat,'longitude':cen_lon,'roi_valid_percentage':PPP,'autoRIFT_software_version':version} error_vector = np.array([25.5,25.5])