diff --git a/geo_autoRIFT/SConscript b/geo_autoRIFT/SConscript index e49f703..97d9cad 100755 --- a/geo_autoRIFT/SConscript +++ b/geo_autoRIFT/SConscript @@ -22,6 +22,7 @@ import os Import('envcontrib') package = 'geo_autoRIFT' envgeoAutorift = envcontrib.Clone() +envgeoAutorift.MergeFlags('-std=c++11') envgeoAutorift['PACKAGE'] = envcontrib['PACKAGE'] + '/' + package install = envcontrib['PRJ_SCONS_INSTALL'] + '/' + envgeoAutorift['PACKAGE'] listFiles = ['__init__.py'] diff --git a/geo_autoRIFT/autoRIFT/bindings/autoriftcoremodule.cpp b/geo_autoRIFT/autoRIFT/bindings/autoriftcoremodule.cpp index 94238d9..ccada72 100755 --- a/geo_autoRIFT/autoRIFT/bindings/autoriftcoremodule.cpp +++ b/geo_autoRIFT/autoRIFT/bindings/autoriftcoremodule.cpp @@ -41,6 +41,7 @@ #include "iostream" #include "numpy/arrayobject.h" #include "opencv2/imgproc/imgproc.hpp" +#include "opencv2/imgproc/types_c.h" #include "opencv2/highgui/highgui.hpp" #include "opencv2/core/core.hpp" diff --git a/geo_autoRIFT/geogrid/GeogridOptical.py b/geo_autoRIFT/geogrid/GeogridOptical.py index 2832532..1ce6f88 100755 --- a/geo_autoRIFT/geogrid/GeogridOptical.py +++ b/geo_autoRIFT/geogrid/GeogridOptical.py @@ -155,10 +155,10 @@ def geogrid(self): # For now print inputs that were obtained - print("Optical Image parameters: ") + 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)) + print("Dimensions: " + str(self.numberOfSamples) + " " + str(self.numberOfLines) + "\n") print("Map inputs: ") print("EPSG: " + str(self.epsgDem)) @@ -167,25 +167,45 @@ def geogrid(self): 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)) - print("DEM: " + str(self.demname)) - print("Velocities: " + str(self.vxname) + " " + str(self.vyname)) - print("Search Range: " + str(self.srxname) + " " + str(self.sryname)) - print("Chip Size Min: " + str(self.csminxname) + " " + str(self.csminyname)) - print("Chip Size Max: " + str(self.csmaxxname) + " " + str(self.csmaxyname)) - print("Stable Surface Mask: " + str(self.ssmname)) - print("Slopes: " + str(self.dhdxname) + " " + str(self.dhdyname)) - print("Output Nodata Value: " + str(self.nodata_out)) - - print("Outputs: ") + 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)) - print("Window offsets: " + str(self.winoffname)) - print("Window search range: " + str(self.winsrname)) - print("Window chip size min: " + str(self.wincsminname)) - print("Window chip size max: " + str(self.wincsmaxname)) - print("Window stable surface mask: " + str(self.winssmname)) - print("Window rdr_off2vel_x vector: " + str(self.winro2vxname)) - print("Window rdr_off2vel_y vector: " + str(self.winro2vyname)) - + + 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 .... ") @@ -277,6 +297,8 @@ def geogrid(self): 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)) + print("Dimensions of geogrid: " + str(pCount) + " x " + str(lCount)) projDem = osr.SpatialReference() @@ -296,12 +318,14 @@ def geogrid(self): if (self.vxname != ""): nodata = vxDS.GetRasterBand(1).GetNoDataValue() + else: + nodata = 0 nodata_out = self.nodata_out pszFormat = "GTiff" - adfGeoTransform = ( self._xlim[0], (self._xlim[1]-self._xlim[0])/pCount, 0, self._ylim[1], 0, (self._ylim[0]-self._ylim[1])/lCount ) + adfGeoTransform = ( geoTrans[0] + pOff * geoTrans[1], geoTrans[1], 0, geoTrans[3] + lOff * geoTrans[5], 0, geoTrans[5] ) oSRS = osr.SpatialReference() pszSRS_WKT = projDem.ExportToWkt() @@ -323,115 +347,113 @@ def geogrid(self): - - 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) - - - - 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) - - - 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) - - - 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.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) - 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 ) + 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) - poBand1Msk = poDstDSMsk.GetRasterBand(1) - poBand1Msk.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) - 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) + 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) + 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) @@ -512,6 +534,8 @@ def geogrid(self): 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]) schrng2 = np.array([-srxLine[jj], sryLine[jj], 0.0]) @@ -576,98 +600,106 @@ def geogrid(self): # else: # raster11[jj] = 0. # raster22[jj] = 0. - if (self.vxname == ""): -# pdb.set_trace() - raster11[jj] = 0. - raster22[jj] = 0. - elif (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) - - - if (self.srxname == ""): - sr_raster11[jj] = nodata_out - sr_raster22[jj] = nodata_out - elif (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.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 - else: - csmin_raster11[jj] = nodata_out - csmin_raster22[jj] = nodata_out + if (self.csmaxxname != ""): csmax_raster11[jj] = csmaxxLine[jj] / self.chipSizeX0 * ChipSizeX0_PIX_X csmax_raster22[jj] = csmaxyLine[jj] / self.chipSizeX0 * ChipSizeX0_PIX_Y - else: - csmax_raster11[jj] = nodata_out - csmax_raster22[jj] = nodata_out + if (self.ssmname != ""): ssm_raster[jj] = ssmLine[jj] - else: - ssm_raster[jj] = nodata_out - + - 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])); + # 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) - 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) - 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) - 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) - 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) - 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) - 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) + 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 - - poDstDSOff = None - - poDstDSSch = None - - poDstDSMin = None - - poDstDSMax = 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 - poDstDSMsk = None - - poDstDSRO2VX = None - - poDstDSRO2VY = None + poDstDSRO2VY = None demDS = None diff --git a/geo_autoRIFT/geogrid/bindings/geogridmodule.cpp b/geo_autoRIFT/geogrid/bindings/geogridmodule.cpp index f005f3e..4a5c3f3 100755 --- a/geo_autoRIFT/geogrid/bindings/geogridmodule.cpp +++ b/geo_autoRIFT/geogrid/bindings/geogridmodule.cpp @@ -403,6 +403,7 @@ PyObject* setRO2VYFilename(PyObject *self, PyObject *args) return Py_BuildValue("i", 0); } + PyObject* setLookSide(PyObject *self, PyObject *args) { uint64_t ptr; diff --git a/geo_autoRIFT/geogrid/src/geogrid.cpp b/geo_autoRIFT/geogrid/src/geogrid.cpp index 578d313..abb2d2f 100755 --- a/geo_autoRIFT/geogrid/src/geogrid.cpp +++ b/geo_autoRIFT/geogrid/src/geogrid.cpp @@ -60,24 +60,68 @@ void geoGrid::geogrid() std::cout << "XLimits: " << xmin << " " << xmax << "\n"; std::cout << "YLimits: " << ymin << " " << ymax << "\n"; std::cout << "Extent in km: " << (xmax - xmin)/1000. << " " << (ymax - ymin)/1000. << "\n"; - std::cout << "DEM: " << demname << "\n"; - std::cout << "Velocities: " << vxname << " " << vyname << "\n"; - std::cout << "Search Range: " << srxname << " " << sryname << "\n"; - std::cout << "Chip Size Min: " << csminxname << " " << csminyname << "\n"; - std::cout << "Chip Size Max: " << csmaxxname << " " << csmaxyname << "\n"; - std::cout << "Slopes: " << dhdxname << " " << dhdyname << "\n"; - std::cout << "Stable Surface Mask: " << ssmname << "\n"; - std::cout << "Output Nodata Value: " << nodata_out << "\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"; - std::cout << "Window offsets: " << offsetname << "\n"; - std::cout << "Window search range: " << searchrangename << "\n"; - std::cout << "Window chip size min: " << chipsizeminname << "\n"; - std::cout << "Window chip size max: " << chipsizemaxname << "\n"; - std::cout << "Window rdr_off2vel_x vector: " << ro2vx_name << "\n"; - std::cout << "Window rdr_off2vel_y vector: " << ro2vy_name << "\n"; - std::cout << "Window stable surface mask: " << stablesurfacemaskname << "\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"; @@ -223,12 +267,15 @@ void geoGrid::geogrid() exit(101); } } - if (ssmDS == NULL) + if (ssmname != "") { - std::cout << "Error opening stable surface mask file { " << ssmname << " }\n"; - std::cout << "Exiting with error code .... (101) \n"; - GDALDestroyDriverManager(); - exit(101); + 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); @@ -251,7 +298,9 @@ void geoGrid::geogrid() std::cout << "Ylimits : " << geoTrans[3] + (lOff + lCount) * geoTrans[5] << " " << geoTrans[3] + lOff * geoTrans[5] << "\n"; - std::cout << "Dimensions of geogrid: " << pCount << " x " << lCount << "\n\n"; + std::cout << "Origin index (in DEM) of geogrid: " << pOff << " " << lOff << "\n"; + + std::cout << "Dimensions of geogrid: " << pCount << " x " << lCount << "\n"; //Create GDAL Transformers @@ -329,6 +378,32 @@ void geoGrid::geogrid() 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; // double nodata_out; @@ -342,7 +417,7 @@ void geoGrid::geogrid() const char *pszFormat = "GTiff"; char **papszOptions = NULL; std::string str = ""; - double adfGeoTransform[6] = { xmin, (xmax-xmin)/pCount, 0, ymax, 0, (ymin-ymax)/lCount }; + 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 ); @@ -353,7 +428,7 @@ void geoGrid::geogrid() poDriver = GetGDALDriverManager()->GetDriverByName(pszFormat); if( poDriver == NULL ) exit(107); - GDALDataset *poDstDS; +// GDALDataset *poDstDS; str = pixlinename; const char * pszDstFilename = str.c_str(); @@ -366,188 +441,210 @@ void geoGrid::geogrid() // CPLFree( pszSRS_WKT ); - GDALRasterBand *poBand1; - GDALRasterBand *poBand2; +// 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); - - - - 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 ); + 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 *poBand1Sch; - GDALRasterBand *poBand2Sch; - poBand1Sch = poDstDSSch->GetRasterBand(1); - poBand2Sch = poDstDSSch->GetRasterBand(2); - poBand1Sch->SetNoDataValue(nodata_out); - poBand2Sch->SetNoDataValue(nodata_out); - - - - 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); +// 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); + } - GDALDriver *poDriverMax; - poDriverMax = GetGDALDriverManager()->GetDriverByName(pszFormat); - if( poDriverMax == NULL ) - exit(107); - GDALDataset *poDstDSMax; + 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); + + } - 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 ); + if (csmaxxname != "") + { - GDALRasterBand *poBand1Max; - GDALRasterBand *poBand2Max; - poBand1Max = poDstDSMax->GetRasterBand(1); - poBand2Max = poDstDSMax->GetRasterBand(2); - poBand1Max->SetNoDataValue(nodata_out); - poBand2Max->SetNoDataValue(nodata_out); + 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); + + } - 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 ); + if (dhdxname != "") + { - poDstDSMsk->SetGeoTransform( adfGeoTransform ); - poDstDSMsk->SetProjection( pszSRS_WKT ); + 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 *poBand1Msk; - poBand1Msk = poDstDSMsk->GetRasterBand(1); - poBand1Msk->SetNoDataValue(nodata_out); - - - - - 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); - - + +// 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 ); + 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); + + } - 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); // ground range and azimuth pixel size @@ -865,6 +962,8 @@ void geoGrid::geogrid() double los[3]; double alt[3]; double normal[3]; + double cross[3]; + double cross_check; double dopfact; double height; @@ -1110,92 +1209,104 @@ void geoGrid::geogrid() raster2a[jj] = nodata_out; raster2b[jj] = nodata_out; // raster2c[jj] = nodata_out; + } else { raster1[jj] = rgind; raster2[jj] = azind; - if (vxname == "") - { - raster11[jj] = 0.; - raster22[jj] = 0.; - } - else if (vel[0] == nodata) - { - raster11[jj] = 0.; - raster22[jj] = 0.; - } - else - { - raster11[jj] = std::round(dot_C(vel,los)*dt/dr/365.0/24.0/3600.0*1); - raster22[jj] = std::round(dot_C(vel,temp)*dt/norm_C(alt)/365.0/24.0/3600.0*1); - } - - if (srxname == "") + if (dhdxname != "") { - sr_raster11[jj] = nodata_out; - sr_raster22[jj] = nodata_out; - } - else if (vel[0] == nodata) - { - sr_raster11[jj] = 0; - sr_raster22[jj] = 0; - } - else - { - sr_raster11[jj] = std::abs(std::round(dot_C(schrng1,los)*dt/dr/365.0/24.0/3600.0*1)); - sr_raster22[jj] = std::abs(std::round(dot_C(schrng1,temp)*dt/norm_C(alt)/365.0/24.0/3600.0*1)); - if (std::abs(std::round(dot_C(schrng2,los)*dt/dr/365.0/24.0/3600.0*1)) > sr_raster11[jj]) + + if (vxname != "") { - sr_raster11[jj] = std::abs(std::round(dot_C(schrng2,los)*dt/dr/365.0/24.0/3600.0*1)); + if (vel[0] == nodata) + { + raster11[jj] = 0.; + raster22[jj] = 0.; + } + else + { + raster11[jj] = std::round(dot_C(vel,los)*dt/dr/365.0/24.0/3600.0*1); + raster22[jj] = std::round(dot_C(vel,temp)*dt/norm_C(alt)/365.0/24.0/3600.0*1); + } + } - if (std::abs(std::round(dot_C(schrng2,temp)*dt/norm_C(alt)/365.0/24.0/3600.0*1)) > sr_raster22[jj]) + + cross_C(los,temp,cross); + unitvec_C(cross, cross); + cross_check = std::abs(std::acos(dot_C(normal,cross))/deg2rad-90.0); + + if (cross_check > 1.0) { - sr_raster22[jj] = std::abs(std::round(dot_C(schrng2,temp)*dt/norm_C(alt)/365.0/24.0/3600.0*1)); + raster1a[jj] = normal[2]/(dt/dr/365.0/24.0/3600.0)*(normal[2]*temp[1]-normal[1]*temp[2])/((normal[2]*los[0]-normal[0]*los[2])*(normal[2]*temp[1]-normal[1]*temp[2])-(normal[2]*temp[0]-normal[0]*temp[2])*(normal[2]*los[1]-normal[1]*los[2])); + raster1b[jj] = -normal[2]/(dt/norm_C(alt)/365.0/24.0/3600.0)*(normal[2]*los[1]-normal[1]*los[2])/((normal[2]*los[0]-normal[0]*los[2])*(normal[2]*temp[1]-normal[1]*temp[2])-(normal[2]*temp[0]-normal[0]*temp[2])*(normal[2]*los[1]-normal[1]*los[2])); + raster2a[jj] = -normal[2]/(dt/dr/365.0/24.0/3600.0)*(normal[2]*temp[0]-normal[0]*temp[2])/((normal[2]*los[0]-normal[0]*los[2])*(normal[2]*temp[1]-normal[1]*temp[2])-(normal[2]*temp[0]-normal[0]*temp[2])*(normal[2]*los[1]-normal[1]*los[2])); + raster2b[jj] = normal[2]/(dt/norm_C(alt)/365.0/24.0/3600.0)*(normal[2]*los[0]-normal[0]*los[2])/((normal[2]*los[0]-normal[0]*los[2])*(normal[2]*temp[1]-normal[1]*temp[2])-(normal[2]*temp[0]-normal[0]*temp[2])*(normal[2]*los[1]-normal[1]*los[2])); } - if (sr_raster11[jj] == 0) + else { - sr_raster11[jj] = 1; + raster1a[jj] = nodata_out; + raster1b[jj] = nodata_out; + raster2a[jj] = nodata_out; + raster2b[jj] = nodata_out; } - if (sr_raster22[jj] == 0) + + if (srxname != "") { - sr_raster22[jj] = 1; + if ((vxname != "")&(vel[0] == nodata)) + { + sr_raster11[jj] = 0; + sr_raster22[jj] = 0; + } + else + { + sr_raster11[jj] = std::abs(std::round(dot_C(schrng1,los)*dt/dr/365.0/24.0/3600.0*1)); + sr_raster22[jj] = std::abs(std::round(dot_C(schrng1,temp)*dt/norm_C(alt)/365.0/24.0/3600.0*1)); + if (std::abs(std::round(dot_C(schrng2,los)*dt/dr/365.0/24.0/3600.0*1)) > sr_raster11[jj]) + { + sr_raster11[jj] = std::abs(std::round(dot_C(schrng2,los)*dt/dr/365.0/24.0/3600.0*1)); + } + if (std::abs(std::round(dot_C(schrng2,temp)*dt/norm_C(alt)/365.0/24.0/3600.0*1)) > sr_raster22[jj]) + { + sr_raster22[jj] = std::abs(std::round(dot_C(schrng2,temp)*dt/norm_C(alt)/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 != "") { csmin_raster11[jj] = csminxLine[jj] / ChipSizeX0 * ChipSizeX0_PIX_grd; csmin_raster22[jj] = csminyLine[jj] / ChipSizeX0 * ChipSizeX0_PIX_azm; } - else - { - csmin_raster11[jj] = nodata_out; - csmin_raster22[jj] = nodata_out; - } + + if (csmaxxname != "") { csmax_raster11[jj] = csmaxxLine[jj] / ChipSizeX0 * ChipSizeX0_PIX_grd; csmax_raster22[jj] = csmaxyLine[jj] / ChipSizeX0 * ChipSizeX0_PIX_azm; } - else - { - csmax_raster11[jj] = nodata_out; - csmax_raster22[jj] = nodata_out; - } + + if (ssmname != "") { ssm_raster[jj] = ssmLine[jj]; } - else - { - ssm_raster[jj] = nodata_out; - } - raster1a[jj] = normal[2]/(dt/dr/365.0/24.0/3600.0)*(normal[2]*temp[1]-normal[1]*temp[2])/((normal[2]*los[0]-normal[0]*los[2])*(normal[2]*temp[1]-normal[1]*temp[2])-(normal[2]*temp[0]-normal[0]*temp[2])*(normal[2]*los[1]-normal[1]*los[2])); - raster1b[jj] = -normal[2]/(dt/norm_C(alt)/365.0/24.0/3600.0)*(normal[2]*los[1]-normal[1]*los[2])/((normal[2]*los[0]-normal[0]*los[2])*(normal[2]*temp[1]-normal[1]*temp[2])-(normal[2]*temp[0]-normal[0]*temp[2])*(normal[2]*los[1]-normal[1]*los[2])); - raster2a[jj] = -normal[2]/(dt/dr/365.0/24.0/3600.0)*(normal[2]*temp[0]-normal[0]*temp[2])/((normal[2]*los[0]-normal[0]*los[2])*(normal[2]*temp[1]-normal[1]*temp[2])-(normal[2]*temp[0]-normal[0]*temp[2])*(normal[2]*los[1]-normal[1]*los[2])); - raster2b[jj] = normal[2]/(dt/norm_C(alt)/365.0/24.0/3600.0)*(normal[2]*los[0]-normal[0]*los[2])/((normal[2]*los[0]-normal[0]*los[2])*(normal[2]*temp[1]-normal[1]*temp[2])-(normal[2]*temp[0]-normal[0]*temp[2])*(normal[2]*los[1]-normal[1]*los[2])); + + // raster1a[jj] = los[0]*dt/dr/365.0/24.0/3600.0; // raster1b[jj] = los[1]*dt/dr/365.0/24.0/3600.0; @@ -1212,68 +1323,112 @@ void geoGrid::geogrid() // 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 ); - 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 ); - 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 ); - 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 ); - 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 ); - poBand1Msk->RasterIO( GF_Write, 0, ii, pCount, 1, - ssm_raster, 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 ); + } - 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 ); - /* Once we're done, close properly the dataset */ - GDALClose( (GDALDatasetH) poDstDSOff ); + if ((dhdxname != "")&(vxname != "")) + { + /* Once we're done, close properly the dataset */ + GDALClose( (GDALDatasetH) poDstDSOff ); + } - /* Once we're done, close properly the dataset */ - GDALClose( (GDALDatasetH) poDstDSSch ); + if ((dhdxname != "")&(srxname != "")) + { + /* Once we're done, close properly the dataset */ + GDALClose( (GDALDatasetH) poDstDSSch ); + } - /* Once we're done, close properly the dataset */ - GDALClose( (GDALDatasetH) poDstDSMin ); + if (csminxname != "") + { + /* Once we're done, close properly the dataset */ + GDALClose( (GDALDatasetH) poDstDSMin ); + } - /* Once we're done, close properly the dataset */ - GDALClose( (GDALDatasetH) poDstDSMax ); + if (csmaxxname != "") + { + /* Once we're done, close properly the dataset */ + GDALClose( (GDALDatasetH) poDstDSMax ); + } - /* Once we're done, close properly the dataset */ - GDALClose( (GDALDatasetH) poDstDSMsk ); + if (ssmname != "") + { + /* Once we're done, close properly the dataset */ + GDALClose( (GDALDatasetH) poDstDSMsk ); + } - /* Once we're done, close properly the dataset */ - GDALClose( (GDALDatasetH) poDstDSRO2VX ); + if (dhdxname != "") + { + /* Once we're done, close properly the dataset */ + GDALClose( (GDALDatasetH) poDstDSRO2VX ); + + /* Once we're done, close properly the dataset */ + GDALClose( (GDALDatasetH) poDstDSRO2VY ); + } - /* Once we're done, close properly the dataset */ - GDALClose( (GDALDatasetH) poDstDSRO2VY ); GDALClose(demDS); diff --git a/testautoRIFT.py b/testautoRIFT.py index e3e2a33..273d699 100755 --- a/testautoRIFT.py +++ b/testautoRIFT.py @@ -69,6 +69,8 @@ def cmdLineParse(): help='Input pixel offsets to vx conversion coefficients file name') parser.add_argument('-vy', '--input_vy', dest='offset2vy', type=str, required=False, help='Input pixel offsets to vy conversion coefficients file name') + parser.add_argument('-ssm', '--input_ssm', dest='stable_surface_mask', type=str, required=False, + help='Input stable surface mask file name') parser.add_argument('-fo', '--flag_optical', dest='optical_flag', type=bool, required=False, default=0, help='flag for reading optical data (e.g. Landsat): use 1 for on and 0 (default) for off') parser.add_argument('-nc', '--sensor_flag_netCDF', dest='nc_sensor', type=str, required=False, default=None, @@ -127,6 +129,7 @@ def runAutorift(I1, I2, xGrid, yGrid, Dx0, Dy0, SRx0, SRy0, CSMINx0, CSMINy0, CS import numpy as np # import isceobj import time + import subprocess obj = autoRIFT() @@ -144,6 +147,20 @@ def runAutorift(I1, I2, xGrid, yGrid, Dx0, Dy0, SRx0, SRy0, CSMINx0, CSMINy0, CS obj.I1 = I1 obj.I2 = I2 + + # test with lena image (533 X 533) +# obj.ChipSizeMinX=16 +# obj.ChipSizeMaxX=32 +# obj.ChipSize0X=16 +# obj.SkipSampleX=16 +# obj.SkipSampleY=16 + + # test with Venus image (407 X 407) +# obj.ChipSizeMinX=8 +# obj.ChipSizeMaxX=16 +# obj.ChipSize0X=8 +# obj.SkipSampleX=8 +# obj.SkipSampleY=8 # create the grid if it does not exist @@ -191,7 +208,6 @@ def runAutorift(I1, I2, xGrid, yGrid, Dx0, Dy0, SRx0, SRy0, CSMINx0, CSMINy0, CS if CSMINx0 is not None: obj.ChipSizeMaxX = CSMAXx0 obj.ChipSizeMinX = CSMINx0 - obj.ChipSize0X = np.min(CSMINx0[CSMINx0!=nodata]) chipsizex0 = float(str.split(subprocess.getoutput('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]) @@ -201,8 +217,11 @@ def runAutorift(I1, I2, xGrid, yGrid, Dx0, Dy0, SRx0, SRy0, CSMINx0, CSMINy0, CS # obj.ChipSize0X = np.min(CSMINx0[CSMINx0!=nodata]) RATIO_Y2X = CSMINy0/CSMINx0 obj.ScaleChipSizeY = np.mean(RATIO_Y2X[(CSMINx0!=nodata)&(CSMINy0!=nodata)]) +# obj.ChipSizeMaxX = obj.ChipSizeMaxX / obj.ChipSizeMaxX * 544 +# obj.ChipSizeMinX = obj.ChipSizeMinX / obj.ChipSizeMinX * 68 else: if ((optflag == 1)&(xGrid is not None)): +# print('test') obj.ChipSizeMaxX = 32 obj.ChipSizeMinX = 16 obj.ChipSize0X = 16 @@ -254,7 +273,7 @@ def runAutorift(I1, I2, xGrid, yGrid, Dx0, Dy0, SRx0, SRy0, CSMINx0, CSMINy0, CS # obj.sparseSearchSampleRate = 16 -# obj.OverSampleRatio = 32 + obj.OverSampleRatio = 64 # 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). @@ -352,6 +371,11 @@ def runAutorift(I1, I2, xGrid, yGrid, Dx0, Dy0, SRx0, SRy0, CSMINx0, CSMINy0, CS if inps.optical_flag == 1: data_m = loadProductOptical(inps.indir_m) data_s = loadProductOptical(inps.indir_s) + # test with lena/Venus image +# import scipy.io as sio +# conts = sio.loadmat(inps.indir_m) +# data_m = conts['I'] +# data_s = conts['I1'] else: data_m = loadProduct(inps.indir_m) data_s = loadProduct(inps.indir_s) @@ -421,6 +445,14 @@ def runAutorift(I1, I2, xGrid, yGrid, Dx0, Dy0, SRx0, SRy0, CSMINx0, CSMINy0, CS band=None ds=None + if inps.stable_surface_mask is not None: + ds = gdal.Open(inps.stable_surface_mask) + band = ds.GetRasterBand(1) + SSM = band.ReadAsArray() + SSM = SSM.astype('bool') + band=None + ds=None + Dx, Dy, InterpMask, ChipSizeX, ScaleChipSizeY, SearchLimitX, SearchLimitY, origSize, noDataMask = runAutorift(data_m, data_s, xGrid, yGrid, Dx0, Dy0, SRx0, SRy0, CSMINx0, CSMINy0, CSMAXx0, CSMAXy0, noDataMask, inps.optical_flag, nodata) @@ -447,6 +479,8 @@ def runAutorift(I1, I2, xGrid, yGrid, Dx0, Dy0, SRx0, SRy0, CSMINx0, CSMINy0, CS CHIPSIZEX[noDataMask] = 0 SEARCHLIMITX[noDataMask] = 0 SEARCHLIMITY[noDataMask] = 0 + if SSM is not None: + SSM[noDataMask] = False import scipy.io as sio sio.savemat('offset.mat',{'Dx':DX,'Dy':DY,'InterpMask':INTERPMASK,'ChipSizeX':CHIPSIZEX}) @@ -491,6 +525,7 @@ def runAutorift(I1, I2, xGrid, yGrid, Dx0, Dy0, SRx0, SRy0, CSMINx0, CSMINy0, CS if inps.offset2vx is not None: + ds = gdal.Open(inps.offset2vx) band = ds.GetRasterBand(1) offset2vx_1 = band.ReadAsArray() @@ -511,6 +546,8 @@ def runAutorift(I1, I2, xGrid, yGrid, Dx0, Dy0, SRx0, SRy0, CSMINx0, CSMINy0, CS VY = offset2vy_1 * DX + offset2vy_2 * DY VX = VX.astype(np.float32) VY = VY.astype(np.float32) + + ############ write velocity output in Geotiff format outRaster = driver.Create("velocity.tif", int(xGrid.shape[1]), int(xGrid.shape[0]), 2, gdal.GDT_Float32) outRaster.SetGeoTransform(tran) @@ -522,97 +559,193 @@ def runAutorift(I1, I2, xGrid, yGrid, Dx0, Dy0, SRx0, SRy0, CSMINx0, CSMINy0, CS outband.WriteArray(VY) outband.FlushCache() - ######################################################################################## - ############ netCDF packaging for Sentinel and Landsat dataset; can add other sensor format as well - if inps.nc_sensor == "S": + ############ prepare for netCDF packaging + + if inps.nc_sensor is not None: - 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]) - epsg = float(str.split(runCmd('fgrep "EPSG:" testGeogrid.txt'))[1]) -# print (str(rangePixelSize)+" "+str(azimuthPixelSize)) - - - runCmd('topsinsar_filename.py') -# import scipy.io as sio - conts = sio.loadmat('topsinsar_filename.mat') - master_filename = conts['master_filename'][0] - slave_filename = conts['slave_filename'][0] - master_split = str.split(master_filename,'_') - slave_split = str.split(slave_filename,'_') - - import netcdf_output as no - version = '1.0.5' - pair_type = 'radar' - detection_method = 'feature' - coordinates = 'radar' -# out_nc_filename = 'Jakobshavn.nc' - out_nc_filename = master_filename[0:-4]+'_'+slave_filename[0:-4]+'.nc' - out_nc_filename = './' + out_nc_filename - roi_valid_percentage = np.sum(CHIPSIZEX!=0)/np.sum(SEARCHLIMITX!=0)*100.0 - CHIPSIZEY = np.round(CHIPSIZEX * ScaleChipSizeY / 2) * 2 + vxrefname = str.split(runCmd('fgrep "Velocities:" testGeogrid.txt'))[1] + vyrefname = str.split(runCmd('fgrep "Velocities:" testGeogrid.txt'))[2] + sxname = str.split(runCmd('fgrep "Slopes:" testGeogrid.txt'))[1][:-4]+"s.tif" + syname = str.split(runCmd('fgrep "Slopes:" testGeogrid.txt'))[2][:-4]+"s.tif" + maskname = str.split(runCmd('fgrep "Slopes:" testGeogrid.txt'))[2][:-8]+"masks.tif" + xoff = int(str.split(runCmd('fgrep "Origin index (in DEM) of geogrid:" testGeogrid.txt'))[6]) + yoff = int(str.split(runCmd('fgrep "Origin index (in DEM) of geogrid:" testGeogrid.txt'))[7]) + xcount = int(str.split(runCmd('fgrep "Dimensions of geogrid:" testGeogrid.txt'))[3]) + ycount = int(str.split(runCmd('fgrep "Dimensions of geogrid:" testGeogrid.txt'))[5]) + + ds = gdal.Open(vxrefname) + band = ds.GetRasterBand(1) + VXref = band.ReadAsArray(xoff, yoff, xcount, ycount) + ds = None + band = None - + ds = gdal.Open(vyrefname) + band = ds.GetRasterBand(1) + VYref = band.ReadAsArray(xoff, yoff, xcount, ycount) + ds = None + band = None - from datetime import date - d0 = date(np.int(master_split[5][0:4]),np.int(master_split[5][4:6]),np.int(master_split[5][6:8])) - d1 = date(np.int(slave_split[5][0:4]),np.int(slave_split[5][4:6]),np.int(slave_split[5][6:8])) - date_dt_base = d1 - d0 - date_dt = np.float64(np.abs(date_dt_base.days)) - if date_dt_base.days < 0: - date_ct = d1 + (d0 - d1)/2 - date_center = date_ct.strftime("%Y%m%d") - else: - 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_split[5][0:8],'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_split[5][0:8],'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,'roi_valid_percentage':roi_valid_percentage,'autoRIFT_software_version':version} - - no.netCDF_packaging(VX, VY, DX, DY, INTERPMASK, CHIPSIZEX, CHIPSIZEY, rangePixelSize, azimuthPixelSize, dt, epsg, srs, tran, out_nc_filename, pair_type, detection_method, coordinates, IMG_INFO_DICT) - - elif inps.nc_sensor == "L": + ds = gdal.Open(sxname) + band = ds.GetRasterBand(1) + SX = band.ReadAsArray(xoff, yoff, xcount, ycount) + ds = None + band = None - 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]) + ds = gdal.Open(syname) + band = ds.GetRasterBand(1) + SY = band.ReadAsArray(xoff, yoff, xcount, ycount) + ds = None + band = None - master_filename = inps.indir_m - slave_filename = inps.indir_s - master_split = str.split(master_filename,'_') - slave_split = str.split(slave_filename,'_') + ds = gdal.Open(maskname) + band = ds.GetRasterBand(1) + MM = band.ReadAsArray(xoff, yoff, xcount, ycount) + ds = None + band = None - import netcdf_output as no - version = '1.0.5' - pair_type = 'optical' - detection_method = 'feature' - coordinates = 'map' -# out_nc_filename = 'Jakobshavn_opt.nc' - out_nc_filename = master_filename[0:-4]+'_'+slave_filename[0:-4]+'.nc' - out_nc_filename = './' + out_nc_filename - roi_valid_percentage = np.sum(CHIPSIZEX!=0)/np.sum(SEARCHLIMITX!=0)*100.0 - CHIPSIZEY = np.round(CHIPSIZEX * ScaleChipSizeY / 2) * 2 - - from datetime import date - d0 = date(np.int(master_split[3][0:4]),np.int(master_split[3][4:6]),np.int(master_split[3][6:8])) - d1 = date(np.int(slave_split[3][0:4]),np.int(slave_split[3][4:6]),np.int(slave_split[3][6:8])) - date_dt_base = d1 - d0 - date_dt = np.float64(np.abs(date_dt_base.days)) - if date_dt_base.days < 0: - date_ct = d1 + (d0 - d1)/2 - date_center = date_ct.strftime("%Y%m%d") - else: - date_ct = d0 + (d1 - d0)/2 - date_center = date_ct.strftime("%Y%m%d") + DXref = offset2vy_2 / (offset2vx_1 * offset2vy_2 - offset2vx_2 * offset2vy_1) * VXref - offset2vx_2 / (offset2vx_1 * offset2vy_2 - offset2vx_2 * offset2vy_1) * VYref + DYref = offset2vx_1 / (offset2vx_1 * offset2vy_2 - offset2vx_2 * offset2vy_1) * VYref - offset2vy_1 / (offset2vx_1 * offset2vy_2 - offset2vx_2 * offset2vy_1) * VXref + + stable_count = np.sum(SSM & np.logical_not(np.isnan(DX))) - 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_split[3][0:8],'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_split[3][0:8],'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,'roi_valid_percentage':roi_valid_percentage,'autoRIFT_software_version':version} + if stable_count == 0: + stable_shift_applied = 0 + else: + stable_shift_applied = 1 - no.netCDF_packaging(VX, VY, DX, DY, INTERPMASK, CHIPSIZEX, CHIPSIZEY, XPixelSize, YPixelSize, None, epsg, srs, tran, out_nc_filename, pair_type, detection_method, coordinates, IMG_INFO_DICT) - - elif inps.nc_sensor is None: - print('netCDF packaging not performed') + if stable_shift_applied == 1: + temp = DX.copy() - DXref.copy() + temp[np.logical_not(SSM)] = np.nan + dx_mean_shift = np.median(temp[(temp > -5)&(temp < 5)]) + DX = DX - dx_mean_shift + + temp = DY.copy() - DYref.copy() + temp[np.logical_not(SSM)] = np.nan + dy_mean_shift = np.median(temp[(temp > -5)&(temp < 5)]) + DY = DY - dy_mean_shift + else: + dx_mean_shift = 0.0 + dy_mean_shift = 0.0 + + VX = offset2vx_1 * DX + offset2vx_2 * DY + VY = offset2vy_1 * DX + offset2vy_2 * DY + VX = VX.astype(np.float32) + VY = VY.astype(np.float32) + ######################################################################################## + ############ netCDF packaging for Sentinel and Landsat dataset; can add other sensor format as well + if inps.nc_sensor == "S": + + 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]) + epsg = float(str.split(runCmd('fgrep "EPSG:" testGeogrid.txt'))[1]) + # print (str(rangePixelSize)+" "+str(azimuthPixelSize)) + + + runCmd('topsinsar_filename.py') + # import scipy.io as sio + conts = sio.loadmat('topsinsar_filename.mat') + master_filename = conts['master_filename'][0] + slave_filename = conts['slave_filename'][0] + master_dt = conts['master_dt'][0] + slave_dt = conts['slave_dt'][0] + master_split = str.split(master_filename,'_') + slave_split = str.split(slave_filename,'_') + + import netcdf_output as no + version = '1.0.7' + pair_type = 'radar' + detection_method = 'feature' + coordinates = 'radar' + # out_nc_filename = 'Jakobshavn.nc' + out_nc_filename = master_filename[0:-4]+'_'+slave_filename[0:-4]+'.nc' + out_nc_filename = './' + out_nc_filename + roi_valid_percentage = int(round(np.sum(CHIPSIZEX!=0)/np.sum(SEARCHLIMITX!=0)*1000.0))/1000 + CHIPSIZEY = np.round(CHIPSIZEX * ScaleChipSizeY / 2) * 2 + + + + from datetime import date + d0 = date(np.int(master_split[5][0:4]),np.int(master_split[5][4:6]),np.int(master_split[5][6:8])) + d1 = date(np.int(slave_split[5][0:4]),np.int(slave_split[5][4:6]),np.int(slave_split[5][6:8])) + date_dt_base = d1 - d0 + date_dt = np.float64(np.abs(date_dt_base.days)) + if date_dt_base.days < 0: + date_ct = d1 + (d0 - d1)/2 + date_center = date_ct.strftime("%Y%m%d") + else: + 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,'roi_valid_percentage':roi_valid_percentage,'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]]) + + no.netCDF_packaging(VX, VY, DX, DY, INTERPMASK, CHIPSIZEX, CHIPSIZEY, SSM, SX, SY, offset2vx_1, offset2vx_2, offset2vy_1, offset2vy_2, MM, VXref, VYref, rangePixelSize, azimuthPixelSize, dt, epsg, srs, tran, out_nc_filename, pair_type, detection_method, coordinates, IMG_INFO_DICT, stable_count, stable_shift_applied, dx_mean_shift, dy_mean_shift, error_vector) + + elif inps.nc_sensor == "L": + + 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]) + + master_path = inps.indir_m + slave_path = inps.indir_s + + import os + master_filename = os.path.basename(master_path) + slave_filename = os.path.basename(slave_path) + + master_split = str.split(master_filename,'_') + slave_split = str.split(slave_filename,'_') + + master_MTL_path = master_path[:-6]+'MTL.txt' + slave_MTL_path = slave_path[:-6]+'MTL.txt' + + master_time = str.split(str.split(runCmd('fgrep "SCENE_CENTER_TIME" '+master_MTL_path))[2][1:-2],':') + slave_time = str.split(str.split(runCmd('fgrep "SCENE_CENTER_TIME" '+slave_MTL_path))[2][1:-2],':') + + from datetime import time as time1 + master_time = time1(int(master_time[0]),int(master_time[1]),int(float(master_time[2]))) + slave_time = time1(int(slave_time[0]),int(slave_time[1]),int(float(slave_time[2]))) + + import netcdf_output as no + version = '1.0.7' + pair_type = 'optical' + detection_method = 'feature' + coordinates = 'map' + # out_nc_filename = 'Jakobshavn_opt.nc' + out_nc_filename = master_filename[0:-4]+'_'+slave_filename[0:-4]+'.nc' + out_nc_filename = './' + out_nc_filename + roi_valid_percentage = int(round(np.sum(CHIPSIZEX!=0)/np.sum(SEARCHLIMITX!=0)*1000.0))/1000 + CHIPSIZEY = np.round(CHIPSIZEX * ScaleChipSizeY / 2) * 2 + + from datetime import date + d0 = date(np.int(master_split[3][0:4]),np.int(master_split[3][4:6]),np.int(master_split[3][6:8])) + d1 = date(np.int(slave_split[3][0:4]),np.int(slave_split[3][4:6]),np.int(slave_split[3][6:8])) + date_dt_base = d1 - d0 + date_dt = np.float64(np.abs(date_dt_base.days)) + if date_dt_base.days < 0: + date_ct = d1 + (d0 - d1)/2 + date_center = date_ct.strftime("%Y%m%d") + else: + date_ct = d0 + (d1 - d0)/2 + date_center = date_ct.strftime("%Y%m%d") + + 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,'roi_valid_percentage':roi_valid_percentage,'autoRIFT_software_version':version} + + error_vector = np.array([57.,57.]) + + no.netCDF_packaging(VX, VY, DX, DY, INTERPMASK, CHIPSIZEX, CHIPSIZEY, SSM, SX, SY, offset2vx_1, offset2vx_2, offset2vy_1, offset2vy_2, MM, VXref, VYref, XPixelSize, YPixelSize, None, epsg, srs, tran, out_nc_filename, pair_type, detection_method, coordinates, IMG_INFO_DICT, stable_count, stable_shift_applied, dx_mean_shift, dy_mean_shift, error_vector) + + elif inps.nc_sensor is None: + print('netCDF packaging not performed') - else: - raise Exception('netCDF packaging not supported for the type "{0}"'.format(inps.nc_sensor)) + else: + raise Exception('netCDF packaging not supported for the type "{0}"'.format(inps.nc_sensor)) print("Write Outputs Done!!!") print(time.time()-t1) diff --git a/testautoRIFT_ISCE.py b/testautoRIFT_ISCE.py index 03313cf..59c3080 100755 --- a/testautoRIFT_ISCE.py +++ b/testautoRIFT_ISCE.py @@ -69,6 +69,8 @@ def cmdLineParse(): help='Input pixel offsets to vx conversion coefficients file name') parser.add_argument('-vy', '--input_vy', dest='offset2vy', type=str, required=False, help='Input pixel offsets to vy conversion coefficients file name') + parser.add_argument('-ssm', '--input_ssm', dest='stable_surface_mask', type=str, required=False, + help='Input stable surface mask file name') parser.add_argument('-fo', '--flag_optical', dest='optical_flag', type=bool, required=False, default=0, help='flag for reading optical data (e.g. Landsat): use 1 for on and 0 (default) for off') parser.add_argument('-nc', '--sensor_flag_netCDF', dest='nc_sensor', type=str, required=False, default=None, @@ -145,14 +147,14 @@ def runAutorift(I1, I2, xGrid, yGrid, Dx0, Dy0, SRx0, SRy0, CSMINx0, CSMINy0, CS obj.I1 = I1 obj.I2 = I2 -# # test with lena image (533 X 533) + # test with lena image (533 X 533) # obj.ChipSizeMinX=16 # obj.ChipSizeMaxX=32 # obj.ChipSize0X=16 # obj.SkipSampleX=16 # obj.SkipSampleY=16 -# # test with Venus image (407 X 407) + # test with Venus image (407 X 407) # obj.ChipSizeMinX=8 # obj.ChipSizeMaxX=16 # obj.ChipSize0X=8 @@ -214,8 +216,11 @@ def runAutorift(I1, I2, xGrid, yGrid, Dx0, Dy0, SRx0, SRy0, CSMINx0, CSMINy0, CS # obj.ChipSize0X = np.min(CSMINx0[CSMINx0!=nodata]) RATIO_Y2X = CSMINy0/CSMINx0 obj.ScaleChipSizeY = np.mean(RATIO_Y2X[(CSMINx0!=nodata)&(CSMINy0!=nodata)]) +# obj.ChipSizeMaxX = obj.ChipSizeMaxX / obj.ChipSizeMaxX * 544 +# obj.ChipSizeMinX = obj.ChipSizeMinX / obj.ChipSizeMinX * 68 else: if ((optflag == 1)&(xGrid is not None)): +# print('test') obj.ChipSizeMaxX = 32 obj.ChipSizeMinX = 16 obj.ChipSize0X = 16 @@ -266,7 +271,7 @@ def runAutorift(I1, I2, xGrid, yGrid, Dx0, Dy0, SRx0, SRy0, CSMINx0, CSMINy0, CS # obj.sparseSearchSampleRate = 16 -# obj.OverSampleRatio = 32 + obj.OverSampleRatio = 64 # 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). @@ -361,7 +366,7 @@ def runAutorift(I1, I2, xGrid, yGrid, Dx0, Dy0, SRx0, SRy0, CSMINx0, CSMINy0, CS if inps.optical_flag == 1: data_m = loadProductOptical(inps.indir_m) data_s = loadProductOptical(inps.indir_s) -# # test with lena/Venus image + # test with lena/Venus image # import scipy.io as sio # conts = sio.loadmat(inps.indir_m) # data_m = conts['I'] @@ -382,6 +387,7 @@ def runAutorift(I1, I2, xGrid, yGrid, Dx0, Dy0, SRx0, SRy0, CSMINx0, CSMINy0, CS CSMINy0 = None CSMAXx0 = None CSMAXy0 = None + SSM = None noDataMask = None nodata = None @@ -435,6 +441,14 @@ def runAutorift(I1, I2, xGrid, yGrid, Dx0, Dy0, SRx0, SRy0, CSMINx0, CSMINy0, CS band=None ds=None + if inps.stable_surface_mask is not None: + ds = gdal.Open(inps.stable_surface_mask) + band = ds.GetRasterBand(1) + SSM = band.ReadAsArray() + SSM = SSM.astype('bool') + band=None + ds=None + Dx, Dy, InterpMask, ChipSizeX, ScaleChipSizeY, SearchLimitX, SearchLimitY, origSize, noDataMask = runAutorift(data_m, data_s, xGrid, yGrid, Dx0, Dy0, SRx0, SRy0, CSMINx0, CSMINy0, CSMAXx0, CSMAXy0, noDataMask, inps.optical_flag, nodata) @@ -462,6 +476,8 @@ def runAutorift(I1, I2, xGrid, yGrid, Dx0, Dy0, SRx0, SRy0, CSMINx0, CSMINy0, CS CHIPSIZEX[noDataMask] = 0 SEARCHLIMITX[noDataMask] = 0 SEARCHLIMITY[noDataMask] = 0 + if SSM is not None: + SSM[noDataMask] = False import scipy.io as sio sio.savemat('offset.mat',{'Dx':DX,'Dy':DY,'InterpMask':INTERPMASK,'ChipSizeX':CHIPSIZEX}) @@ -506,6 +522,7 @@ def runAutorift(I1, I2, xGrid, yGrid, Dx0, Dy0, SRx0, SRy0, CSMINx0, CSMINy0, CS if inps.offset2vx is not None: + ds = gdal.Open(inps.offset2vx) band = ds.GetRasterBand(1) offset2vx_1 = band.ReadAsArray() @@ -521,12 +538,14 @@ def runAutorift(I1, I2, xGrid, yGrid, Dx0, Dy0, SRx0, SRy0, CSMINx0, CSMINy0, CS offset2vy_2 = band.ReadAsArray() band=None ds=None - + VX = offset2vx_1 * DX + offset2vx_2 * DY VY = offset2vy_1 * DX + offset2vy_2 * DY VX = VX.astype(np.float32) VY = VY.astype(np.float32) - + + ############ write velocity output in Geotiff format + outRaster = driver.Create("velocity.tif", int(xGrid.shape[1]), int(xGrid.shape[0]), 2, gdal.GDT_Float32) outRaster.SetGeoTransform(tran) outRaster.SetProjection(proj) @@ -537,96 +556,194 @@ def runAutorift(I1, I2, xGrid, yGrid, Dx0, Dy0, SRx0, SRy0, CSMINx0, CSMINy0, CS outband.WriteArray(VY) outband.FlushCache() - ######################################################################################## - ############ netCDF packaging for Sentinel and Landsat dataset; can add other sensor format as well - if inps.nc_sensor == "S": + ############ prepare for netCDF packaging + + if inps.nc_sensor is not None: + + vxrefname = str.split(runCmd('fgrep "Velocities:" testGeogrid.txt'))[1] + vyrefname = str.split(runCmd('fgrep "Velocities:" testGeogrid.txt'))[2] + sxname = str.split(runCmd('fgrep "Slopes:" testGeogrid.txt'))[1][:-4]+"s.tif" + syname = str.split(runCmd('fgrep "Slopes:" testGeogrid.txt'))[2][:-4]+"s.tif" + maskname = str.split(runCmd('fgrep "Slopes:" testGeogrid.txt'))[2][:-8]+"masks.tif" + xoff = int(str.split(runCmd('fgrep "Origin index (in DEM) of geogrid:" testGeogrid.txt'))[6]) + yoff = int(str.split(runCmd('fgrep "Origin index (in DEM) of geogrid:" testGeogrid.txt'))[7]) + xcount = int(str.split(runCmd('fgrep "Dimensions of geogrid:" testGeogrid.txt'))[3]) + ycount = int(str.split(runCmd('fgrep "Dimensions of geogrid:" testGeogrid.txt'))[5]) + + ds = gdal.Open(vxrefname) + band = ds.GetRasterBand(1) + VXref = band.ReadAsArray(xoff, yoff, xcount, ycount) + ds = None + band = None - 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]) - epsg = float(str.split(runCmd('fgrep "EPSG:" testGeogrid.txt'))[1]) -# print (str(rangePixelSize)+" "+str(azimuthPixelSize)) - - runCmd('topsinsar_filename.py') -# import scipy.io as sio - conts = sio.loadmat('topsinsar_filename.mat') - master_filename = conts['master_filename'][0] - slave_filename = conts['slave_filename'][0] - master_split = str.split(master_filename,'_') - slave_split = str.split(slave_filename,'_') - - import netcdf_output as no - version = '1.0.5' - pair_type = 'radar' - detection_method = 'feature' - coordinates = 'radar' -# out_nc_filename = 'Jakobshavn.nc' - out_nc_filename = master_filename[0:-4]+'_'+slave_filename[0:-4]+'.nc' - out_nc_filename = './' + out_nc_filename - roi_valid_percentage = np.sum(CHIPSIZEX!=0)/np.sum(SEARCHLIMITX!=0)*100.0 - CHIPSIZEY = np.round(CHIPSIZEX * ScaleChipSizeY / 2) * 2 + ds = gdal.Open(vyrefname) + band = ds.GetRasterBand(1) + VYref = band.ReadAsArray(xoff, yoff, xcount, ycount) + ds = None + band = None - + ds = gdal.Open(sxname) + band = ds.GetRasterBand(1) + SX = band.ReadAsArray(xoff, yoff, xcount, ycount) + ds = None + band = None - from datetime import date - d0 = date(np.int(master_split[5][0:4]),np.int(master_split[5][4:6]),np.int(master_split[5][6:8])) - d1 = date(np.int(slave_split[5][0:4]),np.int(slave_split[5][4:6]),np.int(slave_split[5][6:8])) - date_dt_base = d1 - d0 - date_dt = np.float64(np.abs(date_dt_base.days)) - if date_dt_base.days < 0: - date_ct = d1 + (d0 - d1)/2 - date_center = date_ct.strftime("%Y%m%d") - else: - 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_split[5][0:8],'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_split[5][0:8],'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,'roi_valid_percentage':roi_valid_percentage,'autoRIFT_software_version':version} - - no.netCDF_packaging(VX, VY, DX, DY, INTERPMASK, CHIPSIZEX, CHIPSIZEY, rangePixelSize, azimuthPixelSize, dt, epsg, srs, tran, out_nc_filename, pair_type, detection_method, coordinates, IMG_INFO_DICT) - - elif inps.nc_sensor == "L": + ds = gdal.Open(syname) + band = ds.GetRasterBand(1) + SY = band.ReadAsArray(xoff, yoff, xcount, ycount) + ds = None + band = None + + ds = gdal.Open(maskname) + band = ds.GetRasterBand(1) + MM = band.ReadAsArray(xoff, yoff, xcount, ycount) + ds = None + band = None - 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]) + DXref = offset2vy_2 / (offset2vx_1 * offset2vy_2 - offset2vx_2 * offset2vy_1) * VXref - offset2vx_2 / (offset2vx_1 * offset2vy_2 - offset2vx_2 * offset2vy_1) * VYref + DYref = offset2vx_1 / (offset2vx_1 * offset2vy_2 - offset2vx_2 * offset2vy_1) * VYref - offset2vy_1 / (offset2vx_1 * offset2vy_2 - offset2vx_2 * offset2vy_1) * VXref - master_filename = inps.indir_m - slave_filename = inps.indir_s - master_split = str.split(master_filename,'_') - slave_split = str.split(slave_filename,'_') + stable_count = np.sum(SSM & np.logical_not(np.isnan(DX))) - import netcdf_output as no - version = '1.0.5' - pair_type = 'optical' - detection_method = 'feature' - coordinates = 'map' -# out_nc_filename = 'Jakobshavn_opt.nc' - out_nc_filename = master_filename[0:-4]+'_'+slave_filename[0:-4]+'.nc' - out_nc_filename = './' + out_nc_filename - roi_valid_percentage = np.sum(CHIPSIZEX!=0)/np.sum(SEARCHLIMITX!=0)*100.0 - CHIPSIZEY = np.round(CHIPSIZEX * ScaleChipSizeY / 2) * 2 - - from datetime import date - d0 = date(np.int(master_split[3][0:4]),np.int(master_split[3][4:6]),np.int(master_split[3][6:8])) - d1 = date(np.int(slave_split[3][0:4]),np.int(slave_split[3][4:6]),np.int(slave_split[3][6:8])) - date_dt_base = d1 - d0 - date_dt = np.float64(np.abs(date_dt_base.days)) - if date_dt_base.days < 0: - date_ct = d1 + (d0 - d1)/2 - date_center = date_ct.strftime("%Y%m%d") + if stable_count == 0: + stable_shift_applied = 0 else: - 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':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_split[3][0:8],'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_split[3][0:8],'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,'roi_valid_percentage':roi_valid_percentage,'autoRIFT_software_version':version} + stable_shift_applied = 1 + + if stable_shift_applied == 1: + temp = DX.copy() - DXref.copy() + temp[np.logical_not(SSM)] = np.nan + dx_mean_shift = np.median(temp[(temp > -5)&(temp < 5)]) + DX = DX - dx_mean_shift + + temp = DY.copy() - DYref.copy() + temp[np.logical_not(SSM)] = np.nan + dy_mean_shift = np.median(temp[(temp > -5)&(temp < 5)]) + DY = DY - dy_mean_shift + else: + dx_mean_shift = 0.0 + dy_mean_shift = 0.0 - no.netCDF_packaging(VX, VY, DX, DY, INTERPMASK, CHIPSIZEX, CHIPSIZEY, XPixelSize, YPixelSize, None, epsg, srs, tran, out_nc_filename, pair_type, detection_method, coordinates, IMG_INFO_DICT) + VX = offset2vx_1 * DX + offset2vx_2 * DY + VY = offset2vy_1 * DX + offset2vy_2 * DY + VX = VX.astype(np.float32) + VY = VY.astype(np.float32) - elif inps.nc_sensor is None: - print('netCDF packaging not performed') + ######################################################################################## + ############ netCDF packaging for Sentinel and Landsat dataset; can add other sensor format as well + if inps.nc_sensor == "S": + + 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]) + epsg = float(str.split(runCmd('fgrep "EPSG:" testGeogrid.txt'))[1]) + # print (str(rangePixelSize)+" "+str(azimuthPixelSize)) + + runCmd('topsinsar_filename.py') + # import scipy.io as sio + conts = sio.loadmat('topsinsar_filename.mat') + master_filename = conts['master_filename'][0] + slave_filename = conts['slave_filename'][0] + master_dt = conts['master_dt'][0] + slave_dt = conts['slave_dt'][0] + master_split = str.split(master_filename,'_') + slave_split = str.split(slave_filename,'_') + + import netcdf_output as no + version = '1.0.7' + pair_type = 'radar' + detection_method = 'feature' + coordinates = 'radar' + # out_nc_filename = 'Jakobshavn.nc' + out_nc_filename = master_filename[0:-4]+'_'+slave_filename[0:-4]+'.nc' + out_nc_filename = './' + out_nc_filename + roi_valid_percentage = int(round(np.sum(CHIPSIZEX!=0)/np.sum(SEARCHLIMITX!=0)*1000.0))/1000 + CHIPSIZEY = np.round(CHIPSIZEX * ScaleChipSizeY / 2) * 2 + + + + from datetime import date + d0 = date(np.int(master_split[5][0:4]),np.int(master_split[5][4:6]),np.int(master_split[5][6:8])) + d1 = date(np.int(slave_split[5][0:4]),np.int(slave_split[5][4:6]),np.int(slave_split[5][6:8])) + date_dt_base = d1 - d0 + date_dt = np.float64(np.abs(date_dt_base.days)) + if date_dt_base.days < 0: + date_ct = d1 + (d0 - d1)/2 + date_center = date_ct.strftime("%Y%m%d") + else: + 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,'roi_valid_percentage':roi_valid_percentage,'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]]) + + no.netCDF_packaging(VX, VY, DX, DY, INTERPMASK, CHIPSIZEX, CHIPSIZEY, SSM, SX, SY, offset2vx_1, offset2vx_2, offset2vy_1, offset2vy_2, MM, VXref, VYref, rangePixelSize, azimuthPixelSize, dt, epsg, srs, tran, out_nc_filename, pair_type, detection_method, coordinates, IMG_INFO_DICT, stable_count, stable_shift_applied, dx_mean_shift, dy_mean_shift, error_vector) + + elif inps.nc_sensor == "L": + + 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]) + + master_path = inps.indir_m + slave_path = inps.indir_s + + import os + master_filename = os.path.basename(master_path) + slave_filename = os.path.basename(slave_path) + + master_split = str.split(master_filename,'_') + slave_split = str.split(slave_filename,'_') + + master_MTL_path = master_path[:-6]+'MTL.txt' + slave_MTL_path = slave_path[:-6]+'MTL.txt' + + master_time = str.split(str.split(runCmd('fgrep "SCENE_CENTER_TIME" '+master_MTL_path))[2][1:-2],':') + slave_time = str.split(str.split(runCmd('fgrep "SCENE_CENTER_TIME" '+slave_MTL_path))[2][1:-2],':') + + from datetime import time as time1 + master_time = time1(int(master_time[0]),int(master_time[1]),int(float(master_time[2]))) + slave_time = time1(int(slave_time[0]),int(slave_time[1]),int(float(slave_time[2]))) + + import netcdf_output as no + version = '1.0.7' + pair_type = 'optical' + detection_method = 'feature' + coordinates = 'map' + # out_nc_filename = 'Jakobshavn_opt.nc' + out_nc_filename = master_filename[0:-4]+'_'+slave_filename[0:-4]+'.nc' + out_nc_filename = './' + out_nc_filename + roi_valid_percentage = int(round(np.sum(CHIPSIZEX!=0)/np.sum(SEARCHLIMITX!=0)*1000.0))/1000 + CHIPSIZEY = np.round(CHIPSIZEX * ScaleChipSizeY / 2) * 2 + + from datetime import date + d0 = date(np.int(master_split[3][0:4]),np.int(master_split[3][4:6]),np.int(master_split[3][6:8])) + d1 = date(np.int(slave_split[3][0:4]),np.int(slave_split[3][4:6]),np.int(slave_split[3][6:8])) + date_dt_base = d1 - d0 + date_dt = np.float64(np.abs(date_dt_base.days)) + if date_dt_base.days < 0: + date_ct = d1 + (d0 - d1)/2 + date_center = date_ct.strftime("%Y%m%d") + else: + date_ct = d0 + (d1 - d0)/2 + date_center = date_ct.strftime("%Y%m%d") + + 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,'roi_valid_percentage':roi_valid_percentage,'autoRIFT_software_version':version} + + error_vector = np.array([57.,57.]) + + no.netCDF_packaging(VX, VY, DX, DY, INTERPMASK, CHIPSIZEX, CHIPSIZEY, SSM, SX, SY, offset2vx_1, offset2vx_2, offset2vy_1, offset2vy_2, MM, VXref, VYref, XPixelSize, YPixelSize, None, epsg, srs, tran, out_nc_filename, pair_type, detection_method, coordinates, IMG_INFO_DICT, stable_count, stable_shift_applied, dx_mean_shift, dy_mean_shift, error_vector) + + elif inps.nc_sensor is None: + print('netCDF packaging not performed') - else: - raise Exception('netCDF packaging not supported for the type "{0}"'.format(inps.nc_sensor)) + else: + raise Exception('netCDF packaging not supported for the type "{0}"'.format(inps.nc_sensor)) print("Write Outputs Done!!!") print(time.time()-t1)