# coding: utf-8
# The aim of this code is to demonstrate how one might apply an ocean colour algorithm to level 2 reflectances in order to estimate chlorophyll-a concentration.

# Import a number of required modules and functions to aid in the processing
import shutil, os, sys
import argparse
import numpy as np
import matplotlib.pyplot as plt
# snappy modules
import snappy
from snappy import Product, ProductData, ProductIO, ProductUtils, ProgressMonitor


# No we define the input data from level 2 (the example used here shows eddy features in the chl field but you could change this to any other level 2 data file set for a region of your own interest).
# We copy the file to the local /tmp space to avoid overwriting the orignal data.
#example area of the Yellow Sea
indir='./DATA/S3A_OL_2_WFR____20181002T020947_20181002T021247_20181003T112306_0179_036_217_2340_MAR_O_NT_002.SEN3'
COPY_DIR = indir.replace('.SEN3', '_copy.SEN3')
COPY_DIR= './tmp/'+os.path.basename(COPY_DIR)
if os.path.isdir(COPY_DIR):
    print(COPY_DIR+' already exists')
else:
    shutil.copytree(indir, COPY_DIR)

# Here we use the xml file to load all the information about the data so that we can process it with snap but we do not need to load the variables into python memory (hence the commented out section to show that it is not actually required to get the bands if you use the expression correctly).
SOURCE_PRODUCT = ProductIO.readProduct(COPY_DIR+'/xfdumanifest.xml')

#Shen 2010 equations 2,4 and 5
Hchl_expression='(0.74*Oa10_reflectance + 0.26*Oa07_reflectance - Oa08_reflectance)'
newband1 = SOURCE_PRODUCT.addBand('Hchl', Hchl_expression, snappy.ProductData.TYPE_FLOAT32)
newband1.setDescription("Hchl from Shen 2010")

HDelta_expression='(Oa07_reflectance-0.5*(Oa06_reflectance+Oa10_reflectance)'
newband2 = SOURCE_PRODUCT.addBand('HDelta', HDelta_expression, snappy.ProductData.TYPE_FLOAT32)
newband2.setDescription("HDelta from Shen 2010")

SCI_expression='Hchl-HDelta'
newband3 = SOURCE_PRODUCT.addBand('SCI', SCI_expression, snappy.ProductData.TYPE_FLOAT32)
newband3.setDescription("SCI index from shen 2010 MERIS chl algorithm")

# Now we convert SCI to Chl
Chl_expression='179378.0*pow(SCI,2) +92.934*SCI + 0.2736'
newband2 = SOURCE_PRODUCT.addBand('Chl_SCI_spring', Chl_expression)
newband2.setDescription("Trial implimentation of the shen 2010 MERIS chl algorithm")

# We will write the file out with the new band 
product_writer = snappy.ProductIO.getProductWriter("NETCDF4-CF")
SOURCE_PRODUCT.setProductWriter(product_writer)
target_file='./tmp/test_netcdf.nc'
print('about to write')
product_writer.writeProductNodes(SOURCE_PRODUCT, target_file)
print('finished writing')
SOURCE_PRODUCT.closeIO()