Minor changes

.gitignore

@@ -0,0 +1,5 @@
+.Rproj.user
+.Rhistory
+.RData
+.Ruserdata
+._*

00_Image_Application.py

@@ -0,0 +1,72 @@
+
+
+def image_MDN_application(image_NC, image_ref, SAVE_FOLDER, sensor, useRatio):
+
+    from MDN import image_estimates, get_tile_data, get_sensor_bands
+    import matplotlib.pyplot as plt
+    from matplotlib.colors import LogNorm
+    from osgeo import gdal 
+    import numpy as np
+    import pandas as pd
+    import numpy as np
+
+    sensor = sensor
+
+    IMAGE_PATH_NC = image_NC
+    IMAGE_PATH_TIF_REF = image_ref
+
+    SAVE_FOLDER = SAVE_FOLDER+'/secchi.tif'
+
+
+    bands, Rrs = get_tile_data(IMAGE_PATH_NC, sensor, allow_neg=False)
+
+    kwargs={'product':'secchi-Maciel', 'use_ratio':useRatio}
+
+    Secchi_Est, idxs2 = image_estimates(Rrs, sensor=sensor,**kwargs)
+
+
+
+
+
+    print(Secchi_Est, type(Secchi_Est), Secchi_Est.shape, idxs2)
+
+    Secchi = Secchi_Est[:,:,slice(None, None, None)]
+
+
+    ## Apply and convert to tif
+
+    img_prediction = Secchi_Est.reshape((Rrs.shape[0], Rrs.shape[1]))
+
+    reference_tif = gdal.Open(IMAGE_PATH_TIF_REF, gdal.GA_ReadOnly)
+
+
+    #Save
+    gt = reference_tif.GetGeoTransform()
+    proj = reference_tif.GetProjection()
+    driver = gdal.GetDriverByName("Gtiff")
+
+    driver.Register()
+
+    outds = driver.Create(SAVE_FOLDER, xsize = img_prediction.shape[1], 
+                        ysize=img_prediction.shape[0], bands=1, eType = gdal.GDT_Float32)
+
+    outds.SetGeoTransform(gt)
+
+    outds.SetProjection(proj)
+
+    outband = outds.GetRasterBand(1)
+
+    outband.WriteArray(img_prediction)
+
+    outband.SetNoDataValue(np.nan)
+
+    outband.FlushCache()
+    outband = None
+
+    outds = None
+
+    reference_tif = None
+
+
+
+

01_MDN ETM+ MonteCarlo.py

@@ -0,0 +1,83 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Thu Jun 23 09:10:06 2022
+
+@author: damac
+"""
+
+from MDN.parameters import get_args
+from MDN.utils import split_data
+from MDN import get_sensor_bands
+from MDN.product_estimation import get_estimates
+import pandas as pd
+import matplotlib.pyplot as plt
+import MDN.plot_utils
+import numpy as np
+from sklearn.model_selection import train_test_split
+from sklearn.metrics import mean_absolute_error, mean_squared_error, mean_absolute_percentage_error
+from sklearn.ensemble import RandomForestRegressor
+from sklearn.metrics import r2_score
+from sklearn.metrics import mean_squared_error
+from sklearn.metrics import mean_absolute_error
+import math  
+import random
+import tensorflow as tf
+
+# Dataset creation and separation into train / test (70% / 30%)
+dataset = pd.read_csv(r"Data//Simulated Dataset/rrs_sim_ETM_v3.csv", sep = ',')
+
+dataset.set_index("local_year_month",  inplace = True)
+
+local_unique = np.unique(dataset.index)
+
+#Removing some datase
+
+
+for i in range(0,4):
+
+    lake_train, lake_test = train_test_split(local_unique, test_size = 0.3, random_state=i)
+
+
+    train = dataset.loc[lake_train]
+    test = dataset.drop(lake_train, errors='ignore')
+
+    X_train = train[["Rrs479", "Rrs561", "Rrs661"]]
+    X_test = test[["Rrs479",  "Rrs561", "Rrs661"]]
+
+    y_train = train[['secchi_m']]
+    y_test = test[['secchi_m']]
+
+
+    ## MDN Training and validation 
+    #get_args pulls default arguments from parameters.py unless alternative values are given
+
+    args = get_args(no_load=True,sensor='ETM',product='secchi_maciel', 
+                n_rounds = 10, benchmark=False, use_ratio = True,  n_mix = 5)
+
+
+
+
+    # get_estimates trains model on x_train/y_train data and returns estimates predicted from x_test
+    tf.random.set_seed(1)
+
+    estimates, est_slice = get_estimates(args, X_train,
+                                     y_train, 
+                                     X_test,
+                                     y_test)
+
+    # We take the median over the number of MDN models trained (args.n_rounds = 10 by default)
+
+    estimates_MDN = np.median(estimates,  0) 
+
+
+    test['predicted']  = estimates_MDN
+
+    a = f"/Volumes/Crucial X8/Projects/Secchi_Disk_Landsat/MonteCarlo/ETM/MDN_run_number_{i}.csv"
+
+    pd.DataFrame(test).to_csv(a)
+
+    print(a)
+
+
+
+

01_MDN OLI MonteCarlo.py

@@ -0,0 +1,87 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Thu Jun 23 09:10:06 2022
+
+@author: damac
+"""
+
+
+from MDN.parameters import get_args
+from MDN.utils import split_data
+from MDN import get_sensor_bands
+from MDN.product_estimation import get_estimates
+import pandas as pd
+import matplotlib.pyplot as plt
+import MDN.plot_utils
+import numpy as np
+from sklearn.model_selection import train_test_split
+from sklearn.metrics import mean_absolute_error, mean_squared_error, mean_absolute_percentage_error
+from sklearn.ensemble import RandomForestRegressor
+from sklearn.metrics import r2_score
+from sklearn.metrics import mean_squared_error
+from sklearn.metrics import mean_absolute_error
+import math  
+import random
+import tensorflow as tf
+
+# Dataset creation and separation into train / test (70% / 30%)
+dataset = pd.read_csv(r"Data/Simulated Dataset/rrs_sim_OLI_v3.csv", sep = ',')
+
+dataset.set_index("local_year_month",  inplace = True)
+
+local_unique = np.unique(dataset.index)
+
+
+print(dataset)
+#Removing some datase
+
+
+for i in range(0,50):
+
+    lake_train, lake_test = train_test_split(local_unique, test_size = 0.3, random_state=i)
+
+
+    train = dataset.loc[lake_train]
+    test = dataset.drop(lake_train, errors='ignore')
+
+    X_train = train[["Rrs483", "Rrs561", "Rrs655"]]
+    X_test = test[["Rrs483",  "Rrs561", "Rrs655"]]
+
+    y_train = train[['secchi_m']]
+    y_test = test[['secchi_m']]
+
+
+    ## MDN Training and validation 
+    #get_args pulls default arguments from parameters.py unless alternative values are given
+
+    args = get_args(no_load=True,sensor='OLI-NoCoastal',product='secchi_maciel', 
+                n_rounds = 10, benchmark=False, use_ratio = True,  n_mix = 5)
+
+
+
+
+    # get_estimates trains model on x_train/y_train data and returns estimates predicted from x_test
+    tf.random.set_seed(1)
+
+    estimates, est_slice = get_estimates(args, X_train,
+                                     y_train, 
+                                     X_test,
+                                     y_test)
+
+    # We take the median over the number of MDN models trained (args.n_rounds = 10 by default)
+
+    estimates_MDN = np.median(estimates,  0) 
+
+
+    test['predicted']  = estimates_MDN
+
+    #a = f"MonteCarlo/OLI/MDN_run_number_{i}.csv"
+    a = f"/Volumes/Crucial X8/Projects/Secchi_Disk_Landsat/MonteCarlo/OLI/MDN_run_number_{i}.csv"
+
+    pd.DataFrame(test).to_csv(a)
+
+    print(a)
+
+
+
+

01_MDN TM MonteCarlo.py

@@ -0,0 +1,77 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Thu Jun 23 09:10:06 2022
+@author: damac
+"""
+
+from MDN.parameters import get_args
+from MDN.utils import split_data
+from MDN import get_sensor_bands
+from MDN.product_estimation import get_estimates
+import pandas as pd
+import matplotlib.pyplot as plt
+import MDN.plot_utils
+import numpy as np
+from sklearn.model_selection import train_test_split
+from sklearn.metrics import mean_absolute_error, mean_squared_error, mean_absolute_percentage_error
+from sklearn.ensemble import RandomForestRegressor
+from sklearn.metrics import r2_score
+from sklearn.metrics import mean_squared_error
+from sklearn.metrics import mean_absolute_error
+import math  
+import random
+import tensorflow as tf
+
+# Dataset creation and separation into train / test (70% / 30%)
+dataset = pd.read_csv(r"Data//Simulated Dataset/rrs_sim_TM_v3.csv", sep = ',')
+
+dataset.set_index("local_year_month",  inplace = True)
+local_unique = np.unique(dataset.index)
+
+#Removing some datase
+
+
+for i in range(1,50):
+
+    lake_train, lake_test = train_test_split(local_unique, test_size = 0.3, random_state=i)
+
+
+    train = dataset.loc[lake_train]
+    test = dataset.drop(lake_train, errors='ignore')
+
+    X_train = train[["Rrs486", "Rrs571", "Rrs660"]]
+    X_test = test[["Rrs486",  "Rrs571", "Rrs660"]]
+
+    y_train = train[['secchi_m']]
+    y_test = test[['secchi_m']]
+
+
+    ## MDN Training and validation 
+    #get_args pulls default arguments from parameters.py unless alternative values are given
+
+    args = get_args(no_load=True,sensor='TM',product='secchi_maciel', 
+                n_rounds = 10, benchmark=False, use_ratio = True,  n_mix = 5)
+
+    # get_estimates trains model on x_train/y_train data and returns estimates predicted from x_test
+    tf.random.set_seed(1)
+
+    estimates, est_slice = get_estimates(args, X_train,
+                                     y_train, 
+                                     X_test,
+                                     y_test)
+
+    # We take the median over the number of MDN models trained (args.n_rounds = 10 by default)
+
+    estimates_MDN = np.median(estimates,  0)
+
+
+    test['predicted']  = estimates_MDN
+
+    a = f"MonteCarlo//TM/MDN_run_number_{i}.csv"
+    pd.DataFrame(test).to_csv(a)
+
+    print(a)
+
+
+
+

02_MDN ETM Complete.py

@@ -0,0 +1,57 @@
+# -*- coding: utf-8 -*-
+"""
+@author: Daniel Maciel
+
+This script will apply the MDN to Landsat data. It was aplpying for all data in the .CSV file but it will be filtered after that. 
+"""
+
+from MDN.parameters import get_args
+from MDN.product_estimation import get_estimates
+import pandas as pd
+import matplotlib.pyplot as plt
+import MDN.plot_utils
+import numpy as np
+from sklearn.ensemble import RandomForestRegressor
+from sklearn.metrics import mean_absolute_error
+import math  
+import random
+import tensorflow as tf
+import math  
+
+
+#Data for model creation
+dataset = pd.read_csv(r"Data/Simulated Dataset/rrs_sim_ETM_v3.csv", sep = ',')
+
+#Validation datase
+
+
+BANDS = ('blue', 'green', 'red')
+
+dataset.rename(columns = {'B1':'blue',  "B2":'green', "B3":'red', "B4":'nir'}, inplace = True)
+
+X_train =  dataset.loc[:,BANDS]
+
+print(X_train,  'Train Data')
+
+#Select only secchi for train Y
+y_train = dataset[['secchi']]
+
+print(y_train, 'Secchi Train Data')
+
+#Generate the MDN model
+
+args = get_args(no_load=False,sensor='ETM',product='secchi-Maciel', 
+                n_rounds = 10,
+                benchmark=False, 
+                use_ratio = True)
+
+
+#Saving the MDN model
+tf.random.set_seed(1)
+
+get_estimates(args, X_train, y_train)
+
+
+
+
+