update with some aggregate plotting

README.md

@@ -5,12 +5,12 @@ collected by the NEON AOP. This code was created as part of an effort to generat
 
 A full description of the effort can be found at:
 
-> K. Dana Chadwick, Philip Brodrick, Kathleen Grant, Tristan Goulden, Amanda Henderson, Nicola Falco, Haruko Wainwright, Kenneth H. Williams, Markus Bill, Ian Breckheimer, Eoin L. Brodie, Heidi Steltzer, C. F. Rick Williams, Benjamin Blonder, Jiancong Chen, Baptiste Dafflon, Joan Damerow, Matt Hancher, Aizah Khurram, Jack Lamb, Corey Lawrence, Maeve McCormick. John Musinsky, Samuel Pierce, Alexander Polussa, Maceo Hastings Porro, Andea Scott, Hans Wu Singh, Patrick O. Sorensen, Charuleka Varadharajan, Bizuayehu Whitney, Katharine Maher. Integrating airborne remote sensing and field campaigns for ecology and Earth system science. <i>In Review</i>, 2020.
+> K. Dana Chadwick, Philip Brodrick, Kathleen Grant, Tristan Goulden, Amanda Henderson, Nicola Falco, Haruko Wainwright, Kenneth H. Williams, Markus Bill, Ian Breckheimer, Eoin L. Brodie, Heidi Steltzer, C. F. Rick Williams, Benjamin Blonder, Jiancong Chen, Baptiste Dafflon, Joan Damerow, Matt Hancher, Aizah Khurram, Jack Lamb, Corey Lawrence, Maeve McCormick. John Musinsky, Samuel Pierce, Alexander Polussa, Maceo Hastings Porro, Andea Scott, Hans Wu Singh, Patrick O. Sorensen, Charuleka Varadharajan, Bizuayehu Whitney, Katharine Maher. Integrating airborne remote sensing and field campaigns for ecology and Earth system science. Methods in Ecology and Evolution, 2020.
 
 and use of this code should cite that manuscript.
 
 ### Visualization code in GEE for all products in this project can be found here: 
-https://code.earthengine.google.com/5c96bbc96ffd50e3c8b1433b34a0bb86
+https://code.earthengine.google.com/?scriptPath=users%2Fpgbrodrick%2Feast_river%3Aneon_aop_collection_visuals
 <br>
 
 

apply_cover_model.py

@@ -2,82 +2,108 @@
 import gdal
 import os
 import argparse
-import keras
+from tensorflow import keras
 from tqdm import tqdm
-from sklearn.externals import joblib
-
+#from sklearn.externals import joblib
+import joblib
+#check order of bn & standarization
 
 import warnings
 warnings.filterwarnings('ignore')
 
-parser = argparse.ArgumentParser(description='Apply chem equation to BIL')
-parser.add_argument('refl_dat_f')
-parser.add_argument('output_name')
-parser.add_argument('-model_f', default='trained_models/conifer_nn_set_29.h5')
-parser.add_argument('-scaler', default='trained_models/nn_conifer_scaler')
-parser.add_argument('-bn', default=True, type=bool)
-args = parser.parse_args()
-
-if (os.path.isfile(args.output_name)):
-    print('output file: {} found.  terminating'.format(args.output_name))
-    quit()
-
-# make sure these match the settings file corresponding to the coefficient file
-
-# open up chemical equation data
-model = keras.models.load_model(args.model_f)
-
-# open up raster sets
-dataset = gdal.Open(args.refl_dat_f, gdal.GA_ReadOnly)
-data_trans = dataset.GetGeoTransform()
-
-max_y = dataset.RasterYSize
-max_x = dataset.RasterXSize
-
-
-# create blank output file
-driver = gdal.GetDriverByName('GTiff')
-driver.Register()
-
-
-outDataset = driver.Create(args.output_name,
-                           max_x,
-                           max_y,
-                           1,
-                           gdal.GDT_Float32,
-                           options=['COMPRESS=DEFLATE', 'TILED=YES'])
-
-outDataset.SetProjection(dataset.GetProjection())
-outDataset.SetGeoTransform(dataset.GetGeoTransform())
-
-full_bad_bands = np.zeros(426).astype(bool)
-full_bad_bands[:8] = True
-full_bad_bands[192:205] = True
-full_bad_bands[284:327] = True
-full_bad_bands[417:] = True
-
-bad_bands = np.zeros(142).astype(bool)
-bad_bands[:2] = True
-bad_bands[64:68] = True
-bad_bands[95:109] = True
-bad_bands[139:] = True
-output_predictions = np.zeros((max_y, max_x))
-
-
-if args.scaler is not None:
-    scaler = joblib.load(args.scaler)
-
-# loop through lines [y]
-for l in tqdm(range(0, max_y), ncols=80):
-    dat = np.squeeze(dataset.ReadAsArray(0, l, max_x, 1)).astype(np.float32)
-    dat = dat[np.logical_not(full_bad_bands), ...]
-    dat = np.transpose(dat)
-    if (np.nansum(dat) > 0):
-        if (args.bn):
-            dat = dat / np.sqrt(np.nanmean(np.power(dat, 2), axis=1))[:, np.newaxis]
-
-        dat = scaler.transform(dat)
-        output_predictions[l, :] = model.predict(dat)[:, 0]
-
-outDataset.GetRasterBand(1).WriteArray(output_predictions, 0, 0)
-del outDataset
+def main():
+    parser = argparse.ArgumentParser(description='Apply chem equation to BIL')
+    parser.add_argument('refl_dat_f')
+    parser.add_argument('output_name')
+    parser.add_argument('-model_f', default='trained_models/cover_model_nl_2_dr_0.4_nn_550_it_49.h5')
+    parser.add_argument('-scaler', default=None)
+    parser.add_argument('-bn', default=1, type=int, choices=[0, 1])
+    parser.add_argument('-argmax', default=1, type=int, choices=[0, 1])
+    args = parser.parse_args()
+
+    args.bn = args.bn == 1
+    args.argmax = args.argmax == 1
+
+    if os.path.isfile(args.output_name):
+        print('output file: {} found.  terminating'.format(args.output_name))
+        quit()
+
+    apply_model(args.model_f, args.refl_dat_f, args.output_name, args.bn, args.scaler, args.argmax)
+
+
+def apply_model(model_file: str, refl_dat_f: str, output_name: str, bn: bool = True, scaler: str = None, argmax: bool = True):
+    # open up chemical equation data
+    model = keras.models.load_model(model_file)
+
+    # open up raster sets
+    dataset = gdal.Open(refl_dat_f, gdal.GA_ReadOnly)
+    data_trans = dataset.GetGeoTransform()
+
+    max_y = dataset.RasterYSize
+    max_x = dataset.RasterXSize
+
+
+    # create blank output file
+    driver = gdal.GetDriverByName('GTiff')
+    driver.Register()
+
+
+    full_bad_bands = np.zeros(426).astype(bool)
+    full_bad_bands[:8] = True
+    full_bad_bands[192:205] = True
+    full_bad_bands[284:327] = True
+    full_bad_bands[417:] = True
+
+    bad_bands = np.zeros(142).astype(bool)
+    bad_bands[:2] = True
+    bad_bands[64:68] = True
+    bad_bands[95:109] = True
+    bad_bands[139:] = True
+
+
+    if scaler is not None:
+        scaler = joblib.load(scaler)
+
+    n_bands = 1
+    if argmax is False:
+        dat_bands = np.squeeze(dataset.ReadAsArray(0, 0, 5, 1)).astype(np.float32)
+        dat_bands = dat_bands[np.logical_not(full_bad_bands), ...]
+        dat_bands = np.transpose(dat_bands)
+        pred_shape = model.predict(dat_bands)
+        n_bands = len(pred_shape[-1])
+
+    output_predictions = np.zeros((max_y, max_x, n_bands))
+    outDataset = driver.Create(output_name,
+                               max_x,
+                               max_y,
+                               n_bands,
+                               gdal.GDT_Float32,
+                               options=['COMPRESS=DEFLATE', 'TILED=YES'])
+
+    outDataset.SetProjection(dataset.GetProjection())
+    outDataset.SetGeoTransform(dataset.GetGeoTransform())
+
+    # loop through lines [y]
+    #for l in tqdm(range(0, max_y), ncols=80):
+    for l in tqdm(range(0, 200), ncols=80):
+        dat = np.squeeze(dataset.ReadAsArray(0, l, max_x, 1)).astype(np.float32)
+        dat = dat[np.logical_not(full_bad_bands), ...]
+        dat = np.transpose(dat)
+        if np.nansum(dat) > 0:
+            if bn:
+                dat = dat / np.sqrt(np.nanmean(np.power(dat, 2), axis=1))[:, np.newaxis]
+
+            if scaler is not None:
+                dat = scaler.transform(dat)
+
+            if argmax is True:
+                output_predictions[l, :, 0] = np.argmax(model.predict(dat), axis=-1)
+            else:
+                output_predictions[l, ...] = model.predict(dat)
+
+    for _band in range(n_bands):
+        outDataset.GetRasterBand(_band + 1).WriteArray(output_predictions[...,_band], 0, 0)
+    del outDataset
+
+if __name__ == "__main__":
+    main()

cover_class_spec_plotting.py

@@ -0,0 +1,101 @@
+
+import pandas as pd
+import numpy as np
+import matplotlib.pyplot as plt
+import os
+import warnings
+warnings.filterwarnings('ignore')
+
+
+plt.rcParams.update({'font.size': 20})
+plt.rcParams['font.family'] = "serif"
+plt.rcParams['font.serif'] = "Times New Roman"
+plt.rcParams['axes.grid'] = True
+plt.rcParams['axes.axisbelow'] = True
+plt.rcParams['axes.labelpad'] = 6
+
+
+def main():
+    file_path = '~/Google Drive File Stream/My Drive/CB_share/NEON/cover_classification/extraction_output/cover_extraction_20201021.csv'
+    wavelengths_file = 'raster_data/neon_wavelengths.txt'
+    spectra, cover_types, wv = spec_cleaning(file_path, wavelengths_file, 'refl_B_')
+    plot_avg_spectra(spectra, cover_types, wv)
+
+
+def find_nearest(array_like, v):
+    index = np.argmin(np.abs(np.array(array_like) - v))
+    return index
+
+def spec_cleaning(file_path, wavelengths_file, band_preface):
+    extract = pd.read_csv(file_path)
+    headerSpec = list(extract)  # Reflectance bands start at 18 (17 in zero base)
+
+    # defining wavelengths
+    wv = np.genfromtxt(wavelengths_file)
+    bad_bands = []
+    good_band_ranges = []
+
+    bad_band_ranges = [[0, 425], [1345, 1410], [1805, 2020], [2470, 2700]]
+    for _bbr in range(len(bad_band_ranges)):
+        bad_band_ranges[_bbr] = [find_nearest(wv, x) for x in bad_band_ranges[_bbr]]
+        if (_bbr > 0):
+            good_band_ranges.append([bad_band_ranges[_bbr-1][1], bad_band_ranges[_bbr][0]])
+
+        for n in range(bad_band_ranges[_bbr][0], bad_band_ranges[_bbr][1]):
+            bad_bands.append(n)
+    bad_bands.append(len(wv)-1)
+
+    good_bands = np.array([x for x in range(0, 426) if x not in bad_bands])
+
+    # first column of reflectance data
+    rfdat = list(extract).index(band_preface + '1')
+
+    all_band_indices = (np.array(good_bands)+rfdat).tolist()
+    all_band_indices.extend((np.array(bad_bands)+rfdat).tolist())
+    all_band_indices = np.sort(all_band_indices)
+
+    spectra = np.array(extract[np.array(headerSpec)[all_band_indices]])
+    spectra[:, bad_bands] = np.nan
+
+    cover_types = extract['covertype']
+
+    return spectra, cover_types, wv
+
+
+def plot_avg_spectra(spectra, cover_types, wv, brightness_normalize=False):
+
+    c_types = cover_types.unique()
+
+    color_sets = ['navy', 'tan', 'forestgreen', 'royalblue', 'gray', 'darkorange', 'black', 'brown', 'purple']
+
+    # Plot the difference between needles and noneedles in reflectance data
+    figure_export_settings = {'dpi': 200, 'bbox_inches': 'tight'}
+
+    fig = plt.figure(figsize=(8, 5), constrained_layout=True)
+    for _c, cover in enumerate(c_types):
+        c_spectra = spectra[cover_types == cover]
+        print(color_sets[_c])
+        if brightness_normalize:
+            scale_factor = 1.
+        else:
+            scale_factor = 100.
+        plt.plot(wv, np.nanmean(c_spectra, axis=0) / scale_factor, c=color_sets[_c], linewidth=2)
+        plt.fill_between(wv, np.nanmean(c_spectra, axis=0) / scale_factor - np.nanstd(c_spectra, axis=0) / scale_factor,
+                         np.nanmean(
+                             c_spectra, axis=0) / scale_factor + np.nanstd(c_spectra, axis=0) / scale_factor, alpha=.15,
+                         facecolor=color_sets[_c])
+
+    plt.legend(c_types, prop={'size':10})
+    plt.ylabel('Reflectance (%)')
+    if brightness_normalize:
+        plt.ylabel('Brightness Norm. Reflectance')
+    else:
+        plt.ylabel('Reflectance (%)')
+    plt.xlabel('Wavelength (nm)')
+
+    plt.savefig(os.path.join('figs', 'class_spectra.png'), **figure_export_settings)
+    del fig
+
+
+if __name__ == "__main__":
+    main()