Commit cd1e1163 authored by Gaetano Raffaele's avatar Gaetano Raffaele

merging with branch develop

parents 806b2eb0 f5498fe3
......@@ -6,6 +6,16 @@ BASEFOLDER: /home/raffaele/Data/Koumbia/2016
EXTENT: /home/raffaele/Data/Koumbia/Extent/koumbia_buf.shp
REFERENCE: /home/raffaele/Data/Koumbia/GT/GT_FINAL.shp
CHAINNAME: OUT
SETUPNAME: S2_L8
#Chain mode: 0 - Validation only, 1 - Map only, 2 - Both
CHAINMODE: 1
VALIDATION: /home/raffaele/Data/Koumbia/GT/GT_VALID.shp
#Validation mode: 0 - Pixel, 1 - Vector (surface)
VALIDMODE: 1
#Number of parallel processes (only affects zonal statistics)
NPROCESSES: 4
[SENTINEL-2 CONFIGURATION]
......@@ -115,18 +125,30 @@ SCALE: 250
SHAPE: 0.7
COMPACTNESS: 0.9
MEMLIMIT: OFF
<<<<<<< HEAD:chain_setup_lnx.cfg
=======
#Area threshold for GT sample generation
GTAREATHRESH: 225
>>>>>>> develop:chain_setup.cfg
[VHR-SITS COREGISTRATION]
COREGISTER: YES
FORCE: NO
#Select band to be used for homologous points extraction
VHRBAND: 1
S2BAND: 3
L8BAND: 4
<<<<<<< HEAD:chain_setup_lnx.cfg
#Area threshold for GT sample generation
GTAREATHRESH: 225
=======
#Maximum number of GCPs (-1 unlimited)
MAXPOINTS: -1
>>>>>>> develop:chain_setup.cfg
[TRAINING CONFIGURATION]
......
......@@ -3,6 +3,9 @@ import sys
import subprocess
import platform
import numpy as np
import os
from mtdUtils import queuedProcess
def getFeaturesFields(shp,flds_pref):
......@@ -54,7 +57,7 @@ def roughFix(shp,flds):
return
def training(shp,feat_prefix,code,model_fld,params):
def training(shp,code,model_fld,params,feat,feat_mode = 'list'):
# Platform dependent parameters
if platform.system() == 'Linux':
......@@ -72,9 +75,12 @@ def training(shp,feat_prefix,code,model_fld,params):
confmat_file = model_fld + '/' + classifier + '_' + code + '.confmat.txt'
stat_file = model_fld + '/GT_stats.xml'
flds = getFeaturesFields(shp,feat_prefix)
#Rough Fix made in the zonal statistics step
#roughFix(shp,flds)
if feat_mode == 'prefix':
flds = getFeaturesFields(shp, feat)
elif feat_mode == 'list':
flds = feat
else:
sys.exit('ERROR: mode ' + feat_mode + ' not valid.')
if platform.system() == 'Linux':
cmd = ['otbcli_ComputeVectorFeaturesStatistics','-io.vd',shp,'-io.stats',stat_file,'-feat'] + flds
......@@ -117,7 +123,7 @@ def training(shp,feat_prefix,code,model_fld,params):
return flds
def classify(shp_list,feat_prefix,code,stat_file,model_file):
def classify(shp_list,code,stat_file,model_file,out_fld,out_ext,feat,feat_mode = 'list',Nproc=1):
# Platform dependent parameters
if platform.system() == 'Linux':
......@@ -127,11 +133,18 @@ def classify(shp_list,feat_prefix,code,stat_file,model_file):
else:
sys.exit("Platform not supported!")
flds = getFeaturesFields(shp_list[0], feat_prefix)
if feat_mode == 'prefix':
flds = getFeaturesFields(shp_list[0], feat)
elif feat_mode == 'list':
flds = feat
else:
sys.exit('ERROR: mode ' + feat_mode + ' not valid.')
'''
for shp in shp_list:
#roughFix(shp, flds)
if platform.system() == 'Linux':
cmd = ['otbcli_VectorClassifier','-in',shp,'-instat',stat_file,'-model',model_file,'-cfield',code,'-feat'] + flds
cmd = ['otbcli_VectorClassifier','-in',shp,'-instat',stat_file,'-model',model_file,'-out',out_file,'-cfield',code,'-feat'] + flds
subprocess.call(cmd,shell=sh)
elif platform.system() == 'Windows':
import otbApplication
......@@ -139,6 +152,7 @@ def classify(shp_list,feat_prefix,code,stat_file,model_file):
app.SetParameterString('in',shp)
app.SetParameterString('instat', stat_file)
app.SetParameterString('model', model_file)
app.SetParameterString('out', out_file)
app.SetParameterString('cfield', code)
app.UpdateParameters()
app.SetParameterStringList('feat',flds)
......@@ -148,3 +162,36 @@ def classify(shp_list,feat_prefix,code,stat_file,model_file):
sys.exit('Platform not supported!')
return
'''
if platform.system() == 'Linux':
cmd_list = []
out_file_list = []
for shp in shp_list:
out_file = out_fld + '/' + os.path.basename(shp).replace('.shp', out_ext + '.shp')
cmd = ['otbcli_VectorClassifier', '-in', shp, '-instat', stat_file, '-model', model_file, '-out', out_file,
'-cfield', code, '-feat'] + flds
cmd_list.append(cmd)
out_file_list.append(out_file)
queuedProcess(cmd_list,Nproc,shell=sh)
return out_file_list
elif platform.system() == 'Windows':
out_file_list = []
for shp in shp_list:
out_file = out_fld + '/' + os.path.basename(shp).replace('.shp', out_ext + '.shp')
import otbApplication
app = otbApplication.Registry.CreateApplication('VectorClassifier')
app.SetParameterString('in', shp)
app.SetParameterString('instat', stat_file)
app.SetParameterString('model', model_file)
app.SetParameterString('out', out_file)
app.SetParameterString('cfield', code)
app.UpdateParameters()
app.SetParameterStringList('feat', flds)
app.UpdateParameters()
app.ExecuteAndWriteOutput()
out_file_list.append(out_file)
return out_file_list
else:
sys.exit('Platform not supported!')
import getopt
import glob
import os
import platform
import subprocess
import sys
import gdal
from mtdUtils import setNoDataValue
gdt_to_otb = {2 : 'uint16', 3 : 'int16', 6 : 'float'}
def gatherDates(fld,ofld='.'):
lst = sorted([name for name in os.listdir(fld) if (os.path.isdir(os.path.join(fld, name)) and (name[0:2] in ['S2','L8']))])
ofn = ofld + '/input_dates.txt'
sns = set()
with open(ofn,'w') as fl:
for l in lst:
fl.write(l.split('_')[-2] + '\n')
sns.add(l[0:2])
if len(sns) != 1 or list(sns)[0] not in ['S2','L8']:
sys.exit('Error: directory contains images from multiple sensors or invalid detected folders.')
return lst,list(sns)
def genComposite(fld, od_fn, ofld = '.'):
if platform.system() == 'Linux':
sh = False
elif platform.system() == 'Windows':
sh = True
else:
sys.exit("Platform not supported!")
lst,sns = gatherDates(fld,ofld)
ptrn = '*STACK_*.tif'
if sns == 'L8':
ptrn = '*_PS_TOA_*.tif'
gptrn = '*GAPMASK*.tif'
stacklist = []
gaplist = []
nc = set()
dt = set()
for l in lst:
fn = glob.glob(fld + '/' + l + '/' + ptrn)[0]
ds = gdal.Open(fn)
nc.add(ds.RasterCount)
dt.add(ds.GetRasterBand(1).DataType)
ds = None
stacklist.append(fn)
fn = glob.glob(fld + '/' + l + '/' + gptrn)[0]
gaplist.append(fn)
if len(nc) != 1:
sys.exit('Error: stacks to composite do not have the same number of bands.')
if len(dt) != 1:
sys.exit('Error: stacks to composite do not have the same data format.')
dtout = gdt_to_otb[list(dt)[0]]
ncomp = str(list(nc)[0])
cmd = ['otbcli_ConcatenateImages','-il'] + stacklist + ['-out', ofld + '/TIMESERIES.tif', dtout]
subprocess.call(cmd,shell=sh)
setNoDataValue(ofld + '/TIMESERIES.tif', -10000)
cmd = ['otbcli_ConcatenateImages', '-il'] + gaplist + ['-out', ofld + '/GAPMASKS.tif', 'uint8']
subprocess.call(cmd, shell=sh)
cmd = ['otbcli_ImageTimeSeriesGapFilling', '-in', ofld + '/TIMESERIES.tif', '-mask', ofld + '/GAPMASKS.tif', '-out', ofld + '/TIMESERIES_GAPF.tif', dtout, '-comp', ncomp, '-it', 'linear', '-id', ofld + '/input_dates.txt', '-od', od_fn]
subprocess.call(cmd, shell=sh)
setNoDataValue(ofld + '/TIMESERIES_GAPF.tif', -10000)
return
\ No newline at end of file
......@@ -26,65 +26,75 @@ def readConfigFile(cfg_fn, section):
def L8_readTile(fld):
mdf = glob.glob(fld + '/*_MTL*.txt')[0]
out = 'N.A.'
with open(mdf) as mtl:
for line in mtl:
val = line.strip().split(' = ')
if 'LANDSAT_SCENE_ID' in val[0]:
out = val[1][4:10]
return out
loc_id = os.path.basename(os.path.split(mdf)[0])
if 'MOSAIC' in loc_id:
return 'MOSAIC'
else:
out = 'N.A.'
with open(mdf) as mtl:
for line in mtl:
val = line.strip().split(' = ')
if 'LANDSAT_SCENE_ID' in val[0]:
out = val[1][4:10]
return out
def L8_readDate(fld):
mdf = glob.glob(fld + '/*_MTL*.txt')[0]
out = 'N.A.'
with open(mdf) as mtl:
for line in mtl:
val = line.strip().split(' = ')
if 'FILE_DATE' in val[0]:
out = val[1][0:10].replace('-', '')
return out
loc_id = os.path.basename(os.path.split(mdf)[0])
if 'MOSAIC' in loc_id:
return loc_id[-15:-7]
else:
out = 'N.A.'
with open(mdf) as mtl:
for line in mtl:
val = line.strip().split(' = ')
if 'FILE_DATE' in val[0]:
out = val[1][0:10].replace('-', '')
return out
def S2_readDate(fld):
tli = glob.glob(fld + '/tileInfo.json')[0]
with open(tli) as json_data:
data = json.load(json_data)
out = data['timestamp'][0:10].replace('-', '')
return out
loc_id = os.path.basename(os.path.split(tli)[0])
if 'MOSAIC' in loc_id:
return loc_id[-15:-7]
else:
with open(tli) as json_data:
data = json.load(json_data)
out = data['timestamp'][0:10].replace('-', '')
return out
def S2_readTile(fld):
tli = glob.glob(fld + '/tileInfo.json')[0]
tlidir = os.path.basename(os.path.split(tli)[0])
if 'MOSAIC' in tlidir:
loc_id = os.path.basename(os.path.split(tli)[0])
if 'MOSAIC' in loc_id:
out = 'MOSAIC'
else:
with open(tli) as json_data:
data = json.load(json_data)
out = str(data['utmZone']) + str(data['latitudeBand']) + str(data['gridSquare'])
return out
def S2THEIA_readDate(fld):
mtd = glob.glob(fld + '/*_MTD_ALL.xml')[0]
dte = ET.parse(mtd).getroot().iter('ACQUISITION_DATE').next().text
out = dte[0:4] + dte[5:7] + dte[8:10]
return out
loc_id = os.path.basename(os.path.split(mtd)[0])
if 'MOSAIC' in loc_id:
return loc_id[-15:-7]
else:
dte = ET.parse(mtd).getroot().iter('ACQUISITION_DATE').next().text
out = dte[0:4] + dte[5:7] + dte[8:10]
return out
def S2THEIA_readTile(fld):
mtd = glob.glob(fld + '/*_MTD_ALL.xml')[0]
out = ET.parse(mtd).getroot().iter('GEOGRAPHICAL_ZONE').next().text[1:]
loc_id = os.path.basename(os.path.split(mtd)[0])
if 'MOSAIC' in loc_id:
out = 'MOSAIC'
else:
out = ET.parse(mtd).getroot().iter('GEOGRAPHICAL_ZONE').next().text[1:]
return out
......@@ -306,7 +316,11 @@ def featureComputation(fld, cfg, sensor, out_fld = None, extent = None):
expr = []
feat = cfg['feat'].split(',')
if cfg['feat'] == '':
return
else:
feat = cfg['feat'].split(',')
if cfg['gapfillingmode'] == 'BASE':
gapf = 1 #Not implemented
elif cfg['gapfillingmode'] == 'FEATURE':
......
......@@ -10,7 +10,7 @@ from shutil import rmtree, copyfile
import gdal
import osr
from mtdUtils import homoGeo2Pixel, randomword
from mtdUtils import homoGeo2Pixel, randomword, compareImageGeometries
def main(argv):
......@@ -22,7 +22,7 @@ def main(argv):
sys.exit("Platform not supported!")
try:
opts, args = getopt.getopt(argv, '', ['band-vhr=', 'band-hr=', 'step=', 'minstep=', 'minpoints=', 'prec='])
opts, args = getopt.getopt(argv, '', ['band-ref=', 'band-in=', 'step=', 'minstep=', 'minpoints=', 'maxpoints=', 'prec='])
except getopt.GetoptError as err:
print str(err)
......@@ -43,13 +43,14 @@ def main(argv):
# Internal calibration.. To adapt?
step = 256
minstep = 16
enough = 40
enough = 20
prec = 3
maxgcp = -1
for opt, val in opts:
if opt == '--band-vhr':
if opt == '--band-ref':
vhrb = val
elif opt == '--band-hr':
elif opt == '--band-in':
hrb = val
elif opt == '--step':
step = int(val)
......@@ -57,8 +58,10 @@ def main(argv):
minstep = int(val)
elif opt == '--minpoints':
enough = int(val)
elif opt == '--maxpoints':
maxgcp = int(val)
elif opt == '--prec':
enough = float(val)
prec = float(val)
tmphomo = randomword(16) + '.geom'
tmpfld = os.path.split(args[1])[0] + '/' + randomword(16) + '/'
......@@ -68,15 +71,19 @@ def main(argv):
sys.exit('Oups, temp directory exists!')
ext = os.path.splitext(args[1])[1]
cmd = ['otbcli_Superimpose', '-inr', args[1], '-inm', args[0], '-out', tmpfld + 'ref.tif', 'uint16']
subprocess.call(cmd, shell=sh)
if not compareImageGeometries(args[0],args[1]):
refimg = tmpfld + 'ref.tif'
cmd = ['otbcli_Superimpose', '-inr', args[1], '-inm', args[0], '-out', refimg, 'uint16']
subprocess.call(cmd, shell=sh)
else:
refimg = args[0]
# Internal calibration.. To adapt?
npts = 0
while npts < enough and step >= minstep:
cmd = ['otbcli_HomologousPointsExtraction', '-in1', args[1], '-band1', str(hrb), '-in2', tmpfld + 'ref.tif',
cmd = ['otbcli_HomologousPointsExtraction', '-in1', args[1], '-band1', str(hrb), '-in2', refimg,
'-band2',
str(vhrb), '-mode', 'geobins', '-precision', str(prec), '-mfilter', '1', '-backmatching', '1', '-out',
tmpfld + 'homo.txt',
......@@ -92,17 +99,32 @@ def main(argv):
if npts >= enough:
homoGeo2Pixel(tmpfld + 'homo.txt', args[1], tmpfld + 'homocorr.txt')
homo_sample_rate = None
if maxgcp > 0:
homo_sample_rate = float(maxgcp)/float(npts)
homoGeo2Pixel(tmpfld + 'homo.txt', args[1], tmpfld + 'homocorr.txt', sample_rate=homo_sample_rate)
cmd = ['otbcli_GenerateRPCSensorModel', '-inpoints', tmpfld + 'homocorr.txt', '-outgeom', tmphomo]
subprocess.call(cmd, shell=sh)
cmd = ['otbcli_OrthoRectification', '-io.in \"' + args[1] + '?&skipcarto=true&geom=' + tmphomo + '\"',
'-io.out',
tmpfld + 'sits_band_coreg.tif', 'uint16', '-map', 'epsg', '-map.epsg.code', code, '-opt.gridspacing',
gsp]
# Remove Cartographic Information from source (temporary)
ds = gdal.Open(args[1],1)
geoT = ds.GetGeoTransform()
proj = ds.GetProjection()
ds.SetGeoTransform(tuple([0.0,1.0,0.0,0.0,0.0,1.0]))
ds.SetProjection('')
ds = None
cmd = ['otbcli_OrthoRectification', '-io.in \"' + args[1] + '?&geom=' + tmphomo + '\"', '-io.out', tmpfld + 'sits_band_coreg.tif', 'uint16', '-map', 'epsg', '-map.epsg.code', code, '-opt.gridspacing', gsp]
cmdex = ' '.join(cmd)
subprocess.call(cmdex, shell=True)
# Restore Cartographic Information on source
ds = gdal.Open(args[1], 1)
ds.SetGeoTransform(geoT)
ds.SetProjection(proj)
ds = None
cmd = ['otbcli_Superimpose', '-inr', args[1], '-inm', tmpfld + 'sits_band_coreg.tif', '-out',
args[1].replace(ext, ext.replace('.', '_COREG.')), 'uint16']
subprocess.call(cmd, shell=sh)
......@@ -112,13 +134,25 @@ def main(argv):
ff = os.path.split(f)[-1]
fn = os.path.splitext(ff)
if fn[1] in ['.tif', '.TIF']:
cmd = ['otbcli_OrthoRectification', '-io.in \"' + f + '?&skipcarto=true&geom=' + tmphomo + '\"',
'-io.out',
tmpfld + fn[0] + '_coreg.tif', 'uint8', '-map', 'epsg', '-map.epsg.code', code,
'-opt.gridspacing',
gsp, '-interpolator nn']
# Remove Cartographic Information from source (temporary)
ds = gdal.Open(f, 1)
geoT = ds.GetGeoTransform()
proj = ds.GetProjection()
ds.SetGeoTransform(tuple([0.0, 1.0, 0.0, 0.0, 0.0, 1.0]))
ds.SetProjection('')
ds = None
cmd = ['otbcli_OrthoRectification', '-io.in \"' + f + '?&geom=' + tmphomo + '\"', '-io.out', tmpfld + fn[0] + '_coreg.tif', 'uint8', '-map', 'epsg', '-map.epsg.code', code, '-opt.gridspacing', gsp, '-interpolator nn']
cmdex = ' '.join(cmd)
subprocess.call(cmdex, shell=True)
# Restore Cartographic Information on source
ds = gdal.Open(f, 1)
ds.SetGeoTransform(geoT)
ds.SetProjection(proj)
ds = None
cmd = ['otbcli_Superimpose', '-inr', args[1], '-inm', tmpfld + fn[0] + '_coreg.tif', '-out',
f.replace(fn[1], fn[1].replace('.', '_COREG.')), 'uint8', '-interpolator', 'nn']
subprocess.call(cmd, shell=sh)
......@@ -136,6 +170,6 @@ def main(argv):
if __name__ == '__main__':
if len(sys.argv) < 2:
sys.exit(
'Usage: python coregister.py [--band-vhr <target vhr band>] [--band-hr <source hr band>] [--step <geobins initial step def. 256>] [--minstep <geobins minimum step def. 16>] [--minpoints <minimum # of homologous points def. 40>] [--prec <colocalisation precision in pixels] <vhsr scene file> <stack to coregister>')
'Usage: python coregister.py [--band-ref <target ref band>] [--band-in <source band>] [--step <geobins initial step def. 256>] [--minstep <geobins minimum step def. 16>] [--minpoints <minimum # of homologous points def. 40>] [--maxpoints <maximum # of GCP points def. -1 unlimited> [--prec <colocalisation precision in pixels] <vhsr scene file> <stack to coregister>')
else:
main(sys.argv[1:])
......@@ -108,8 +108,8 @@ def main(argv):
for im in frdr:
dt = im[2][0:4]+im[2][5:7]+im[2][8:10]
pathrow = (str(int(im[5])).zfill(3) + str(int(im[6])).zfill(3))
if float(im[3]) <= maxcc and dt not in found_dates:
found_dates.append(dt)
if float(im[3]) <= maxcc and pathrow+dt not in found_dates:
found_dates.append(pathrow+dt)
f.write('python getLandsat8.py' + opt_str + ' ' + '--url ' + os.path.dirname(im[-1]) + '/ ' + dt + ' ' + pathrow + '\n')
g.write('Adding date ' + dt + ' for path/row ' + pathrow[0:3] + '/' + pathrow[3:6] + ' with ' + str(float(im[3])) + '% cloudy pixels (Collection 1)\n')
else:
......@@ -129,8 +129,8 @@ def main(argv):
for im in frdr:
dt = im[1][0:4] + im[1][5:7] + im[1][8:10]
pathrow = (str(int(im[4])).zfill(3) + str(int(im[5])).zfill(3))
if float(im[2]) <= maxcc and dt not in found_dates:
found_dates.append(dt)
if float(im[2]) <= maxcc and pathrow+dt not in found_dates:
found_dates.append(pathrow+dt)
f.write('python getLandsat8.py' + opt_str + ' ' + '--url ' + os.path.dirname(
im[-1]) + '/ ' + dt + ' ' + pathrow + '\n')
g.write('Adding date ' + dt + ' for path/row ' + pathrow[0:3] + '/' + pathrow[3:6] + ' with ' + str(
......
......@@ -5,9 +5,11 @@ import platform
import subprocess
import sys
import ogr
import gdal
import warnings
from shutil import copyfile
from datetime import datetime,timedelta
from mtdUtils import checkSRS
from mtdUtils import checkSRS, setNoDataValue, getNoDataValue
def main(argv):
......@@ -19,7 +21,7 @@ def main(argv):
sys.exit("Platform not supported!")
try:
opts, args = getopt.getopt(argv, '', ['onstack','extent='])
opts, args = getopt.getopt(argv, '', ['onstack','extent=','force'])
except getopt.GetoptError as err:
print str(err)
......@@ -27,13 +29,16 @@ def main(argv):
dirs = [imroot + '/' + n for n in os.listdir(imroot) if os.path.isdir(imroot + '/' + n)]
ext_fn = None
force = False
onstack = False
for opt,val in opts:
if opt == '--onstack':
onstack = True
if opt == '--extent':
ext_fn = val
if opt == '--force':
force = True
if onstack and (ext_fn is None or not os.path.exists(ext_fn)):
sys.exit("Missing or invalid extent file.")
......@@ -47,10 +52,10 @@ def main(argv):
for d in dirs:
qq = os.path.basename(d).split('_')
if qq[0] in ["S2","L8"]:
if qq[0][0:2] in ["S2","L8"]:
if sen is None:
sen = qq[0]
sdt.append(tuple(qq[1:] + [d]))
sen = qq[0][0:2]
sdt.append(tuple(qq[-2:] + [d]))
sdt.sort(key=lambda l: l[0])
......@@ -83,17 +88,33 @@ def main(argv):
if not os.path.exists(mosd):
os.mkdir(mosd)
for mos in toMosaic:
od = mosd + os.path.basename(mos[0])[0:-5] + 'MOSAIC'
if not os.path.exists(od):
os.mkdir(od)
if sen == 'S2':
obd = os.path.basename(mos[0])[0:-5] + 'MOSAIC'
if 'THEIA' in obd:
obd = obd[0:13] + obd[15:]
od = mosd + obd
if not os.path.exists(od):
os.mkdir(od)
else:
if not force:
warnings.warn('S2: directory ' + od + 'exists, skipping.')
continue
if onstack:
fns = []
for d in mos:
if not os.path.exists(od + '/tileInfo.json'):
if os.path.exists(d + '/tileInfo.json') and not os.path.exists(od + '/tileInfo.json'):
copyfile(d + '/tileInfo.json', od + '/tileInfo.json')
msk_fn = glob.glob(d + '/' + 'MASK*.tif')[0]
stack_fn = [ff for ff in glob.glob(d + '/' + 'STACK*.tif') if '_MASKED' not in ff][0]
else:
mtd = glob.glob(d + '/*_MTD_ALL.xml')
if len(mtd) == 1 and not os.path.exists(od + '/' + 'MOSAIC_MTD_ALL.xml'):
copyfile(mtd[0], od + '/' + 'MOSAIC_MTD_ALL.xml')
#msk_fn = glob.glob(d + '/' + '*MASK*.tif')[0]
lst = glob.glob(d + '/' + '*MASK*.tif')
if 'GAPMASK' in lst[0]:
msk_fn = lst[1]
else:
msk_fn = lst[0]
stack_fn = [ff for ff in glob.glob(d + '/' + '*STACK*.tif') if '_MASKED' not in ff][0]
#stack_fn = glob.glob(d + '/' + 'STACK*.tif')[0]
mstack_fn = stack_fn.replace('STACK','STACK_MASKED')
subprocess.call(['otbcli_ManageNoData', '-in', stack_fn, '-mode', 'apply', '-mode.apply.mask', msk_fn, '-out', mstack_fn, 'uint16'],shell=sh)
......@@ -111,7 +132,7 @@ def main(argv):
gmfn = od + '/' + os.path.basename(fns[0]).replace('STACK_MASKED', 'GAPMASK')
cmd = ['otbcli_Mosaic', '-il'] + fns + ['-out', ofn, 'uint16']
subprocess.call(cmd, shell=sh)
cmd = ['otbcli_ManageNoData', '-in', ofn, '-out', mfn, '-mode', 'buildmask']
cmd = ['otbcli_ManageNoData', '-in', ofn, '-out', mfn, 'uint8', '-mode', 'buildmask']
subprocess.call(cmd, shell=sh)
if not os.path.exists(env_fn):
cmd = ['otbcli_Rasterization', '-in', ext_fn, '-out', env_fn, 'uint8',
......@@ -122,6 +143,7 @@ def main(argv):
for f in fns:
os.remove(f)
else:
#Following does not work for THEIA images!!
bnds = ['B01','B02','B03','B04','B05','B06','B07','B08','B8A','B09','B10','B11','B12']
for b in bnds:
fns = []
......@@ -133,12 +155,75 @@ def main(argv):
ofn = os.path.basename(fns[0])
cmd = ['otbcli_Mosaic', '-il'] + fns + ['-out', od + '/' + ofn, 'uint16']
subprocess.call(cmd, shell=sh)
elif sen == "L8":
od = mosd + os.path.basename(mos[0])[0:-6] + 'MOSAIC'
if not os.path.exists(od):
os.mkdir(od)
else:
if not force:
warnings.warn('L8: directory ' + od + 'exists, skipping.')
continue
if onstack:
fns = []
for d in mos:
if not os.path.exists(od + '/MOSAIC_MTL.txt'):
mtl = glob.glob(d + '/*_MTL*.txt')[0]
copyfile(mtl, od + '/MOSAIC_MTL.txt')
msk_fn = glob.glob(d + '/' + '*_MASK_HR_*.TIF')[0]
stack_fn = [ff for ff in glob.glob(d + '/' + '*_PS_TOA_*.TIF') if '_MASKED' not in ff][0]
mstack_fn = stack_fn.replace('PS_TOA', 'PS_TOA_MASKED')
ndv = getNoDataValue(stack_fn)
setNoDataValue(stack_fn,0)
cmd = ['otbcli_ManageNoData', '-in', stack_fn, '-mode', 'apply', '-mode.apply.mask', msk_fn,
'-out', mstack_fn, 'uint16']
subprocess.call(cmd, shell=sh)
setNoDataValue(stack_fn,ndv)
if ext_srs_epsg is None:
ext_srs_epsg = checkSRS(mstack_fn)
elif checkSRS(mstack_fn) != ext_srs_epsg:
mstack_fn_rep = mstack_fn.replace('.TIF', '_rep.TIF')
cmd = ['gdalwarp', '-t_srs', ext_srs_epsg, mstack_fn, mstack_fn_rep]
subprocess.call(cmd, shell=sh)
os.remove(mstack_fn)
os.rename(mstack_fn_rep, mstack_fn)
fns.append(mstack_fn)
ofn = od + '/' + os.path.basename(fns[0]).replace('PS_TOA_MASKED', 'PS_TOA')
mfn = od + '/' + os.path.basename(fns[0]).replace('PS_TOA_MASKED', 'MASK_HR')
gmfn = od + '/' + os.path.basename(fns[0]).replace('PS_TOA_MASKED', 'GAPMASK')
cmd = ['otbcli_Mosaic', '-il'] + fns + ['-out', ofn, 'uint16']
subprocess.call(cmd, shell=sh)
#cmd = ['otbcli_ManageNoData', '-in', ofn, '-out', mfn, '-mode', 'buildmask']
#No-data checked in B4 only (avoid problems with B9 in TOA)
cmd = ['otbcli_BandMathX', '-il', ofn, '-out', mfn, 'uint8', '-exp', '{im1b4 != ' + str(getNoDataValue(ofn)) + '}']
subprocess.call(cmd, shell=sh)
if not os.path.exists(env_fn):
cmd = ['otbcli_Rasterization', '-in', ext_fn, '-out', env_fn, 'uint8',
'-im', ofn, '-mode.binary.foreground', '1']
subprocess.call(cmd, shell=sh)
cmd = ['otbcli_BandMathX', '-il', env_fn, mfn, '-exp', '{im1b1 != 0 && im2b1 == 0}', '-out', gmfn,
'uint8']
subprocess.call(cmd, shell=sh)
for f in fns:
os.remove(f)
else:
bnds = ['B1', 'B2', 'B3', 'B4', 'B5', 'B6', 'B7', 'B8', 'B9', 'B10', 'B11']
for b in bnds:
fns = []
for d in mos:
# TO DO: manage different projections for band-by-band mode
if not os.path.exists(od + '/MOSAIC_MTL.txt'):
mtl = glob.glob(d + '/*_MTL*.txt')[0]
copyfile(mtl, od + '/MOSAIC_MTL.txt')
fns.append(glob.glob(d + '/' + b + '*.TIF')[0])
ofn = os.path.basename(fns[0])
cmd = ['otbcli_Mosaic', '-il'] + fns + ['-out', od + '/' + ofn, 'uint16']
subprocess.call(cmd, shell=sh)
<