#Sentinel 3 - SLSTR (Sea and Land Surface Temperature Radiometer) import
#Sentinel 3 data from Seninel 3 Pre-Operations Data Hub, available at https://scihub.copernicus.eu/s3/#/home
#Works with 1 km resolution nadir data for both Level 1 RBT and Level 2 LST and SST

from tkinter import *
from tkinter.filedialog import askdirectory
from tkinter.filedialog import askopenfilename
from tkinter.ttk import Progressbar, Style
import os
import sys
import numpy as np
from netCDF4 import Dataset
from osgeo import gdal
import ilwis
import shutil
from threading import Thread

#passing parameters from sub menu

def execute(SentDir, InFileName, WorkDir, OutFileName, exportFormat):
    print(exportFormat)
    if len(SentDir) > 0 and len(InFileName) > 0 and len(WorkDir) > 0 and len(OutFileName) > 0:

        progress['value'] = 0

        os.chdir(WorkDir)

        #process geodetic information
        #using geodetic_in.nc

        nc_f = os.path.join(SentDir,'geodetic_in.nc')
        shutil.copy(SentDir+'\\'+InFileName, WorkDir)

        nc_fid = Dataset(nc_f, 'r')
        lats = nc_fid.variables['latitude_in'][:]
        lons = nc_fid.variables['longitude_in'][:]

        height, width = lats.shape
        offset = 283.73
        scale = 0.01
        step = 20

        #create gdal vrt file
        vrt_file = open(OutFileName + '.vrt', "w")
        vrt_file.write('<VRTDataset rasterXSize="' + str(width) + '" rasterYSize="' + str(height) + '">\n')
        vrt_file.write('  <Metadata>\n')
        vrt_file.write('    <MDI key="resolution">[ 1000 1000 ]</MDI>\n')
        vrt_file.write('    <MDI key="title">S3 SLSTR L1 Radiance and Brightness Temperatures Product. (Measurements + Annotations)</MDI>\n')
        vrt_file.write('  </Metadata>\n')
        vrt_file.write('  <GCPList Projection="GEOGCS[&quot;WGS 84&quot;,DATUM[&quot;WGS_1984&quot;,SPHEROID[&quot;WGS 84&quot;,6378137,298.257223563,AUTHORITY[&quot;EPSG&quot;,&quot;7030&quot;]],AUTHORITY[&quot;EPSG&quot;,&quot;6326&quot;]],PRIMEM[&quot;Greenwich&quot;,0,AUTHORITY[&quot;EPSG&quot;,&quot;8901&quot;]],UNIT[&quot;degree&quot;,0.0174532925199433,AUTHORITY[&quot;EPSG&quot;,&quot;9122&quot;]],AUTHORITY[&quot;EPSG&quot;,&quot;4326&quot;]]">\n')
        for i in range(0,width,step):
            for j in range(0,height,step):
                pixel = i + 0.5
                line = j + 0.5
                X = lons[j,i]
                Y = lats[j,i]
                vrt_file.write('    <GCP Id="" Pixel="' + str(pixel) + '" Line="' + str(line) + '" X="{:.16f}" Y="{:.16f}" />\n'.format(X,Y))
        vrt_file.write('  </GCPList>\n')
        vrt_file.write('  <VRTRasterBand dataType="Int16" band="1">\n')
        vrt_file.write('    <Metadata>\n')
        vrt_file.write('      <MDI key="'+(InFileName[:-3])+'_add_offset">' + str(offset) + ' </MDI>\n')
        vrt_file.write('      <MDI key="'+(InFileName[:-3])+'_long_name">Gridded pixel brightness temperature for channel '+(InFileName)+' (1km TIR grid, nadir view)</MDI>\n')
        vrt_file.write('      <MDI key="'+(InFileName[:-3])+'_scale_factor">' + str(scale) + ' </MDI>\n')
        vrt_file.write('      <MDI key="'+(InFileName[:-3])+'_standard_name">toa_brightness_temperature</MDI>\n')
        vrt_file.write('      <MDI key="'+(InFileName[:-3])+'_units">K</MDI>\n')
        vrt_file.write('      <MDI key="'+(InFileName[:-3])+'n__FillValue">-32768 </MDI>\n')
        vrt_file.write('    </Metadata>\n')
        vrt_file.write('    <SimpleSource>\n')
        vrt_file.write('      <SourceFilename relativeToVRT="1">HDF5:'+(InFileName)+'://'+(InFileName[:-3])+'</SourceFilename>\n')
        vrt_file.write('      <SourceBand>1</SourceBand>\n')
        vrt_file.write('      <SourceProperties RasterXSize="' + str(width) + '" RasterYSize="' + str(height) + '" DataType="Int16" BlockXSize="' + str(width) + '" BlockYSize="' + str(height) + '" />\n')
        vrt_file.write('      <SrcRect xOff="0" yOff="0" xSize="' + str(width) + '" ySize="' + str(height) + '" />\n')
        vrt_file.write('      <DstRect xOff="0" yOff="0" xSize="' + str(width) + '" ySize="' + str(height) + '" />\n')
        vrt_file.write('    </SimpleSource>\n')
        vrt_file.write('  </VRTRasterBand>\n')
        vrt_file.write('</VRTDataset>\n')
        vrt_file.close()


        progress['value'] = 40

        #conduct gdal warp operation using vrt file created and save file as ilwis mpr
        input_raster = gdal.Open(OutFileName + '.vrt')
        output_raster = (OutFileName+'.mpr')
        gdal.Warp(output_raster,input_raster,format='ILWIS')
        os.remove ((WorkDir.replace('\\','/')) + '/' + OutFileName +'.vrt')

        progress['value'] = 65 

        #final radiometric correction applying ilwis-objects - loaded as plugin under pyhton
        ilwis.Engine.setWorkingCatalog('file:///'+ WorkDir)
        rc = ilwis.RasterCoverage(output_raster)
        print (rc.name())
        if exportFormat == 'ilwis':
            rc2=ilwis.Engine.do('mapcalc','iff((@1 != -32768) and (@1 != 0),@1*'+str(scale)+'+'+str(offset)+',?)',rc)
            rc2.store((OutFileName) +'.mpr', 'map', 'ilwis3')
            os.remove ((WorkDir.replace('\\','/')) + '/' + OutFileName +'.mpr.aux.xml')

        elif exportFormat == 'qgis':
        #print ('Output file in QGIS-Geotiff format created = '+ OutFileName +'.tif')
        #set nodata (NaN) (of surrounding area) to -999 for display in QGIS
            rc2=ilwis.Engine.do('mapcalc','iff((@1 != -32768) and (@1 != 0),@1*'+str(scale)+'+'+str(offset)+',-999)',rc)
            rc2.store((OutFileName) +'.tif', "GTiff", "gdal")
            os.remove ((WorkDir.replace('\\','/')) + '/' + OutFileName +'.mpr.aux.xml')
            os.remove((WorkDir.replace('\\', '/'))+'/'+ OutFileName +'.csy')
            os.remove((WorkDir.replace('\\', '/'))+'/'+ OutFileName +'.grf')
            os.remove((WorkDir.replace('\\', '/'))+'/'+ OutFileName +'.mpr')
            os.remove((WorkDir.replace('\\', '/'))+'/'+ OutFileName +'.mp#')
            


        progress['value'] = 100
    button0.config(state='normal')
    button1.config(state='normal')

#redundant files deleted in main menu script

def executeInThread(SentDir, InFileName, WorkDir, OutFileName, exportFormat):
    button0.config(state='disabled')
    button1.config(state='disabled')
    thread = Thread(target = execute, args = (SentDir, InFileName, WorkDir, OutFileName, exportFormat))
    thread.start()

#main()

sent_path = os.path.realpath('.')
print (sent_path)
os.chdir(sent_path)
os.environ["sent_lib"] = sent_path

app = Tk()
app.title("Import S3 SLSTR-Land-BT channel F")
app.iconbitmap(sent_path +'\pics\esa.ico')

w = 650 # width for the Tk app
h = 210 # height for the Tk app

app.geometry('%dx%d' % (w, h))
app.resizable(0,0)

labelText1 = StringVar()
labelText1.set ("""S3 SLSTR L1 Radiance and Brightness Temperatures Product, brightness temperature for channel F
Typical filename: S3A_SL_1_RBT____*_SVL_O_NR_002.SEN3\F1_BT_in.nc and also geodetic_in.nc required""")

label1 = Label(app, textvariable=labelText1, anchor='w', justify = LEFT)
label1.pack(side = LEFT)
label1.place(x = 10, y = 5, width=650, height=35)

L1 = Label(app, text="Input Directory:", anchor='w')
L1.pack(side = LEFT)
L1.place(x = 10, y = 50, width=150, height=25)

inputDir = StringVar()
E1 = Entry(app, textvariable=inputDir, bd =2)
E1.pack()
E1.place(x = 165, y = 50, width=445, height=25)

button_dirin = Button(app, text="....", width=2, command = lambda: inputDir.set(askdirectory(initialdir = inputDir.get())))
button_dirin.pack()
button_dirin.place(x = 615, y = 49)

L2 = Label(app, text="Input Filename:", anchor='w')
L2.pack(side = LEFT)
L2.place(x = 10, y = 80, width=150, height=25)

inputFile = StringVar()
E2 = Entry(app, textvariable=inputFile, bd =2)
E2.pack()
E2.place(x = 165, y = 80, width=445, height=25)

button_filein = Button(app, text="....", width=2, command = lambda: inputFile.set(askopenfilename(initialdir = inputDir.get()).split('/')[-1]))
button_filein.pack()
button_filein.place(x = 615, y = 79)

L3 = Label(app, text="Output Directory:", anchor='w')
L3.pack(side = LEFT)
L3.place(x = 10, y = 110, width=150, height=25)

outputDir = StringVar()
E3 = Entry(app, textvariable=outputDir, bd =2)
E3.pack()
E3.place(x = 165, y = 110, width=445, height=25)

button_dirout = Button(app, text="....", width=2, command = lambda: outputDir.set(askdirectory(initialdir = outputDir.get())))
button_dirout.pack()
button_dirout.place(x = 615, y = 109) 

L4 = Label(app, text="Output Filename:", anchor='w')
L4.pack(side = LEFT)
L4.place(x = 10, y = 140, width=150, height=25)

outputFile = StringVar()
E4 = Entry(app, textvariable=outputFile, bd =2)
E4.pack()
E4.place(x = 165, y = 140, width=445, height=25)

button0 = Button(app, text="Quit", width=15, command=app.destroy)
button0.pack()
button0.place(x = 525, y = 180)

exportFormat = StringVar()
   
R1 = Radiobutton(app, text="ILWIS", variable=exportFormat, value='ilwis')
R1.pack()
R1.place(x = 165, y = 180)

R2 = Radiobutton(app, text="QGIS", variable=exportFormat, value='qgis')
R2.pack()
R2.place(x = 285, y = 180)

R1.select()
R2.deselect()

button1 = Button(app, text="Execute", width=15, command=lambda: executeInThread(inputDir.get(), inputFile.get(), outputDir.get(), outputFile.get(), exportFormat.get()))
button1.pack()
button1.place(x = 405, y = 180)

style = Style()
style.theme_use('clam')
style.configure('blue.Horizontal.TProgressbar', foreground='blue', background='blue')

progress = Progressbar(app, style='blue.Horizontal.TProgressbar', orient=HORIZONTAL, length=100, mode='determinate')
progress.pack()
progress.place(x = 10, y = 180)
app.mainloop()

print("inputFile=" + inputFile.get())
print("outputDir=" + outputDir.get())
print("outputFile=" + outputFile.get())


