Commit c88d4fde authored by Lozac'h Loic's avatar Lozac'h Loic
Browse files

sync repo

parent 8651b22e
#!/usr/bin/python
# To change this license header, choose License Headers in Project Properties.
# To change this template file, choose Tools | Templates
# and open the template in the editor.
import os
import argparse
import otbApplication
def search_files(directory='.', resolution='NDVI', extension='tif', fictype='f'):
images=[]
extension = extension.lower()
resolution = resolution.lower()
for dirpath, dirnames, files in os.walk(directory):
if fictype == 'f':
for name in files:
# print(os.path.join(dirpath, name) + " test")
if extension and name.lower().endswith(extension) and name.lower().find(resolution) >= 0 :
# print(os.path.join(dirpath, name) + " OK")
abspath = os.path.abspath(os.path.join(dirpath, name))
images.append(abspath)
elif fictype == 'd':
for dirname in dirnames:
# print(os.path.join(dirpath, name) + " test")
if extension and dirname.lower().endswith(extension) and dirname.lower().find(resolution) >= 0 :
# print(os.path.join(dirpath, name) + " OK")
abspath = os.path.abspath(os.path.join(dirpath, dirname))
images.append(abspath)
else:
print("search_files type error")
exit
return images
def ex_pipeline(argsformat, extr, ref, out):
print("processing.......................................................")
print("image = "+extr)
print("output = "+out)
app0 = otbApplication.Registry.CreateApplication("ExtractROI")
app0.SetParameterString("in",extr)
app0.SetParameterString("out", "temp0.tif")
app0.SetParameterString("mode","fit")
app0.SetParameterString("mode.fit.im",ref)
app0.Execute()
app3 = otbApplication.Registry.CreateApplication("ManageNoData")
app3.SetParameterInputImage("in", app0.GetParameterOutputImage("out"))
app3.SetParameterString("out", out+"?gdal:co:COMPRESS=DEFLATE&gdal:co:BIGTIFF=YES")
if argsformat == 'L8-NDVI':
app3.SetParameterOutputImagePixelType("out", otbApplication.ImagePixelType_uint8)
elif argsformat == 'S1-SNAP':
app3.SetParameterOutputImagePixelType("out", otbApplication.ImagePixelType_float)
else:
print("error!")
exit
app3.SetParameterString("mode", "changevalue")
app3.ExecuteAndWriteOutput()
print("done.............................................................")
def exmos_pipeline(argsformat, mos1, mos2, ref, out):
print("processing.......................................................")
print("image1 = "+mos1)
print("image2 = "+mos2)
print("output = "+out)
app0 = otbApplication.Registry.CreateApplication("Superimpose")
app0.SetParameterString("inm",mos1)
app0.SetParameterString("out", "temp0.tif")
app0.SetParameterString("interpolator","nn")
app0.SetParameterString("inr",ref)
app0.Execute()
app1 = otbApplication.Registry.CreateApplication("Superimpose")
app1.SetParameterString("inm",mos2)
app1.SetParameterString("out", "temp1.tif")
app1.SetParameterString("interpolator","nn")
app1.SetParameterString("inr",ref)
app1.Execute()
app2 = otbApplication.Registry.CreateApplication("Mosaic")
app2.AddImageToParameterInputImageList("il",app0.GetParameterOutputImage("out"))
app2.AddImageToParameterInputImageList("il", app1.GetParameterOutputImage("out"))
app2.SetParameterString("out", "temp2.tif")
app2.Execute()
app3 = otbApplication.Registry.CreateApplication("ManageNoData")
app3.SetParameterInputImage("in", app2.GetParameterOutputImage("out"))
app3.SetParameterString("out", out+"?gdal:co:COMPRESS=DEFLATE&gdal:co:BIGTIFF=YES")
if argsformat == 'L8-NDVI':
app3.SetParameterOutputImagePixelType("out", otbApplication.ImagePixelType_uint8)
elif argsformat == 'S1-SNAP':
app3.SetParameterOutputImagePixelType("out", otbApplication.ImagePixelType_float)
else:
print("error!")
exit
app3.SetParameterString("mode", "changevalue")
app3.ExecuteAndWriteOutput()
print("done.............................................................")
if __name__ == "__main__":
# Make parser object
parser = argparse.ArgumentParser(description=
"""
Extract and mosaic Sentinel-1 or Landsat8 to fit Sentinel-2 Tile
""")
parser.add_argument('-format', choices=['S1-SNAP','L8-NDVI'], required=True)
parser.add_argument('-indir', action='store', required=True, help='Directory containing Sentinel1 SNAP Calibrated or Landsat8 NDVI')
parser.add_argument('-inref', action='store', required=True, help='Image references Sentinel-2 Tile, example any S2 NDVI')
parser.add_argument('-outdir', action='store', required=True, help='Output directory')
# parser.add_argument('-user', action='store', required=True, help='SciHub Copernicus user name')
# parser.add_argument('-psswd', action='store', required=True, help='SciHub Copernicus user password')
# parser
args=parser.parse_args()
if not os.path.isdir(args.outdir):
os.mkdir(args.outdir)
abspathout = os.path.abspath(args.outdir)
splitinref = os.path.basename(args.inref).split("_")
tileid = splitinref[1]
if args.format == "S1-SNAP" :
indir=[]
indir=search_files(args.indir, 'S1', 'data', 'd')
dejafait=[]
for file in indir:
if file in dejafait:
continue
bname = os.path.basename(file).split(".")[0]
splitbname = bname.split("_")
acqdate = splitbname[4].split("T")[0]
acqtime = splitbname[4].split("T")[1]
outfilebase = os.path.join(abspathout,"_".join(splitbname[:3])+"_"+tileid+"_"+acqdate+"T"+acqtime)
l = [indir.index(i) for i in indir if acqdate in i]
if len(l) == 2 :
in1 = os.path.join(indir[l[0]], "Sigma0_VV.img")
in2 = os.path.join(indir[l[1]], "Sigma0_VV.img")
outfile = outfilebase +"_VV.TIF"
if os.path.exists(outfile) :
continue
exmos_pipeline(args.format, in1, in2, args.inref, outfile)
in1 = os.path.join(indir[l[0]], "incidenceAngleFromEllipsoid.img")
in2 = os.path.join(indir[l[1]], "incidenceAngleFromEllipsoid.img")
outfile = outfilebase +"_THETA.TIF"
exmos_pipeline(args.format, in1, in2, args.inref, outfile)
dejafait.append(indir[l[0]])
dejafait.append(indir[l[1]])
elif len(l) == 1 :
in1 = os.path.join(file, "Sigma0_VV.img")
outfile = outfilebase +"_VV.TIF"
if os.path.exists(outfile) :
continue
ex_pipeline(args.format, in1, args.inref, outfile)
in1 = os.path.join(file, "incidenceAngleFromEllipsoid.img")
outfile = outfilebase +"_THETA.TIF"
ex_pipeline(args.format, in1, args.inref, outfile)
else :
print("Error on index :"+l)
elif args.format == "L8-NDVI" :
indir=[]
indir=search_files(args.indir, 'NDVI_LC08', 'tif', 'f')
for file in indir:
bname = os.path.basename(file).split(".")[0]
splitbname = bname.split("_")
acqdate = splitbname[-1]
outfile = os.path.join(abspathout,"NDVI_"+tileid+"_"+acqdate+"T000000_"+splitbname[1]+".TIF")
if os.path.exists(outfile) :
continue
l = [indir.index(i) for i in indir if acqdate in i]
if len(l) == 2 :
exmos_pipeline(args.format, indir[l[0]], indir[l[1]], args.inref, outfile)
elif len(l) == 1 :
ex_pipeline(args.format, file, args.inref, outfile)
else :
print("Error on index :"+l)
else:
print("S2 format not recognized!")
exit
\ No newline at end of file
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Wed Aug 28 13:55:34 2019
@author: llozach
"""
import os,datetime, math
import argparse
import otbApplication
def search_files(directory='.', extension='tif', resolution='NDVI'):
images=[]
extension = extension.lower()
resolution = resolution.lower()
for dirpath, dirnames, files in os.walk(directory):
for name in files:
if extension and name.lower().endswith(extension) and name.lower().find(resolution) >= 0 :
#print(os.path.join(dirpath, name))
abspath = os.path.abspath(os.path.join(dirpath, name))
images.append(abspath)
print(str(len(images)) + " image(s) found")
return images
def create_mask(ndvi,mask):
app1 = otbApplication.Registry.CreateApplication("BandMath")
app1.AddParameterStringList("il",ndvi)
app1.SetParameterString("out", mask)
app1.SetParameterString("exp", "im1b1==0?1:0")
app1.SetParameterOutputImagePixelType("out", otbApplication.ImagePixelType_uint8)
app1.ExecuteAndWriteOutput()
def concate_images(listinputimg, pixeltype, output):
apps=[]
for img in listinputimg:
app0 = otbApplication.Registry.CreateApplication("Superimpose")
app0.SetParameterString("inm", img)
app0.SetParameterString("out", "supimp.tif")
app0.SetParameterString("inr", listinputimg[0])
app0.SetParameterString("interpolator", "nn")
app0.Execute()
apps.append(app0)
print('---\nConcatenation to {}...'.format(output))
appConcat = otbApplication.Registry.CreateApplication("ConcatenateImages")
for a in apps:
appConcat.AddImageToParameterInputImageList("il",a.GetParameterOutputImage("out"))
appConcat.SetParameterString("out", output)
#appConcat.SetParameterOutputImagePixelType("out", pixeltype)
appConcat.Execute()
return appConcat
def monthrange(start_date, end_date):
# print("MonthRange: "+str (math.floor((end_date - start_date).days/31)))
# for n in range(int (math.floor((end_date - start_date).days/31))+1):
# yield start_date + datetime.timedelta(days=n*31)
qannee = end_date.year - start_date.year
dmois = end_date.month +12*(qannee) - start_date.month
yield start_date
anneesuiv = start_date.year
for m in range(dmois-1):
m+=1
moissuiv = (start_date.month + m) % 12
if moissuiv == 0 :
moissuiv = 12
if moissuiv == 1 :
anneesuiv = anneesuiv+1
yield datetime.date(anneesuiv,moissuiv,1)
yield end_date
def create_outputdates(listinput_date):
# stdate = datetime.date(int(startdate[:4]),int(startdate[4:6]),1)
# eddate = datetime.date(int(enddate[:4]),int(enddate[4:6]),15)
inst = datetime.date(int(listinput_date[0][:4]),int(listinput_date[0][4:6]),int(listinput_date[0][6:8]))
ined = datetime.date(int(listinput_date[-1][:4]),int(listinput_date[-1][4:6]),int(listinput_date[-1][6:8]))
if inst.day > 15:
stdate = datetime.date(inst.year,inst.month +1,1)
else:
stdate = datetime.date(inst.year,inst.month,15)
if ined.day >= 15:
eddate = datetime.date(ined.year,ined.month,15)
else:
eddate = datetime.date(ined.year,ined.month,1)
exformat = "%Y%m%d"
print("Start date: "+stdate.strftime(exformat))
print("End date: "+eddate.strftime(exformat))
outputdates=[]
#outputdates.append(stdate.strftime(exformat))
procd = stdate
for dstep in monthrange(stdate, eddate):
procd = dstep
if dstep == stdate:
if stdate.day == 1:
outputdates.append(procd.strftime(exformat))
procd = datetime.date(dstep.year,dstep.month,15)
outputdates.append(procd.strftime(exformat))
else:
outputdates.append(procd.strftime(exformat))
elif dstep == eddate:
if eddate.day == 1:
outputdates.append(procd.strftime(exformat))
else:
procd = datetime.date(dstep.year,dstep.month,1)
outputdates.append(procd.strftime(exformat))
procd = datetime.date(dstep.year,dstep.month,15)
outputdates.append(procd.strftime(exformat))
else:
procd = datetime.date(dstep.year,dstep.month,1)
outputdates.append(procd.strftime(exformat))
procd = datetime.date(dstep.year,dstep.month,15)
outputdates.append(procd.strftime(exformat))
return outputdates
if __name__ == "__main__":
# Make parser object
parser = argparse.ArgumentParser(description=
"""
Perform Gapfilling over NDVI time serie.
""")
parser.add_argument('-ndvidir', action='store', required=True, help='Directory containing time serie NDVI to be gapfilled')
parser.add_argument('-outmaskdir', action='store', required=True, help='Output directory for NDVIs mask used for gapfilling')
parser.add_argument('-interpolation', choices=['linear', 'spline'], default='linear', required=False, help='Interpolation mode for OTB ImageTimeSeriesGapFilling')
parser.add_argument('-outdir', action='store', required=True, help='Output directory for time serie gapfilled NDVIs image')
# parser.add_argument('-user', action='store', required=True, help='SciHub Copernicus user name')
# parser.add_argument('-psswd', action='store', required=True, help='SciHub Copernicus user password')
args=parser.parse_args()
outdir = args.outdir
ndvifiles=search_files(args.ndvidir)
parameters=[]
for ndvi in ndvifiles:
maskname = os.path.join(args.outmaskdir,os.path.basename(ndvi)[:-4]+"_MASK.TIF")
parameters.append([os.path.basename(ndvi).split("_")[2].split("T")[0],ndvi,maskname])
parameters.sort(key=lambda tup: tup[0])
print("Creating MASK files from input NDVIs")
for param in parameters:
if os.path.exists(param[2]):
continue
create_mask(param[1],param[2])
listinput_date = [row[0] for row in parameters]
listinput_img = [row[1] for row in parameters]
listinput_mask = [row[2] for row in parameters]
listoutput_date = create_outputdates(listinput_date)
inputdatesfile = os.path.join(args.ndvidir,"InputDates.txt")
outputdatesfile = os.path.join(args.ndvidir,"OutputDates.txt")
with open(inputdatesfile,'w') as indatefile:
for text in listinput_date:
indatefile.write(text+"\n")
with open(outputdatesfile,'w') as outdatefile:
for text in listoutput_date:
outdatefile.write(text+"\n")
# Concatenate
print("Concatenate NDVI files : ")
[print(" "+nm) for nm in listinput_img]
#appConcatImg = concate_images(listinput_img, otbApplication.ImagePixelType_uint8, "tempimg.tif")
apps=[]
for img in listinput_img:
app0 = otbApplication.Registry.CreateApplication("Superimpose")
app0.SetParameterString("inm", img)
app0.SetParameterString("out", "supimp.tif")
app0.SetParameterString("inr", listinput_img[0])
app0.SetParameterString("interpolator", "nn")
app0.Execute()
apps.append(app0)
appConcatImg = otbApplication.Registry.CreateApplication("ConcatenateImages")
for a in apps:
appConcatImg.AddImageToParameterInputImageList("il",a.GetParameterOutputImage("out"))
# appConcatImg.SetParameterStringList("il", listinput_img)
appConcatImg.SetParameterString("out", "tempimg.tif")
appConcatImg.Execute()
print("Concatenate MASK files : ")
[print(" "+nm) for nm in listinput_mask]
#appConcatMask = concate_images(listinput_mask, otbApplication.ImagePixelType_uint8, "tempmask.tif")
appsm=[]
for img in listinput_mask:
app0 = otbApplication.Registry.CreateApplication("Superimpose")
app0.SetParameterString("inm", img)
app0.SetParameterString("out", "supimp.tif")
app0.SetParameterString("inr", listinput_img[0])
app0.SetParameterString("interpolator", "nn")
app0.Execute()
appsm.append(app0)
appConcatMask = otbApplication.Registry.CreateApplication("ConcatenateImages")
for a in appsm:
appConcatMask.AddImageToParameterInputImageList("il",a.GetParameterOutputImage("out"))
# appConcatMask.SetParameterStringList("il", listinput_mask)
appConcatMask.SetParameterString("out", "tempmask.tif")
appConcatMask.Execute()
# Gapfill
print("Gapfilled NDVIs from input dates in "+inputdatesfile+" : ")
[print(" "+nm) for nm in listinput_date]
print("Gapfilled NDVIs to output dates in "+outputdatesfile+" : ")
[print(" "+nm) for nm in listoutput_date]
appGapFill = otbApplication.Registry.CreateApplication("ImageTimeSeriesGapFilling")
appGapFill.SetParameterInputImage("in", appConcatImg.GetParameterOutputImage("out"))
appGapFill.SetParameterInputImage("mask", appConcatMask.GetParameterOutputImage("out"))
appGapFill.SetParameterInt("comp", 1)
appGapFill.SetParameterString("it", args.interpolation)
appGapFill.SetParameterString("id", inputdatesfile)
appGapFill.SetParameterString("od", outputdatesfile)
appGapFill.SetParameterString("out", "tempgapfill.tif")
appGapFill.Execute()
prefsp = os.path.basename(ndvifiles[0]).split('_')
outprefix = os.path.join(args.outdir,prefsp[0]+"_"+prefsp[1])
# Splitimage
print("Spliting gapfilled image to monoband images")
SplitImage = otbApplication.Registry.CreateApplication("SplitImage")
SplitImage.SetParameterInputImage("in", appGapFill.GetParameterOutputImage("out"))
SplitImage.SetParameterString("out", outprefix+".TIF")
SplitImage.SetParameterOutputImagePixelType("out", otbApplication.ImagePixelType_uint8)
SplitImage.ExecuteAndWriteOutput()
# Rename image
print("Renaming files with dates")
i=0
for gapdate in listoutput_date:
os.rename(outprefix+"_"+str(i)+".TIF",outprefix+"_"+gapdate+".TIF")
i+=1
print("Done")
\ No newline at end of file
#!/usr/bin/python
# To change this license header, choose License Headers in Project Properties.
# To change this template file, choose Tools | Templates
# and open the template in the editor.
__author__ = "Loic Lozach"
__date__ = "$Dec 14, 2018 12:40:21 PM$"
import os, argparse, subprocess, math, datetime,glob
import otbApplication
from osgeo import gdal, osr, ogr
from subprocess import Popen, PIPE
maskedlabels, labels, vectorlabels, labelsfield, lulc, agrivalues, ndvi, sarvv, \
slope, slopevalue, sarth, modeldir, output, outformat, sarmode = "","","","","","","","","","","","","","",""
def search_files(directory='.', resolution='NDVI', extension='tif', fictype='f'):
images=[]
extension = extension.lower()
resolution = resolution.lower()
for dirpath, dirnames, files in os.walk(directory):
if fictype == 'f':
for name in files:
# print(os.path.join(dirpath, name) + " test")
if extension and name.lower().endswith(extension) and name.lower().find(resolution) >= 0 :
# print(os.path.join(dirpath, name) + " OK")
abspath = os.path.abspath(os.path.join(dirpath, name))
images.append(abspath)
elif fictype == 'd':
for dirname in dirnames:
# print(os.path.join(dirpath, name) + " test")
if extension and dirname.lower().endswith(extension) and dirname.lower().find(resolution) >= 0 :
# print(os.path.join(dirpath, name) + " OK")
abspath = os.path.abspath(os.path.join(dirpath, dirname))
images.append(abspath)
else:
print("search_files type error")
exit()
return images
def reprojectRaster(absfile, indeximg):
global sarvv, ndvi, maskedlabels, slope
ext=absfile.split(".")[1]
longname=absfile.split(".")[0]
if ext.lower() == "img" :
temp_file = absfile
absfile = longname+".tif"
else :
temp_file = longname + ".old." + ext
os.rename(absfile, temp_file)
p = Popen(['gdalwarp', '-t_srs', 'EPSG:4326', temp_file, absfile], stdout=PIPE)
# p.wait()
output = p.communicate()[0]
if p.returncode != 0:
print("gdalwarp failed %d : %s" % (p.returncode, output))
# with open("gdalwrap_log.err",'a') as err:
# err.write("######################################################################################################\n")
# err.write(absfile)
# err.write(output)
# err.write("\n######################################################################################################")
if ext.lower() == "img" :
os.remove(absfile)
else :
os.remove(absfile)
os.rename(temp_file, absfile)
return None, 1
if ext.lower() == "img" :
if indeximg == 0 :
sarvv = absfile
elif indeximg == 1 :
ndvi = absfile
elif indeximg == 2 :
maskedlabels = absfile
elif indeximg == 3 :
slope = absfile
print("gdalwarp succeeded on : "+absfile)
os.remove(temp_file)
return absfile, 0