#Sentinel 5 - Level 1B import 
#Sentinel 5 data from Sentinel 5 Pre-Operations Data Hub, available at https://s5phub.copernicus.eu/dhus/#/home


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
import struct
from time import time

def writeIlwis(path, name, data, xsize, ysize):
    f=open(path + '\\' + name + '.mp#','wb')
    f.write(data)
    f.flush()
    f.close()
    f=open(path + '\\' + name + '.mpr','w')
    f.write('[Ilwis]\n')
    f.write('Description=Map ' + name + '\n')
    f.write('Time=' + str(int(time())) + '\n')
    f.write('Version=3.1\n')
    f.write('Class=Raster Map\n')
    f.write('Type=BaseMap\n')
    f.write('[BaseMap]\n')
    f.write('CoordSystem=unknown.csy\n')
    f.write('CoordBounds=0 -' + str(ysize) + ' ' + str(xsize) + ' -0\n')
    f.write('Domain=value.dom\n')
    f.write('DomainInfo=value.dom;Long;value;0;-9999999.9:9999999.9:0.1:offset=0;\n')
    f.write('Range=-1e+038:1e+038:0.000000:offset=0\n')
    f.write('ThreeDCoordinates=No\n')
    f.write('HistogramSize=?\n')
    f.write('Type=Map\n')
    f.write('Proximity=2.250000000000\n')
    f.write('[Map]\n')
    f.write('GeoRef=none.grf\n')
    f.write('Size=' + str(ysize) + ' ' + str(xsize) + '\n')
    f.write('Type=MapStore\n')
    f.write('[MapStore]\n')
    f.write('StoreTime=' + str(int(time())) + '\n')
    f.write('Data=' + name + '.mp#\n')
    f.write('Structure=Line\n')
    f.write('StartOffset=0\n')
    f.write('RowLength=' + str(xsize) + '\n')
    f.write('PixelInterLeaved=No\n')
    f.write('SwapBytes=No\n')
    f.write('Type=Float\n')
    f.write('UseAs=No\n')
    f.write('Format=2\n')
    f.flush()
    f.close()

def ephoton(wavelength):
    A         = 6.02214E23 # [mol-1]       Constant of Avogadro
    h         = 6.6262E-34 # [J s]         Planck's constant
    c         = 299792458 # [m s-1]        Speed of light
    return A*h*c/wavelength # [J]          Energy of 1 photon

#passing parameters from sub menu
def execute(SentDir, InFileName, inputInterval, 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,InFileName)
        shutil.copy(SentDir+'\\'+InFileName, WorkDir)
        
        wv=[(266.8,299.1),(300,320),(320,405),(397.8,498.6),(657,725),(725,775),(2305,2345),(2345,2385)]
        bandId=InFileName[18:19]
        #print('+bandId'+bandId)
                       

        nc_fid = Dataset(nc_f, 'r')
        lats = nc_fid.groups['BAND'+bandId+'_RADIANCE']['STANDARD_MODE']['GEODATA']['latitude'][:][0]
        rads = nc_fid.groups['BAND'+bandId+'_RADIANCE']['STANDARD_MODE']['OBSERVATIONS']['radiance'][:][0]
        height, width = lats.shape
        zsize=rads.shape[2]
        offset = 0
        scale = 1
        step = int(inputInterval)
        ilwis.Engine.setWorkingCatalog('file:///'+ WorkDir)
        progress['value'] = 10
		
        for z in range(0,zsize,step):
            print('Extracting band', int(z/step)+1, 'of', int(zsize/step))
            rads_band=rads[:,:,z]
            data=struct.pack('<%sf' % rads_band.size, *rads_band.flatten())
            bandFilename = OutFileName + '_Radiance_' + str(z)
            writeIlwis(WorkDir, bandFilename, data, width, height)

            #create gdal vrt file
            #check assignment of coordinate origin - here centre of pixel is assumed and NOT the corner
            vrt_file = open(bandFilename + '.vrt', "w")
            vrt_file.write('<VRTDataset rasterXSize="' + str(width) + '" rasterYSize="' + str(height) + '">\n')
            vrt_file.write('  <Metadata domain="GEOLOCATION">\n')
            vrt_file.write('    <MDI key="spatial_resolution">7x7km2</MDI>\n')
            vrt_file.write('    <MDI key="title">TROPOMI/S5P L1B Swath 7x7km</MDI>\n')
            vrt_file.write('    <MDI key="SRS">WGS84</MDI>\n')
            vrt_file.write('    <MDI key="X_DATASET">HDF5:"'+(InFileName)+'"://BAND'+bandId+'_RADIANCE/STANDARD_MODE/GEODATA/longitude</MDI>\n')
            vrt_file.write('    <MDI key="X_BAND">1</MDI>\n')
            vrt_file.write('    <MDI key="Y_DATASET">HDF5:"'+(InFileName)+'"://BAND'+bandId+'_RADIANCE/STANDARD_MODE/GEODATA/latitude</MDI>\n')
            vrt_file.write('    <MDI key="Y_BAND">1</MDI>\n')
            vrt_file.write('    <MDI key="PIXEL_OFFSET">0.5</MDI>\n')
            vrt_file.write('    <MDI key="PIXEL_STEP">1</MDI>\n')
            vrt_file.write('    <MDI key="LINE_OFFSET">0.5</MDI>\n')
            vrt_file.write('    <MDI key="LINE_STEP">1</MDI>\n')
            vrt_file.write('  </Metadata>\n')
            vrt_file.write('  <VRTRasterBand dataType="Float32" band="1">\n')
            vrt_file.write('    <Metadata>\n')
            vrt_file.write('      <MDI key="'+(InFileName[:-3])+'_long_name">S5P - Radiance Values</MDI>\n')
            vrt_file.write('      <MDI key="'+(InFileName[:-3])+'_add_offset">' + str(offset) + ' </MDI>\n')
            vrt_file.write('      <MDI key="'+(InFileName[:-3])+'_scale_factor">' + str(scale) + ' </MDI>\n')
            vrt_file.write('      <MDI key="'+(InFileName[:-3])+'_standard_name">Radiance</MDI>\n')
            vrt_file.write('      <MDI key="'+(InFileName[:-3])+'_units">mol.m-2.nm-1.sr-1.s-1</MDI>\n')
            vrt_file.write('      <MDI key="'+(InFileName[:-3])+'__FillValue">9.96921e+36</MDI>\n')
            vrt_file.write('    </Metadata>\n')
            vrt_file.write('    <SimpleSource>\n')
            vrt_file.write('      <SourceFilename relativeToVRT="1">' + bandFilename + '.mpr</SourceFilename>\n')
            vrt_file.write('      <SourceBand>1</SourceBand>\n')
            vrt_file.write('      <SourceProperties RasterXSize="' + str(width) + '" RasterYSize="' + str(height) + '" DataType="Float32" 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'] = 10 + 40 * (z+1)/zsize

            #conduct gdal warp operation using vrt file created and save file as ilwis mpr
            input_raster = gdal.Open(bandFilename + '.vrt')
            output_raster = (bandFilename+'_Resampled.mpr')
            gdal.Warp(output_raster,input_raster,format='ILWIS',srcSRS='WGS84',outputType=gdal.GDT_Float32,resampleAlg='cubic',geoloc=True,srcNodata='0',dstNodata='-999')
            #os.remove ((WorkDir.replace('\\','/')) + '/' + OutFileName +'.vrt')

            progress['value'] = 10 + 65 * (z+1)/zsize
            
            wavelengthInfo = wv[int(bandId)-1]
            step = (wavelengthInfo[1]-wavelengthInfo[0])/zsize
            wavelength = wavelengthInfo[0] + step * (z + 0.5)
            scale = ephoton(wavelength * 1e-9)

            #final radiometric correction applying ilwis-objects - loaded as plugin under python
            rc = ilwis.RasterCoverage(output_raster)
            print (rc.name())
            if exportFormat == 'ilwis':
                rc2=ilwis.Engine.do('mapcalc','iff((@1 >0),@1*'+str(scale)+'+'+str(offset)+',?)',rc)
                rc2.store(bandFilename +'_Rad.mpr', 'map', 'ilwis3')
                rc2=None
                rc=None
                os.remove ((WorkDir.replace('\\','/')) + '/' + bandFilename +'_Resampled.mpr.aux.xml')
                os.remove ((WorkDir.replace('\\','/')) + '/' + bandFilename +'_Resampled.mpr')
                os.remove ((WorkDir.replace('\\','/')) + '/' + bandFilename +'_Resampled.mp#')
                os.remove ((WorkDir.replace('\\','/')) + '/' + bandFilename + '.vrt')
                os.remove ((WorkDir.replace('\\','/')) + '/' + bandFilename + '.mpr')
                #os.remove ((WorkDir.replace('\\','/')) + '/' + bandFilename + '.mp#')

            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 != -999),@1*'+str(scale)+'+'+str(offset)+',-999)',rc)
                rc2.store(bandFilename +'.tif', "GTiff", "gdal")
                os.remove ((WorkDir.replace('\\','/')) + '/' + bandFilename +'_Resampled.mpr.aux.xml')
                os.remove((WorkDir.replace('\\', '/'))+'/'+ bandFilename +'_Resampled.csy')
                os.remove((WorkDir.replace('\\', '/'))+'/'+ bandFilename +'_Resampled.grf')
                os.remove((WorkDir.replace('\\', '/'))+'/'+ bandFilename +'_Resampled.mpr')
                os.remove((WorkDir.replace('\\', '/'))+'/'+ bandFilename +'_Resampled.mp#')
                os.remove ((WorkDir.replace('\\','/')) + '/' + bandFilename + '.vrt')
                os.remove ((WorkDir.replace('\\','/')) + '/' + bandFilename + '.mpr')
                
            progress['value'] = 100 * (z+1)/zsize
    button0.config(state='normal')
    button1.config(state='normal')

    progress['value'] = 100

    for z in range(0,zsize,step):
        bandFilename = OutFileName + '_Radiance_' + str(z)
        os.remove ((WorkDir.replace('\\','/')) + '/' + bandFilename + '.mp#')



#redundant files deleted in main menu script

def executeInThread(SentDir, InFileName, inputInterval, WorkDir, OutFileName, exportFormat):
    button0.config(state='disabled')
    button1.config(state='disabled')
    thread = Thread(target = execute, args = (SentDir, InFileName, inputInterval, 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 S5P L1B Bands")
app.iconbitmap(sent_path +'\pics\esa.ico')

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

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

labelText1 = StringVar()
labelText1.set ("""S5P Sentinel-5 precursor/TROPOMI Level 1B 'BAND 1-8'
Typical filename: S5P_OFFL_L1B_RA_"BAND"_Start of granule in UTC_End of granule in UTC_Orbit number_*_*_*.nc""")

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)

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

inputInterval = StringVar()
E5 = Entry(app, textvariable=inputInterval, bd =2)
E5.pack()
E5.place(x = 165, y = 110, width=45, height=25)

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

outputDir = StringVar()
E3 = Entry(app, textvariable=outputDir, bd =2)
E3.pack()
E3.place(x = 165, y = 140, 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 = 139) 

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

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

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

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

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

R1.select()
R2.deselect()

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

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 = 210)
app.mainloop()

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


