diff --git a/geo_autoRIFT/geogrid/GeogridOptical.py b/geo_autoRIFT/geogrid/GeogridOptical.py index c09ea09..54aea4c 100755 --- a/geo_autoRIFT/geogrid/GeogridOptical.py +++ b/geo_autoRIFT/geogrid/GeogridOptical.py @@ -50,11 +50,11 @@ def runGeogrid(self): ##Determine appropriate EPSG system self.epsgDem = self.getProjectionSystem(self.demname) self.epsgDat = self.getProjectionSystem(self.dat1name) - + ###Determine extent of data needed bbox = self.determineBbox() - + ##Run self.geogrid() @@ -102,7 +102,7 @@ def determineBbox(self, zrange=[-200,4000]): import numpy as np import datetime from osgeo import osr - + # import pdb # pdb.set_trace() @@ -115,18 +115,18 @@ def determineBbox(self, zrange=[-200,4000]): coordDat.ImportFromEPSG(self.epsgDat) else: raise Exception('EPSG code does not exist for image data') - + coordDem = osr.SpatialReference() if self.epsgDem: coordDem.ImportFromEPSG(self.epsgDem) else: raise Exception('EPSG code does not exist for DEM') - - + + trans = osr.CoordinateTransformation(coordDat, coordDem) - - + + utms = [] xyzs = [] @@ -145,21 +145,21 @@ def determineBbox(self, zrange=[-200,4000]): self._xlim = [np.min(xyzs[:,0]), np.max(xyzs[:,0])] self._ylim = [np.min(xyzs[:,1]), np.max(xyzs[:,1])] - - - - - + + + + + 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)) @@ -190,10 +190,10 @@ def geogrid(self): 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)) @@ -206,53 +206,53 @@ def geogrid(self): print("Output Nodata Value: " + str(self.nodata_out) + "\n") - - + + print("Starting processing .... ") - - - - + + + + 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)) @@ -276,7 +276,7 @@ def geogrid(self): 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)) @@ -284,8 +284,8 @@ def geogrid(self): 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() @@ -297,26 +297,30 @@ def geogrid(self): 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) - + if (self.vxname != ""): nodata = vxDS.GetRasterBand(1).GetNoDataValue() else: @@ -352,12 +356,12 @@ 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) @@ -368,12 +372,12 @@ def geogrid(self): 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) @@ -383,27 +387,27 @@ def geogrid(self): 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) @@ -414,12 +418,12 @@ def geogrid(self): 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) @@ -430,12 +434,12 @@ def geogrid(self): 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) @@ -445,12 +449,12 @@ def geogrid(self): 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) @@ -474,18 +478,20 @@ def geogrid(self): 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 - - + + @@ -493,31 +499,31 @@ def geogrid(self): 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) @@ -543,13 +549,13 @@ def geogrid(self): 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 @@ -568,7 +574,7 @@ def geogrid(self): 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)): @@ -577,7 +583,7 @@ def geogrid(self): 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 @@ -585,7 +591,7 @@ def geogrid(self): 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 @@ -614,7 +620,7 @@ def geogrid(self): 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])); @@ -662,7 +668,7 @@ def geogrid(self): # 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 != "")): @@ -685,7 +691,7 @@ def geogrid(self): 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 @@ -699,23 +705,23 @@ def geogrid(self): poDstDSMsk = None if (self.dhdxname != ""): poDstDSRO2VX = None - + poDstDSRO2VY = None - + demDS = None - + if (self.dhdxname != ""): sxDS = None syDS = None - + if (self.vxname != ""): vxDS = None vyDS = None - + if (self.srxname != ""): srxDS = None sryDS = None - + if (self.csminxname != ""): csminxDS = None csminyDS = None @@ -723,7 +729,7 @@ def geogrid(self): if (self.csmaxxname != ""): csmaxxDS = None csmaxyDS = None - + if (self.ssmname != ""): ssmDS = None @@ -735,7 +741,7 @@ def geogrid(self): def coregister(self,in1,in2): import os import numpy as np - + from osgeo import gdal, osr import struct @@ -748,22 +754,22 @@ def coregister(self,in1,in2): trans2 = DS2.GetGeoTransform() xsize2 = DS2.RasterXSize ysize2 = DS2.RasterYSize - + W = np.max([trans1[0],trans2[0]]) N = np.min([trans1[3],trans2[3]]) E = np.min([trans1[0]+xsize1*trans1[1],trans2[0]+xsize2*trans2[1]]) S = np.max([trans1[3]+ysize1*trans1[5],trans2[3]+ysize2*trans2[5]]) - + x1a = int(np.round((W-trans1[0])/trans1[1])) x1b = int(np.round((E-trans1[0])/trans1[1])) y1a = int(np.round((N-trans1[3])/trans1[5])) y1b = int(np.round((S-trans1[3])/trans1[5])) - + x2a = int(np.round((W-trans2[0])/trans2[1])) x2b = int(np.round((E-trans2[0])/trans2[1])) y2a = int(np.round((N-trans2[3])/trans2[5])) y2b = int(np.round((S-trans2[3])/trans2[5])) - + x1a = np.min([x1a, xsize1-1]) x1b = np.min([x1b, xsize1-1]) y1a = np.min([y1a, ysize1-1]) @@ -772,7 +778,7 @@ def coregister(self,in1,in2): x2b = np.min([x2b, xsize2-1]) y2a = np.min([y2a, ysize2-1]) y2b = np.min([y2b, ysize2-1]) - + x1a = np.max([x1a, 0]) x1b = np.max([x1b, 0]) y1a = np.max([y1a, 0]) @@ -781,7 +787,7 @@ def coregister(self,in1,in2): x2b = np.max([x2b, 0]) y2a = np.max([y2a, 0]) y2b = np.max([y2b, 0]) - + x1a = int(x1a) x1b = int(x1b) y1a = int(y1a) @@ -798,9 +804,9 @@ def coregister(self,in1,in2): - - - + + + def __init__(self): super(GeogridOptical, self).__init__() @@ -829,7 +835,7 @@ def __init__(self): self.csmaxxname = None self.csmaxyname = None self.ssmname = None - + ##Output related parameters self.winlocname = None self.winoffname = None @@ -847,4 +853,10 @@ def __init__(self): self._ylim = None self.nodata_out = None - + ##parameters for autoRIFT + self.pOff = None + self.lOff = None + self.pCount = None + self.lCount = None + self.X_res = None + self.Y_res = None diff --git a/testautoRIFT.py b/testautoRIFT.py index 1bd59ef..6352a70 100755 --- a/testautoRIFT.py +++ b/testautoRIFT.py @@ -92,7 +92,7 @@ def loadProduct(filename): import isce import logging from imageMath import IML - + IMG = IML.mmapFromISCE(filename, logging) img = IMG.bands[0] # pdb.set_trace() @@ -104,30 +104,30 @@ def loadProductOptical(file_m, file_s): ''' Load the product using Product Manager. ''' - from geogrid import GeogridOptical # from components.contrib.geo_autoRIFT.geogrid import GeogridOptical obj = GeogridOptical() - + x1a, y1a, xsize1, ysize1, x2a, y2a, xsize2, ysize2, trans = obj.coregister(file_m, file_s) - + DS1 = gdal.Open(file_m) DS2 = gdal.Open(file_s) - + I1 = DS1.ReadAsArray(xoff=x1a, yoff=y1a, xsize=xsize1, ysize=ysize1) I2 = DS2.ReadAsArray(xoff=x2a, yoff=y2a, xsize=xsize2, ysize=ysize2) - + I1 = I1.astype(np.float32) I2 = I2.astype(np.float32) - + DS1=None DS2=None - + return I1, I2 -def runAutorift(I1, I2, xGrid, yGrid, Dx0, Dy0, SRx0, SRy0, CSMINx0, CSMINy0, CSMAXx0, CSMAXy0, noDataMask, optflag, nodata, mpflag): +def runAutorift(I1, I2, xGrid, yGrid, Dx0, Dy0, SRx0, SRy0, CSMINx0, CSMINy0, CSMAXx0, CSMAXy0, noDataMask, optflag, + nodata, mpflag, geogrid_run_info=None): ''' Wire and run geogrid. ''' @@ -139,8 +139,8 @@ def runAutorift(I1, I2, xGrid, yGrid, Dx0, Dy0, SRx0, SRy0, CSMINx0, CSMINy0, CS # import isceobj import time import subprocess - - + + obj = autoRIFT() # obj.configure() @@ -149,7 +149,7 @@ def runAutorift(I1, I2, xGrid, yGrid, Dx0, Dy0, SRx0, SRy0, CSMINx0, CSMINy0, CS # I2 = I2.astype(np.uint8) obj.MultiThread = mpflag - + # take the amplitude only for the radar images if optflag == 0: I1 = np.abs(I1) @@ -171,8 +171,11 @@ def runAutorift(I1, I2, xGrid, yGrid, Dx0, Dy0, SRx0, SRy0, CSMINx0, CSMINy0, CS # obj.ChipSize0X=8 # obj.SkipSampleX=8 # obj.SkipSampleY=8 - - + + # test with small tiff images +# obj.SkipSampleX=16 +# obj.SkipSampleY=16 + # create the grid if it does not exist if xGrid is None: m,n = obj.I1.shape @@ -188,7 +191,7 @@ def runAutorift(I1, I2, xGrid, yGrid, Dx0, Dy0, SRx0, SRy0, CSMINx0, CSMINy0, CS obj.yGrid = yGrid - + # generate the nodata mask where offset searching will be skipped based on 1) imported nodata mask and/or 2) zero values in the image for ii in range(obj.xGrid.shape[0]): for jj in range(obj.xGrid.shape[1]): @@ -218,11 +221,15 @@ def runAutorift(I1, I2, xGrid, yGrid, Dx0, Dy0, SRx0, SRy0, CSMINx0, CSMINy0, CS if CSMINx0 is not None: obj.ChipSizeMaxX = CSMAXx0 obj.ChipSizeMinX = CSMINx0 - 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]) - except: - pixsizex = float(str.split(subprocess.getoutput('fgrep "X-direction pixel size:" testGeogrid.txt'))[-1]) + if geogrid_run_info is None: + 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]) + except: + pixsizex = float(str.split(subprocess.getoutput('fgrep "X-direction pixel size:" testGeogrid.txt'))[-1]) + else: + chipsizex0 = geogrid_run_info['chipsizex0'] + pixsizex = geogrid_run_info['XPixelSize'] obj.ChipSize0X = int(np.ceil(chipsizex0/pixsizex/4)*4) # obj.ChipSize0X = np.min(CSMINx0[CSMINx0!=nodata]) RATIO_Y2X = CSMINy0/CSMINx0 @@ -278,7 +285,6 @@ def runAutorift(I1, I2, xGrid, yGrid, Dx0, Dy0, SRx0, SRy0, CSMINx0, CSMINy0, CS print("Uniform Data Type Done!!!") print(time.time()-t1) - # pdb.set_trace() # obj.sparseSearchSampleRate = 16 @@ -296,9 +302,6 @@ def runAutorift(I1, I2, xGrid, yGrid, Dx0, Dy0, SRx0, SRy0, CSMINx0, CSMINy0, CS - - - # ########## export preprocessed images to files; can be commented out if not debugging # # t1 = time.time() @@ -366,9 +369,6 @@ def runAutorift(I1, I2, xGrid, yGrid, Dx0, Dy0, SRx0, SRy0, CSMINx0, CSMINy0, CS - - - def main(): ''' Main driver. @@ -384,7 +384,8 @@ def main(): def generateAutoriftProduct(indir_m, indir_s, grid_location, init_offset, search_range, chip_size_min, chip_size_max, - offset2vx, offset2vy, stable_surface_mask, optical_flag, nc_sensor, mpflag): + offset2vx, offset2vy, stable_surface_mask, optical_flag, nc_sensor, mpflag, + geogrid_run_info=None): import numpy as np import time @@ -480,7 +481,8 @@ def generateAutoriftProduct(indir_m, indir_s, grid_location, init_offset, search 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, optical_flag, nodata, mpflag + data_m, data_s, xGrid, yGrid, Dx0, Dy0, SRx0, SRy0, CSMINx0, CSMINy0, CSMAXx0, CSMAXy0, + noDataMask, optical_flag, nodata, mpflag, geogrid_run_info=geogrid_run_info, ) if optical_flag == 0: @@ -492,7 +494,7 @@ def generateAutoriftProduct(indir_m, indir_s, grid_location, init_offset, search CHIPSIZEX = np.zeros(origSize,dtype=np.float32) SEARCHLIMITX = np.zeros(origSize,dtype=np.float32) SEARCHLIMITY = np.zeros(origSize,dtype=np.float32) - + DX[0:Dx.shape[0],0:Dx.shape[1]] = Dx DY[0:Dy.shape[0],0:Dy.shape[1]] = Dy INTERPMASK[0:InterpMask.shape[0],0:InterpMask.shape[1]] = InterpMask @@ -524,6 +526,7 @@ def generateAutoriftProduct(indir_m, indir_s, grid_location, init_offset, search # SEARCHLIMITY = conts['SearchLimitY'] # ##################### + netcdf_file = None if grid_location is not None: @@ -573,7 +576,7 @@ def generateAutoriftProduct(indir_m, indir_s, grid_location, init_offset, search 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) @@ -585,20 +588,30 @@ def generateAutoriftProduct(indir_m, indir_s, grid_location, init_offset, search outband = outRaster.GetRasterBand(2) outband.WriteArray(VY) outband.FlushCache() - + ############ prepare for netCDF packaging if 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]) + if geogrid_run_info is 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]+"sp.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]) + else: + vxrefname = geogrid_run_info['vxname'] + vyrefname = geogrid_run_info['vyname'] + sxname = geogrid_run_info['sxname'] + syname = geogrid_run_info['syname'] + maskname = geogrid_run_info['maskname'] + xoff = geogrid_run_info['xoff'] + yoff = geogrid_run_info['yoff'] + xcount = geogrid_run_info['xcount'] + ycount = geogrid_run_info['ycount'] ds = gdal.Open(vxrefname) band = ds.GetRasterBand(1) @@ -629,23 +642,23 @@ def generateAutoriftProduct(indir_m, indir_s, grid_location, init_offset, search MM = band.ReadAsArray(xoff, yoff, xcount, ycount) ds = None band = None - + 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)) & (DX-DXref > -5) & (DX-DXref < 5) & (DY-DYref > -5) & (DY-DYref < 5)) if stable_count == 0: stable_shift_applied = 0 else: 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)]) @@ -653,21 +666,28 @@ def generateAutoriftProduct(indir_m, indir_s, grid_location, init_offset, search 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 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)) - + if geogrid_run_info is None: + chipsizex0 = float(str.split(runCmd('fgrep "Smallest Allowable Chip Size 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]) + epsg = float(str.split(runCmd('fgrep "EPSG:" testGeogrid.txt'))[1]) + # print (str(rangePixelSize)+" "+str(azimuthPixelSize)) + else: + chipsizex0 = geogrid_run_info['chipsizex0'] + rangePixelSize = geogrid_run_info['XPixelSize'] + azimuthPixelSize = geogrid_run_info['YPixelSize'] + dt = geogrid_run_info['dt'] + epsg = geogrid_run_info['epsg'] runCmd('topsinsar_filename.py') # import scipy.io as sio @@ -686,9 +706,8 @@ def generateAutoriftProduct(indir_m, indir_s, grid_location, init_offset, search roi_valid_percentage = int(round(np.sum(CHIPSIZEX!=0)/np.sum(SEARCHLIMITX!=0)*1000.0))/1000 # out_nc_filename = 'Jakobshavn.nc' PPP = roi_valid_percentage * 100 - out_nc_filename = master_filename[0:-4]+'_X_'+slave_filename[0:-4]+'_G0240V02_P{0}{1}{2}.nc'.format(int(PPP/10),int((PPP-int(PPP/10)*10)),int((PPP-int(PPP/10)*10-int((PPP-int(PPP/10)*10)))*10)) - out_nc_filename = './' + out_nc_filename - + 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" CHIPSIZEY = np.round(CHIPSIZEX * ScaleChipSizeY / 2) * 2 @@ -709,13 +728,25 @@ def generateAutoriftProduct(indir_m, indir_s, grid_location, init_offset, search 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) + netcdf_file = 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 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]) + if geogrid_run_info is None: + chipsizex0 = float(str.split(runCmd('fgrep "Smallest Allowable Chip Size 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'] + XPixelSize = geogrid_run_info['XPixelSize'] + YPixelSize = geogrid_run_info['YPixelSize'] + epsg = geogrid_run_info['epsg'] master_path = indir_m slave_path = indir_s @@ -723,10 +754,10 @@ def generateAutoriftProduct(indir_m, indir_s, grid_location, init_offset, search 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' # @@ -738,7 +769,7 @@ def generateAutoriftProduct(indir_m, indir_s, grid_location, init_offset, search 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 pair_type = 'optical' detection_method = 'feature' @@ -746,8 +777,9 @@ def generateAutoriftProduct(indir_m, indir_s, grid_location, init_offset, search roi_valid_percentage = int(round(np.sum(CHIPSIZEX!=0)/np.sum(SEARCHLIMITX!=0)*1000.0))/1000 # out_nc_filename = 'Jakobshavn_opt.nc' PPP = roi_valid_percentage * 100 - out_nc_filename = master_filename[0:-8]+'_X_'+slave_filename[0:-8]+'_G0240V02_P{0}{1}{2}.nc'.format(int(PPP/10),int((PPP-int(PPP/10)*10)),int((PPP-int(PPP/10)*10-int((PPP-int(PPP/10)*10)))*10)) - out_nc_filename = './' + out_nc_filename + out_nc_filename = f"./{master_filename[0:-7]}_X_{slave_filename[0:-7]}" \ + f"_G{chipsizex0:04.0f}V02_P{np.floor(PPP):03.0f}.nc" + CHIPSIZEY = np.round(CHIPSIZEX * ScaleChipSizeY / 2) * 2 from datetime import date @@ -761,50 +793,62 @@ def generateAutoriftProduct(indir_m, indir_s, grid_location, init_offset, search 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 nc_sensor == "S2": + netcdf_file = 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 + ) - 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]) + 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]) + 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'] + XPixelSize = geogrid_run_info['XPixelSize'] + YPixelSize = geogrid_run_info['YPixelSize'] + epsg = geogrid_run_info['epsg'] master_path = indir_m slave_path = indir_s master_split = master_path.split('_') slave_split = slave_path.split('_') - + import os - + master_filename = master_split[0][-3:]+'_'+master_split[2]+'_'+master_split[4][:3]+'_'+os.path.basename(master_path) slave_filename = slave_split[0][-3:]+'_'+slave_split[2]+'_'+slave_split[4][:3]+'_'+os.path.basename(slave_path) - + master_time = ['00','00','00'] slave_time = ['00','00','00'] - + 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 pair_type = 'optical' detection_method = 'feature' coordinates = 'map' roi_valid_percentage = int(round(np.sum(CHIPSIZEX!=0)/np.sum(SEARCHLIMITX!=0)*1000.0))/1000 PPP = roi_valid_percentage * 100 - out_nc_filename = master_filename[0:-8]+'_X_'+slave_filename[0:-8]+'_G0240V02_P{0}{1}{2}.nc'.format(int(PPP/10),int((PPP-int(PPP/10)*10)),int((PPP-int(PPP/10)*10-int((PPP-int(PPP/10)*10)))*10)) - out_nc_filename = './' + out_nc_filename + 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" CHIPSIZEY = np.round(CHIPSIZEX * ScaleChipSizeY / 2) * 2 - + from datetime import date d0 = date(np.int(master_split[2][0:4]),np.int(master_split[2][4:6]),np.int(master_split[2][6:8])) d1 = date(np.int(slave_split[2][0:4]),np.int(slave_split[2][4:6]),np.int(slave_split[2][6:8])) @@ -824,7 +868,13 @@ def generateAutoriftProduct(indir_m, indir_s, grid_location, init_offset, search 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) + netcdf_file = 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 nc_sensor is None: print('netCDF packaging not performed') @@ -835,6 +885,8 @@ def generateAutoriftProduct(indir_m, indir_s, grid_location, init_offset, search print("Write Outputs Done!!!") print(time.time()-t1) + return netcdf_file + if __name__ == '__main__': main() diff --git a/testautoRIFT_ISCE.py b/testautoRIFT_ISCE.py index 9b0a96f..6b17f55 100755 --- a/testautoRIFT_ISCE.py +++ b/testautoRIFT_ISCE.py @@ -77,6 +77,7 @@ def cmdLineParse(): help='flag for packaging output formatted for Sentinel ("S") and Landsat ("L") dataset; default is None') parser.add_argument('-mpflag', '--mpflag', dest='mpflag', type=int, required=False, default=0, help='number of threads for multiple threading (default is specified by 0, which uses the original single-core version and surpasses the multithreading routine)') + return parser.parse_args() class Dummy(object): @@ -91,7 +92,7 @@ def loadProduct(filename): import isce import logging from imageMath import IML - + IMG = IML.mmapFromISCE(filename, logging) img = IMG.bands[0] # pdb.set_trace() @@ -103,29 +104,31 @@ def loadProductOptical(file_m, file_s): ''' Load the product using Product Manager. ''' + import isce from components.contrib.geo_autoRIFT.geogrid import GeogridOptical # from geogrid import GeogridOptical obj = GeogridOptical() - + x1a, y1a, xsize1, ysize1, x2a, y2a, xsize2, ysize2, trans = obj.coregister(file_m, file_s) - + DS1 = gdal.Open(file_m) DS2 = gdal.Open(file_s) - + I1 = DS1.ReadAsArray(xoff=x1a, yoff=y1a, xsize=xsize1, ysize=ysize1) I2 = DS2.ReadAsArray(xoff=x2a, yoff=y2a, xsize=xsize2, ysize=ysize2) - + I1 = I1.astype(np.float32) I2 = I2.astype(np.float32) - + DS1=None DS2=None - + return I1, I2 -def runAutorift(I1, I2, xGrid, yGrid, Dx0, Dy0, SRx0, SRy0, CSMINx0, CSMINy0, CSMAXx0, CSMAXy0, noDataMask, optflag, nodata, mpflag): +def runAutorift(I1, I2, xGrid, yGrid, Dx0, Dy0, SRx0, SRy0, CSMINx0, CSMINy0, CSMAXx0, CSMAXy0, noDataMask, optflag, + nodata, mpflag, geogrid_run_info=None): ''' Wire and run geogrid. ''' @@ -136,12 +139,12 @@ def runAutorift(I1, I2, xGrid, yGrid, Dx0, Dy0, SRx0, SRy0, CSMINx0, CSMINy0, CS import isceobj import time import subprocess - - + + obj = autoRIFT_ISCE() obj.configure() - ########## uncomment if starting from preprocessed images +########## uncomment if starting from preprocessed images # I1 = I1.astype(np.uint8) # I2 = I2.astype(np.uint8) @@ -172,7 +175,7 @@ def runAutorift(I1, I2, xGrid, yGrid, Dx0, Dy0, SRx0, SRy0, CSMINx0, CSMINy0, CS # test with small tiff images # obj.SkipSampleX=16 # obj.SkipSampleY=16 - + # create the grid if it does not exist if xGrid is None: m,n = obj.I1.shape @@ -188,7 +191,7 @@ def runAutorift(I1, I2, xGrid, yGrid, Dx0, Dy0, SRx0, SRy0, CSMINx0, CSMINy0, CS obj.yGrid = yGrid - + # generate the nodata mask where offset searching will be skipped based on 1) imported nodata mask and/or 2) zero values in the image for ii in range(obj.xGrid.shape[0]): for jj in range(obj.xGrid.shape[1]): @@ -218,11 +221,15 @@ def runAutorift(I1, I2, xGrid, yGrid, Dx0, Dy0, SRx0, SRy0, CSMINx0, CSMINy0, CS if CSMINx0 is not None: obj.ChipSizeMaxX = CSMAXx0 obj.ChipSizeMinX = CSMINx0 - 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]) - except: - pixsizex = float(str.split(subprocess.getoutput('fgrep "X-direction pixel size:" testGeogrid.txt'))[-1]) + if geogrid_run_info is None: + 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]) + except: + pixsizex = float(str.split(subprocess.getoutput('fgrep "X-direction pixel size:" testGeogrid.txt'))[-1]) + else: + chipsizex0 = geogrid_run_info['chipsizex0'] + pixsizex = geogrid_run_info['XPixelSize'] obj.ChipSize0X = int(np.ceil(chipsizex0/pixsizex/4)*4) # obj.ChipSize0X = np.min(CSMINx0[CSMINx0!=nodata]) RATIO_Y2X = CSMINy0/CSMINx0 @@ -265,7 +272,10 @@ def runAutorift(I1, I2, xGrid, yGrid, Dx0, Dy0, SRx0, SRy0, CSMINx0, CSMINy0, CS ######## preprocessing t1 = time.time() print("Pre-process Start!!!") +# obj.zeroMask = 1 obj.preprocess_filt_hps() +# obj.I1 = np.abs(I1) +# obj.I2 = np.abs(I2) print("Pre-process Done!!!") print(time.time()-t1) @@ -374,11 +384,13 @@ def main(): def generateAutoriftProduct(indir_m, indir_s, grid_location, init_offset, search_range, chip_size_min, chip_size_max, - offset2vx, offset2vy, stable_surface_mask, optical_flag, nc_sensor, mpflag): + offset2vx, offset2vy, stable_surface_mask, optical_flag, nc_sensor, mpflag, + geogrid_run_info=None): import numpy as np import time + import isce from components.contrib.geo_autoRIFT.autoRIFT import __version__ as version # from autoRIFT import __version__ as version @@ -470,7 +482,8 @@ def generateAutoriftProduct(indir_m, indir_s, grid_location, init_offset, search 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, optical_flag, nodata, mpflag + data_m, data_s, xGrid, yGrid, Dx0, Dy0, SRx0, SRy0, CSMINx0, CSMINy0, CSMAXx0, CSMAXy0, + noDataMask, optical_flag, nodata, mpflag, geogrid_run_info=geogrid_run_info, ) if optical_flag == 0: @@ -482,7 +495,7 @@ def generateAutoriftProduct(indir_m, indir_s, grid_location, init_offset, search CHIPSIZEX = np.zeros(origSize,dtype=np.float32) SEARCHLIMITX = np.zeros(origSize,dtype=np.float32) SEARCHLIMITY = np.zeros(origSize,dtype=np.float32) - + DX[0:Dx.shape[0],0:Dx.shape[1]] = Dx DY[0:Dy.shape[0],0:Dy.shape[1]] = Dy INTERPMASK[0:InterpMask.shape[0],0:InterpMask.shape[1]] = InterpMask @@ -514,6 +527,7 @@ def generateAutoriftProduct(indir_m, indir_s, grid_location, init_offset, search # SEARCHLIMITY = conts['SearchLimitY'] # ##################### + netcdf_file = None if grid_location is not None: @@ -558,14 +572,14 @@ def generateAutoriftProduct(indir_m, indir_s, grid_location, init_offset, search 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) @@ -575,20 +589,30 @@ def generateAutoriftProduct(indir_m, indir_s, grid_location, init_offset, search outband = outRaster.GetRasterBand(2) outband.WriteArray(VY) outband.FlushCache() - + ############ prepare for netCDF packaging if 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]) + if geogrid_run_info is 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]+"sp.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]) + else: + vxrefname = geogrid_run_info['vxname'] + vyrefname = geogrid_run_info['vyname'] + sxname = geogrid_run_info['sxname'] + syname = geogrid_run_info['syname'] + maskname = geogrid_run_info['maskname'] + xoff = geogrid_run_info['xoff'] + yoff = geogrid_run_info['yoff'] + xcount = geogrid_run_info['xcount'] + ycount = geogrid_run_info['ycount'] ds = gdal.Open(vxrefname) band = ds.GetRasterBand(1) @@ -619,23 +643,23 @@ def generateAutoriftProduct(indir_m, indir_s, grid_location, init_offset, search MM = band.ReadAsArray(xoff, yoff, xcount, ycount) ds = None band = None - + 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)) & (DX-DXref > -5) & (DX-DXref < 5) & (DY-DYref > -5) & (DY-DYref < 5)) - + if stable_count == 0: stable_shift_applied = 0 else: 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)]) @@ -643,7 +667,7 @@ def generateAutoriftProduct(indir_m, indir_s, grid_location, init_offset, search 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) @@ -652,12 +676,19 @@ def generateAutoriftProduct(indir_m, indir_s, grid_location, init_offset, search ######################################################################################## ############ netCDF packaging for Sentinel and Landsat dataset; can add other sensor format as well if 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)) + if geogrid_run_info is None: + chipsizex0 = float(str.split(runCmd('fgrep "Smallest Allowable Chip Size 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]) + epsg = float(str.split(runCmd('fgrep "EPSG:" testGeogrid.txt'))[1]) + # print (str(rangePixelSize)+" "+str(azimuthPixelSize)) + else: + chipsizex0 = geogrid_run_info['chipsizex0'] + rangePixelSize = geogrid_run_info['XPixelSize'] + azimuthPixelSize = geogrid_run_info['YPixelSize'] + dt = geogrid_run_info['dt'] + epsg = geogrid_run_info['epsg'] runCmd('topsinsar_filename.py') # import scipy.io as sio @@ -676,8 +707,8 @@ def generateAutoriftProduct(indir_m, indir_s, grid_location, init_offset, search roi_valid_percentage = int(round(np.sum(CHIPSIZEX!=0)/np.sum(SEARCHLIMITX!=0)*1000.0))/1000 # out_nc_filename = 'Jakobshavn.nc' PPP = roi_valid_percentage * 100 - out_nc_filename = master_filename[0:-4]+'_X_'+slave_filename[0:-4]+'_G0240V02_P{0}{1}{2}.nc'.format(int(PPP/10),int((PPP-int(PPP/10)*10)),int((PPP-int(PPP/10)*10-int((PPP-int(PPP/10)*10)))*10)) - out_nc_filename = './' + out_nc_filename + 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" CHIPSIZEY = np.round(CHIPSIZEX * ScaleChipSizeY / 2) * 2 @@ -693,19 +724,30 @@ def generateAutoriftProduct(indir_m, indir_s, grid_location, init_offset, search 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) + netcdf_file = 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 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]) + if geogrid_run_info is None: + chipsizex0 = float(str.split(runCmd('fgrep "Smallest Allowable Chip Size 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'] + XPixelSize = geogrid_run_info['XPixelSize'] + YPixelSize = geogrid_run_info['YPixelSize'] + epsg = geogrid_run_info['epsg'] master_path = indir_m slave_path = indir_s @@ -713,10 +755,10 @@ def generateAutoriftProduct(indir_m, indir_s, grid_location, init_offset, search 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' # @@ -724,11 +766,11 @@ def generateAutoriftProduct(indir_m, indir_s, grid_location, init_offset, search # slave_time = str.split(str.split(runCmd('fgrep "SCENE_CENTER_TIME" '+slave_MTL_path))[2][1:-2],':') master_time = ['00','00','00'] slave_time = ['00','00','00'] - + 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 pair_type = 'optical' detection_method = 'feature' @@ -736,8 +778,9 @@ def generateAutoriftProduct(indir_m, indir_s, grid_location, init_offset, search roi_valid_percentage = int(round(np.sum(CHIPSIZEX!=0)/np.sum(SEARCHLIMITX!=0)*1000.0))/1000 # out_nc_filename = 'Jakobshavn_opt.nc' PPP = roi_valid_percentage * 100 - out_nc_filename = master_filename[0:-8]+'_X_'+slave_filename[0:-8]+'_G0240V02_P{0}{1}{2}.nc'.format(int(PPP/10),int((PPP-int(PPP/10)*10)),int((PPP-int(PPP/10)*10-int((PPP-int(PPP/10)*10)))*10)) - out_nc_filename = './' + out_nc_filename + out_nc_filename = f"./{master_filename[0:-7]}_X_{slave_filename[0:-7]}" \ + f"_G{chipsizex0:04.0f}V02_P{np.floor(PPP):03.0f}.nc" + CHIPSIZEY = np.round(CHIPSIZEX * ScaleChipSizeY / 2) * 2 from datetime import date @@ -751,21 +794,33 @@ def generateAutoriftProduct(indir_m, indir_s, grid_location, init_offset, search 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 nc_sensor == "S2": + netcdf_file = 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 + ) - 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]) + 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]) + 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'] + XPixelSize = geogrid_run_info['XPixelSize'] + YPixelSize = geogrid_run_info['YPixelSize'] + epsg = geogrid_run_info['epsg'] master_path = indir_m slave_path = indir_s @@ -774,27 +829,27 @@ def generateAutoriftProduct(indir_m, indir_s, grid_location, init_offset, search slave_split = slave_path.split('_') import os - + master_filename = master_split[0][-3:]+'_'+master_split[2]+'_'+master_split[4][:3]+'_'+os.path.basename(master_path) slave_filename = slave_split[0][-3:]+'_'+slave_split[2]+'_'+slave_split[4][:3]+'_'+os.path.basename(slave_path) - + master_time = ['00','00','00'] slave_time = ['00','00','00'] - + 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 pair_type = 'optical' detection_method = 'feature' coordinates = 'map' roi_valid_percentage = int(round(np.sum(CHIPSIZEX!=0)/np.sum(SEARCHLIMITX!=0)*1000.0))/1000 PPP = roi_valid_percentage * 100 - out_nc_filename = master_filename[0:-8]+'_X_'+slave_filename[0:-8]+'_G0240V02_P{0}{1}{2}.nc'.format(int(PPP/10),int((PPP-int(PPP/10)*10)),int((PPP-int(PPP/10)*10-int((PPP-int(PPP/10)*10)))*10)) - out_nc_filename = './' + out_nc_filename + 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" CHIPSIZEY = np.round(CHIPSIZEX * ScaleChipSizeY / 2) * 2 - + from datetime import date d0 = date(np.int(master_split[2][0:4]),np.int(master_split[2][4:6]),np.int(master_split[2][6:8])) d1 = date(np.int(slave_split[2][0:4]),np.int(slave_split[2][4:6]),np.int(slave_split[2][6:8])) @@ -806,15 +861,21 @@ def generateAutoriftProduct(indir_m, indir_s, grid_location, init_offset, search else: date_ct = d0 + (d1 - d0)/2 date_center = date_ct.strftime("%Y%m%d") - + 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,'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) + + netcdf_file = 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 nc_sensor is None: print('netCDF packaging not performed') @@ -825,6 +886,8 @@ def generateAutoriftProduct(indir_m, indir_s, grid_location, init_offset, search print("Write Outputs Done!!!") print(time.time()-t1) + return netcdf_file + if __name__ == '__main__': main()