From f4d98f6606a201c168e273c4f8bf55e6d557de84 Mon Sep 17 00:00:00 2001 From: Yang Lei Date: Fri, 30 Oct 2020 18:22:21 -0700 Subject: [PATCH] update readers for optical data with gdal vsicurl access --- geo_autoRIFT/geogrid/GeogridOptical.py | 162 ++++++++++++++++++++++--- testautoRIFT.py | 133 +++++++++++++++++--- testautoRIFT_ISCE.py | 129 ++++++++++++++++++-- 3 files changed, 379 insertions(+), 45 deletions(-) diff --git a/geo_autoRIFT/geogrid/GeogridOptical.py b/geo_autoRIFT/geogrid/GeogridOptical.py index 1ce6f88..3559d7e 100755 --- a/geo_autoRIFT/geogrid/GeogridOptical.py +++ b/geo_autoRIFT/geogrid/GeogridOptical.py @@ -48,8 +48,8 @@ def runGeogrid(self): ''' ##Determine appropriate EPSG system - self.epsgDem = self.getProjectionSystem(self.demname) - self.epsgDat = self.getProjectionSystem(self.dat1name) + self.epsgDem = self.getProjectionSystem(self.demname, self.urlflag) + self.epsgDat = self.getProjectionSystem(self.dat1name, self.urlflag) ###Determine extent of data needed bbox = self.determineBbox() @@ -59,7 +59,7 @@ def runGeogrid(self): self.geogrid() - def getProjectionSystem(self, filename): + def getProjectionSystem(self, filename, urlflag): ''' Testing with Greenland. ''' @@ -67,7 +67,10 @@ def getProjectionSystem(self, filename): raise Exception('File {0} does not exist'.format(filename)) from osgeo import gdal, osr - ds = gdal.Open(filename, gdal.GA_ReadOnly) + if urlflag is 1: + ds = gdal.Open('/vsicurl/%s' %(filename)) + else: + ds = gdal.Open(filename, gdal.GA_ReadOnly) srs = osr.SpatialReference() srs.ImportFromWkt(ds.GetProjection()) srs.AutoIdentifyEPSG() @@ -83,7 +86,10 @@ def getProjectionSystem(self, filename): else: raise Exception('Non-standard coordinate system encountered') if not epsgstr: #Empty string->use shell command gdalsrsinfo for last trial - cmd = 'gdalsrsinfo -o epsg {0}'.format(filename) + if urlflag is 1: + cmd = 'gdalsrsinfo -o epsg /vsicurl/{0}'.format(filename) + else: + cmd = 'gdalsrsinfo -o epsg {0}'.format(filename) epsgstr = subprocess.check_output(cmd, shell=True) # pdb.set_trace() epsgstr = re.findall("EPSG:(\d+)", str(epsgstr))[0] @@ -154,6 +160,13 @@ def determineBbox(self, zrange=[-200,4000]): def geogrid(self): # For now print inputs that were obtained + + urlflag = self.urlflag + + if urlflag is 1: + print("\nReading input images into memory directly from URL's") + else: + print("\nReading input images locally from files") print("\nOptical Image parameters: ") print("X-direction coordinate: " + str(self.startingX) + " " + str(self.XSize)) @@ -218,30 +231,56 @@ def geogrid(self): import struct # pdb.set_trace() - demDS = gdal.Open(self.demname, gdal.GA_ReadOnly) + if urlflag is 1: + demDS = gdal.Open('/vsicurl/%s' %(self.demname)) + else: + 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 urlflag is 1: + sxDS = gdal.Open('/vsicurl/%s' %(self.dhdxname)) + syDS = gdal.Open('/vsicurl/%s' %(self.dhdyname)) + else: + 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 urlflag is 1: + vxDS = gdal.Open('/vsicurl/%s' %(self.vxname)) + vyDS = gdal.Open('/vsicurl/%s' %(self.vyname)) + else: + 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 urlflag is 1: + srxDS = gdal.Open('/vsicurl/%s' %(self.srxname)) + sryDS = gdal.Open('/vsicurl/%s' %(self.sryname)) + else: + 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 urlflag is 1: + csminxDS = gdal.Open('/vsicurl/%s' %(self.csminxname)) + csminyDS = gdal.Open('/vsicurl/%s' %(self.csminyname)) + else: + 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 urlflag is 1: + csmaxxDS = gdal.Open('/vsicurl/%s' %(self.csmaxxname)) + csmaxyDS = gdal.Open('/vsicurl/%s' %(self.csmaxyname)) + else: + 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 urlflag is 1: + ssmDS = gdal.Open('/vsicurl/%s' %(self.ssmname)) + else: + ssmDS = gdal.Open(self.ssmname, gdal.GA_ReadOnly) if demDS is None: raise Exception('Error opening DEM file {0}'.format(self.demname)) @@ -731,8 +770,94 @@ def geogrid(self): + def coregister(self,in1,in2,urlflag): + import os + import numpy as np + + from osgeo import gdal, osr + import struct + + if urlflag is 1: + DS1 = gdal.Open('/vsicurl/%s' %(in1)) + else: + DS1 = gdal.Open(in1, gdal.GA_ReadOnly) + trans1 = DS1.GetGeoTransform() + xsize1 = DS1.RasterXSize + ysize1 = DS1.RasterYSize + + if urlflag is 1: + DS2 = gdal.Open('/vsicurl/%s' %(in2)) + else: + DS2 = gdal.Open(in2, gdal.GA_ReadOnly) + trans2 = DS2.GetGeoTransform() + xsize2 = DS2.RasterXSize + ysize2 = DS2.RasterYSize + + W = np.max([trans1[0],trans2[0]]) + N = np.min([trans1[3],trans2[3]]) + E = np.min([trans1[0]+xsize1*trans1[1],trans2[0]+xsize2*trans2[1]]) + S = np.max([trans1[3]+ysize1*trans1[5],trans2[3]+ysize2*trans2[5]]) + + x1a = int(np.round((W-trans1[0])/trans1[1])) + x1b = int(np.round((E-trans1[0])/trans1[1])) + y1a = int(np.round((N-trans1[3])/trans1[5])) + y1b = int(np.round((S-trans1[3])/trans1[5])) + + x2a = int(np.round((W-trans2[0])/trans2[1])) + x2b = int(np.round((E-trans2[0])/trans2[1])) + y2a = int(np.round((N-trans2[3])/trans2[5])) + y2b = int(np.round((S-trans2[3])/trans2[5])) + + x1a = np.min([x1a, xsize1-1]) + x1b = np.min([x1b, xsize1-1]) + y1a = np.min([y1a, ysize1-1]) + y1b = np.min([y1b, ysize1-1]) + x2a = np.min([x2a, xsize2-1]) + x2b = np.min([x2b, xsize2-1]) + y2a = np.min([y2a, ysize2-1]) + y2b = np.min([y2b, ysize2-1]) + + x1a = np.max([x1a, 0]) + x1b = np.max([x1b, 0]) + y1a = np.max([y1a, 0]) + y1b = np.max([y1b, 0]) + x2a = np.max([x2a, 0]) + x2b = np.max([x2b, 0]) + y2a = np.max([y2a, 0]) + y2b = np.max([y2b, 0]) + + x1a = int(x1a) + x1b = int(x1b) + y1a = int(y1a) + y1b = int(y1b) + x2a = int(x2a) + x2b = int(x2b) + y2a = int(y2a) + y2b = int(y2b) + + trans = (W, trans1[1], 0.0, N, 0.0, trans1[5]) + + if urlflag is 0: + + I1 = DS1.ReadAsArray(xoff=x1a, yoff=y1a, xsize=x1b-x1a+1, ysize=y1b-y1a+1) + I2 = DS2.ReadAsArray(xoff=x2a, yoff=y2a, xsize=x2b-x2a+1, ysize=y2b-y2a+1) - + fileformat = "GTiff" + driver = gdal.GetDriverByName(fileformat) + + DST1 = driver.Create(os.path.basename(in1), xsize=(x1b-x1a+1), ysize=(y1b-y1a+1), bands=1, eType=gdal.GDT_UInt16) + DST1.SetGeoTransform(trans) + DST1.SetProjection(DS1.GetProjectionRef()) + DST1.GetRasterBand(1).WriteArray(I1) + DST1 = None + + DST2 = driver.Create(os.path.basename(in2), xsize=(x2b-x2a+1), ysize=(y2b-y2a+1), bands=1, eType=gdal.GDT_UInt16) + DST2.SetGeoTransform(trans) + DST2.SetProjection(DS2.GetProjectionRef()) + DST2.GetRasterBand(1).WriteArray(I2) + DST2 = None + + return x1a, y1a, x1b-x1a+1, y1b-y1a+1, x2a, y2a, x2b-x2a+1, y2b-y2a+1, trans @@ -753,6 +878,7 @@ def __init__(self): self.numberOfLines = None self.repeatTime = None self.chipSizeX0 = None + self.urlflag = None ##Input related parameters self.dat1name = None diff --git a/testautoRIFT.py b/testautoRIFT.py index 273d699..97c77fb 100755 --- a/testautoRIFT.py +++ b/testautoRIFT.py @@ -75,7 +75,8 @@ def cmdLineParse(): 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, help='flag for packaging output formatted for Sentinel ("S") and Landsat ("L") dataset; default is None') - + parser.add_argument('-urlflag', '--urlflag', dest='urlflag', type=int, required=False, default=0, + help='flag for reading and coregistering optical data (GeoTIFF images, e.g. Landsat): use 1 for url read and 0 for local machine read; if not specified (i.e. None; default), will just read from local machine without coregistration') return parser.parse_args() @@ -116,6 +117,32 @@ def loadProductOptical(filename): return img +def loadProductOpticalURL(file_m, file_s): + import numpy as np + ''' + 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, 1) + + DS1 = gdal.Open('/vsicurl/%s' %(file_m)) + DS2 = gdal.Open('/vsicurl/%s' %(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): @@ -368,9 +395,14 @@ def runAutorift(I1, I2, xGrid, yGrid, Dx0, Dy0, SRx0, SRy0, CSMINx0, CSMINy0, CS inps = cmdLineParse() + urlflag = inps.urlflag + if inps.optical_flag == 1: - data_m = loadProductOptical(inps.indir_m) - data_s = loadProductOptical(inps.indir_s) + if urlflag is 1: + data_m, data_s = loadProductOpticalURL(inps.indir_m, inps.indir_s) + else: + 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) @@ -573,31 +605,46 @@ def runAutorift(I1, I2, xGrid, yGrid, Dx0, Dy0, SRx0, SRy0, CSMINx0, CSMINy0, CS 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) + if urlflag is 1: + ds = gdal.Open('/vsicurl/%s' %(vxrefname)) + else: + ds = gdal.Open(vxrefname) band = ds.GetRasterBand(1) VXref = band.ReadAsArray(xoff, yoff, xcount, ycount) ds = None band = None - ds = gdal.Open(vyrefname) + if urlflag is 1: + ds = gdal.Open('/vsicurl/%s' %(vyrefname)) + else: + ds = gdal.Open(vyrefname) band = ds.GetRasterBand(1) VYref = band.ReadAsArray(xoff, yoff, xcount, ycount) ds = None band = None - ds = gdal.Open(sxname) + if urlflag is 1: + ds = gdal.Open('/vsicurl/%s' %(sxname)) + else: + ds = gdal.Open(sxname) band = ds.GetRasterBand(1) SX = band.ReadAsArray(xoff, yoff, xcount, ycount) ds = None band = None - ds = gdal.Open(syname) + if urlflag is 1: + ds = gdal.Open('/vsicurl/%s' %(syname)) + else: + ds = gdal.Open(syname) band = ds.GetRasterBand(1) SY = band.ReadAsArray(xoff, yoff, xcount, ycount) ds = None band = None - ds = gdal.Open(maskname) + if urlflag is 1: + ds = gdal.Open('/vsicurl/%s' %(maskname)) + else: + ds = gdal.Open(maskname) band = ds.GetRasterBand(1) MM = band.ReadAsArray(xoff, yoff, xcount, ycount) ds = None @@ -699,12 +746,14 @@ def runAutorift(I1, I2, xGrid, yGrid, Dx0, Dy0, SRx0, SRy0, CSMINx0, CSMINy0, CS 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],':') - +# 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],':') + 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]))) @@ -740,6 +789,62 @@ def runAutorift(I1, I2, xGrid, yGrid, Dx0, Dy0, SRx0, SRy0, CSMINx0, CSMINy0, CS 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 == "S2": + + 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 + + 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 + version = '1.0.7' + pair_type = 'optical' + detection_method = 'feature' + coordinates = 'map' + + 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[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])) + 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[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) elif inps.nc_sensor is None: print('netCDF packaging not performed') diff --git a/testautoRIFT_ISCE.py b/testautoRIFT_ISCE.py index 59c3080..ba79211 100755 --- a/testautoRIFT_ISCE.py +++ b/testautoRIFT_ISCE.py @@ -75,7 +75,8 @@ def cmdLineParse(): 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, help='flag for packaging output formatted for Sentinel ("S") and Landsat ("L") dataset; default is None') - + parser.add_argument('-urlflag', '--urlflag', dest='urlflag', type=int, required=False, default=0, + help='flag for reading and coregistering optical data (GeoTIFF images, e.g. Landsat): use 1 for url read and 0 for local machine read; if not specified (i.e. None; default), will just read from local machine without coregistration') return parser.parse_args() @@ -115,7 +116,31 @@ def loadProductOptical(filename): return img +def loadProductOpticalURL(file_m, file_s): + import numpy as np + ''' + Load the product using Product Manager. + ''' + 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, 1) + + DS1 = gdal.Open('/vsicurl/%s' %(file_m)) + DS2 = gdal.Open('/vsicurl/%s' %(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): @@ -363,9 +388,14 @@ def runAutorift(I1, I2, xGrid, yGrid, Dx0, Dy0, SRx0, SRy0, CSMINx0, CSMINy0, CS inps = cmdLineParse() + urlflag = inps.urlflag + if inps.optical_flag == 1: - data_m = loadProductOptical(inps.indir_m) - data_s = loadProductOptical(inps.indir_s) + if urlflag is 1: + data_m, data_s = loadProductOpticalURL(inps.indir_m, inps.indir_s) + else: + 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) @@ -570,31 +600,46 @@ def runAutorift(I1, I2, xGrid, yGrid, Dx0, Dy0, SRx0, SRy0, CSMINx0, CSMINy0, CS 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) + if urlflag is 1: + ds = gdal.Open('/vsicurl/%s' %(vxrefname)) + else: + ds = gdal.Open(vxrefname) band = ds.GetRasterBand(1) VXref = band.ReadAsArray(xoff, yoff, xcount, ycount) ds = None band = None - ds = gdal.Open(vyrefname) + if urlflag is 1: + ds = gdal.Open('/vsicurl/%s' %(vyrefname)) + else: + ds = gdal.Open(vyrefname) band = ds.GetRasterBand(1) VYref = band.ReadAsArray(xoff, yoff, xcount, ycount) ds = None band = None - ds = gdal.Open(sxname) + if urlflag is 1: + ds = gdal.Open('/vsicurl/%s' %(sxname)) + else: + ds = gdal.Open(sxname) band = ds.GetRasterBand(1) SX = band.ReadAsArray(xoff, yoff, xcount, ycount) ds = None band = None - ds = gdal.Open(syname) + if urlflag is 1: + ds = gdal.Open('/vsicurl/%s' %(syname)) + else: + ds = gdal.Open(syname) band = ds.GetRasterBand(1) SY = band.ReadAsArray(xoff, yoff, xcount, ycount) ds = None band = None - ds = gdal.Open(maskname) + if urlflag is 1: + ds = gdal.Open('/vsicurl/%s' %(maskname)) + else: + ds = gdal.Open(maskname) band = ds.GetRasterBand(1) MM = band.ReadAsArray(xoff, yoff, xcount, ycount) ds = None @@ -697,11 +742,13 @@ def runAutorift(I1, I2, xGrid, yGrid, Dx0, Dy0, SRx0, SRy0, CSMINx0, CSMINy0, CS 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],':') +# 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],':') + 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]))) @@ -738,6 +785,62 @@ def runAutorift(I1, I2, xGrid, yGrid, Dx0, Dy0, SRx0, SRy0, CSMINx0, CSMINy0, CS 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 == "S2": + + 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 + + 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 + version = '1.0.7' + pair_type = 'optical' + detection_method = 'feature' + coordinates = 'map' + + 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[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])) + 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[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) elif inps.nc_sensor is None: print('netCDF packaging not performed')