diff --git a/aggregated_combined_data.py b/aggregated_combined_data.py index bc14810..ba2f167 100644 --- a/aggregated_combined_data.py +++ b/aggregated_combined_data.py @@ -22,14 +22,15 @@ import matplotlib.pyplot as plt import matplotlib.gridspec as gridspec import os -from isofit.radiative_transfer.modtran import ModtranRT -from isofit.radiative_transfer.six_s import SixSRT +from isofit.radiative_transfer.engines.modtran import ModtranRT +from isofit.radiative_transfer.engines.six_s import SixSRT from isofit.configs import configs import argparse from sklearn.linear_model import LinearRegression from sklearn.ensemble import RandomForestRegressor import sklearn.metrics import ray +from isofit.radiative_transfer.radiative_transfer import confPriority @@ -39,163 +40,195 @@ def d2_subset(data,ranges): a = a[:,ranges[1]] return a +def readModtran(modtran_obj, filename): + try: + solzen = modtran_obj.load_tp6(f"{filename}.tp6") + coszen = np.cos(solzen * np.pi / 180.0) + params = modtran_obj.load_chn(f"{filename}.chn", coszen) + + # Remove thermal terms in VSWIR runs to avoid incorrect usage + if modtran_obj.treat_as_emissive is False: + for key in ["thermal_upwelling", "thermal_downwelling"]: + if key in params: + Logger.debug( + f"Deleting key because treat_as_emissive is False: {key}" + ) + del params[key] + + params["solzen"] = solzen + params["coszen"] = coszen + + return params + except: + return None + +def readSixs(sixs_obj, filename, wl, multipart_transmittance=False, wl_size=0): + + return sixs_obj.parse_file( + file=file, + wl=sixs_obj.wl, + multipart_transmittance=sixs_obj.multipart_transmittance, + wl_size=sixs_obj.wl.size, + ) def main(): # Parse arguments parser = argparse.ArgumentParser(description="built luts for emulation.") parser.add_argument('-config_file', type=str, default='templates/isofit_template.json') - parser.add_argument('-keys', type=str, default=['transm', 'rhoatm', 'sphalb'], nargs='+') + parser.add_argument('-keys', type=str, default=['rhoatm', 'sphalb', 'transm_down_dir', 'transm_down_dif', 'transm_up_dir', 'transm_up_dif' ], nargs='+') parser.add_argument('-munge_dir', type=str, default='munged') + parser.add_argument('-figs_dir', type=str, default=None) + parser.add_argument('-unstruct', type=int, default=0, choices=[0,1]) args = parser.parse_args() np.random.seed(13) - for key_ind, key in enumerate(args.keys): - munge_file = os.path.join(args.munge_dir, key + '.npz') - - if os.path.isfile(munge_file) is False: - config = configs.create_new_config(args.config_file) - - # Note - this goes way faster if you comment out the Vector Interpolater build section in each of these - isofit_modtran = ModtranRT(config.forward_model.radiative_transfer.radiative_transfer_engines[0], - config, build_lut = False) - isofit_sixs = SixSRT(config.forward_model.radiative_transfer.radiative_transfer_engines[1], - config, build_lut = False) - - sixs_results = get_obj_res(isofit_sixs, key, resample=False) - modtran_results = get_obj_res(isofit_modtran, key) - - if os.path.isdir(os.path.dirname(munge_file) is False): - os.mkdir(os.path.dirname(munge_file)) - - for fn in isofit_modtran.files: - mod_output = isofit_modtran.load_rt(fn) - sol_irr = mod_output['sol'] - if np.all(np.isfinite(sol_irr)): - break - - np.savez(munge_file, modtran_results=modtran_results, sixs_results=sixs_results, sol_irr=sol_irr) - - modtran_results = None - sixs_results = None - for key_ind, key in enumerate(args.keys): - munge_file = os.path.join(args.munge_dir, key + '.npz') - - npzf = np.load(munge_file) - - dim1 = int(np.product(np.array(npzf['modtran_results'].shape)[:-1])) - dim2 = npzf['modtran_results'].shape[-1] - if modtran_results is None: - modtran_results = np.zeros((dim1,dim2*len(args.keys))) - modtran_results[:,dim2*key_ind:dim2*(key_ind+1)] = npzf['modtran_results'] - dim1 = int(np.product(np.array(npzf['sixs_results'].shape)[:-1])) - dim2 = npzf['sixs_results'].shape[-1] - if sixs_results is None: - sixs_results = np.zeros((dim1,dim2*len(args.keys))) - sixs_results[:,dim2*key_ind:dim2*(key_ind+1)] = npzf['sixs_results'] - - sol_irr = npzf['sol_irr'] config = configs.create_new_config(args.config_file) - isofit_modtran = ModtranRT(config.forward_model.radiative_transfer.radiative_transfer_engines[0], - config, build_lut=False) - isofit_sixs = SixSRT(config.forward_model.radiative_transfer.radiative_transfer_engines[1], - config, build_lut=False) - sixs_names = isofit_sixs.lut_names - modtran_names = isofit_modtran.lut_names - - if 'elev' in sixs_names: - sixs_names[sixs_names.index('elev')] = 'GNDALT' - if 'viewzen' in sixs_names: - sixs_names[sixs_names.index('viewzen')] = 'OBSZEN' - if 'viewaz' in sixs_names: - sixs_names[sixs_names.index('viewaz')] = 'TRUEAZ' - if 'alt' in sixs_names: - sixs_names[sixs_names.index('alt')] = 'H1ALT' - if 'AOT550' in sixs_names: - sixs_names[sixs_names.index('AOT550')] = 'AERFRAC_2' - - reorder_sixs = [sixs_names.index(x) for x in modtran_names] - - points = isofit_modtran.points.copy() - points_sixs = isofit_sixs.points.copy()[:,reorder_sixs] - - if 'OBSZEN' in modtran_names: - print('adjusting') - points_sixs[:, modtran_names.index('OBSZEN')] = 180 - points_sixs[:, modtran_names.index('OBSZEN')] - - - ind = np.lexsort(tuple([points[:,x] for x in range(points.shape[-1])])) - points = points[ind,:] - modtran_results = modtran_results[ind,:] - ind_sixs = np.lexsort(tuple([points_sixs[:,x] for x in range(points_sixs.shape[-1])])) - points_sixs = points_sixs[ind_sixs,:] - sixs_results = sixs_results[ind_sixs,:] + # Note - this goes way faster if you comment out the Vector Interpolater build section in each of these + rt_config = config.forward_model.radiative_transfer + instrument_config = config.forward_model.instrument - good_data = np.all(np.isnan(modtran_results) == False,axis=1) - good_data[np.any(np.isnan(sixs_results),axis=1)] = False + + params = {'engine_config': rt_config.radiative_transfer_engines[0]} + + params['lut_grid'] = confPriority('lut_grid', [params['engine_config'], instrument_config, rt_config]) + params['lut_grid'] = {key: params['lut_grid'][key] for key in params['engine_config'].lut_names.keys()} + params['wavelength_file'] = confPriority('wavelength_file', [params['engine_config'], instrument_config, rt_config]) + if args.unstruct: + params['engine_config'].rte_configure_and_exit = True + else: + params['engine_config'].rte_configure_and_exit = False + #params['engine_config'].rt_mode = 'rdn' + isofit_modtran = ModtranRT(**params) + + params = {'engine_config' : rt_config.radiative_transfer_engines[1]} + + params['lut_grid'] = confPriority('lut_grid', [params['engine_config'], instrument_config, rt_config]) + params['lut_grid'] = {key: params['lut_grid'][key] for key in params['engine_config'].lut_names.keys()} + if args.unstruct: + params['engine_config'].rte_configure_and_exit = True + else: + params['engine_config'].rte_configure_and_exit = False + #params['engine_config'].rt_mode = 'rdn' - modtran_results = modtran_results[good_data,:] - sixs_results = sixs_results[good_data,:] - points = points[good_data,...] - points_sixs = points_sixs[good_data,...] + # Get raw 6s return, not wavelength convolved (this is what we'll use for inference too) + sixs_wl = np.arange(350, 2500 + 2.5, 2.5) + sixs_fwhm = np.full(sixs_wl.size, 2.0) + #params['wavelength_file'] = confPriority('wavelength_file', [params['engine_config'], instrument_config, rt_config]) + params['wl'] = sixs_wl + params['fwhm'] = sixs_fwhm + isofit_sixs = SixSRT(**params) - print(sixs_results.shape) - print(modtran_results.shape) - tmp = isofit_sixs.load_rt(isofit_sixs.files[0]) - np.savez(os.path.join(args.munge_dir, 'combined_training_data.npz'), sixs_results=sixs_results, modtran_results=modtran_results, - points=points, points_sixs=points_sixs, keys=args.keys, point_names=modtran_names, modtran_wavelengths=isofit_modtran.wl, - sixs_wavelengths=isofit_sixs.grid, - sol_irr=sol_irr) -@ray.remote -def read_data_piece(ind, maxind, point, fn, key, resample, obj): - if ind % 100 == 0: - print('{}: {}/{}'.format(key, ind, maxind)) - try: - if resample is False: - mod_output = obj.load_rt(fn, resample=False) - else: - mod_output = obj.load_rt(fn) - res = mod_output[key] - except: - res = None - return ind, res - + outdict_sixs, outdict_modtran = {}, {} + munge_file = os.path.join(args.munge_dir, 'combined_data.npz') + good_points = np.ones(isofit_sixs.lut['rhoatm'].shape[0]).astype(bool) + for key_ind, key in enumerate(args.keys): -def get_obj_res(obj, key, resample=True): + if args.figs_dir is not None: + point_dims = list(isofit_modtran.lut_grid.keys()) + point_dims_s = list(isofit_sixs.lut_grid.keys()) + rtm_key = 'transm_down_dir' + for _ind in range(isofit_sixs.lut['rhoatm'].shape[0]): + + sp = isofit_sixs.points[_ind,:] + mp = isofit_modtran.points[_ind,:] + + std_dir = isofit_sixs.lut['transm_down_dir'][_ind,:] + mtd_dir = isofit_modtran.lut['transm_down_dir'][_ind,:] + std_dif = isofit_sixs.lut['transm_down_dif'][_ind,:] + mtd_dif = isofit_modtran.lut['transm_down_dif'][_ind,:] + stu_dir = isofit_sixs.lut['transm_up_dir'][_ind,:] + mtu_dir = isofit_modtran.lut['transm_up_dir'][_ind,:] + stu_dif = isofit_sixs.lut['transm_up_dif'][_ind,:] + mtu_dif = isofit_modtran.lut['transm_up_dif'][_ind,:] + s_r = isofit_sixs.lut['rhoatm'][_ind,:] + m_r = isofit_modtran.lut['rhoatm'][_ind,:] + + name='_'.join([f'{point_dims_s[x]}_{sp[x]}' for x in range(len(sp))]) + plt.plot(isofit_sixs.wl, std_dir, color='black', label='t_down_dir') + plt.plot(isofit_sixs.wl, std_dif, color='grey', label='t_down_dif') + plt.plot(isofit_sixs.wl, stu_dir, color='red', label='t_up_dir') + plt.plot(isofit_sixs.wl, stu_dif, color='purple', label='t_up_dif') + plt.plot(isofit_sixs.wl, s_r, color='green', label='rhoatm') + + name2='_'.join([f'{point_dims[x]}_{mp[x]}' for x in range(len(mp))]) + plt.plot(isofit_modtran.wl, mtd_dir, color='black', ls = '--') + plt.plot(isofit_modtran.wl, mtd_dif, color='grey', ls = '--') + plt.plot(isofit_modtran.wl, mtu_dir, color='red', ls = '--') + plt.plot(isofit_modtran.wl, mtu_dif, color='purple', ls = '--') + plt.plot(isofit_modtran.wl, m_r, color='green', ls = '--') + + plt.legend(fontsize=4, loc='lower right') + + plt.title(f'S: {name}\n M: {name2}',fontsize=6) + + plt.savefig(f'{args.figs_dir}/{name.replace("\n","_")}.png',dpi=100) + plt.clf() + + point_names = isofit_sixs.lut.point.to_index().names + bad_points = np.zeros(isofit_sixs.lut[key].shape[0],dtype=bool) + if 'surface_elevation_km' in point_names and 'observer_altitude_km' in point_names: + bad_points = isofit_sixs.lut['surface_elevation_km'] >= isofit_sixs.lut['observer_altitude_km'] -2 + bad_points[np.any(isofit_sixs.lut['transm_down_dif'].data[:,:] > 10,axis=1)] = True + bad_points[np.any(isofit_modtran.lut['transm_down_dif'].data[:,:] > 10,axis=1)] = True + good_points = np.logical_not(bad_points) + + outdict_sixs[key] = np.array(isofit_sixs.lut[key])[good_points,:] + outdict_modtran[key] = np.array(isofit_modtran.lut[key])[good_points,:] + + + # hack at some things to prevent runaway values + outdict_sixs['transm_up_dif'][outdict_sixs['transm_up_dif'][:,:] > 1] = 0 + outdict_sixs['transm_up_dif'][outdict_sixs['transm_up_dif'][:,:] < 0] = 0 + outdict_modtran['transm_up_dif'][outdict_modtran['transm_up_dif'][:,:] > 1] = 0 + outdict_modtran['transm_up_dif'][outdict_modtran['transm_up_dif'][:,:] < 0] = 0 + + outdict_sixs['sphalb'][:,isofit_sixs.wl > 1200][outdict_sixs['sphalb'][:,isofit_sixs.wl > 1200] > 0.1] = 0 + outdict_sixs['sphalb'][outdict_sixs['sphalb'][:,:] < 0] = 0 + + subs = outdict_modtran['sphalb'][:,isofit_modtran.wl > 1200].copy() + subs[subs > 0.1] = 0 + outdict_modtran['sphalb'][:,isofit_modtran.wl > 1200] = subs + print(np.sum(outdict_modtran['sphalb'][:,isofit_modtran.wl > 1200] > 0.1)) + + outdict_modtran['sphalb'][outdict_modtran['sphalb'][:,:] < 0] = 0 + + + + #keys=list(isofit_modtran.lut_grid.keys()) + keys=['rhoatm','sphalb','transm_down_dir','transm_down_dif', 'transm_up_dir','transm_up_dif'] + #keys=list(isofit_sixs.lut.point.to_index().names) + stacked_modtran = np.zeros((outdict_modtran[keys[0]].shape[0], outdict_modtran[keys[0]].shape[1] * len(keys))) + stacked_sixs = np.zeros((outdict_sixs[keys[0]].shape[0], outdict_sixs[keys[0]].shape[1] * len(keys))) + n_bands_modtran = int(stacked_modtran.shape[-1]/len(keys)) + n_bands_sixs = int(stacked_sixs.shape[-1]/len(keys)) + for n in range(len(keys)): + stacked_modtran[:,n*n_bands_modtran:(n+1)*n_bands_modtran] = outdict_modtran[keys[n]] + stacked_sixs[:,n*n_bands_sixs:(n+1)*n_bands_sixs] = outdict_sixs[keys[n]] + + + np.savez(munge_file, modtran_results=stacked_modtran, + sixs_results=stacked_sixs, + points=isofit_sixs.points[good_points,:], + sol_irr=isofit_modtran.lut.solar_irr, + sixs_wavelengths=isofit_sixs.wl, + modtran_wavelengths=isofit_modtran.wl, + point_names=list(isofit_sixs.lut.point.to_index().names), + keys=keys + ) - # We don't want the VectorInterpolator, but rather the raw inputs - ray.init(temp_dir='/tmp/ray/brodrick/') - if hasattr(obj,'sixs_ngrid_init'): - results = np.zeros((obj.points.shape[0],obj.sixs_ngrid_init), dtype=float) - else: - results = np.zeros((obj.points.shape[0],obj.n_chan), dtype=float) - objid = ray.put(obj) - jobs = [] - for ind, (point, fn) in enumerate(zip(obj.points, obj.files)): - jobs.append(read_data_piece.remote(ind, results.shape[0], point, fn, key, resample, objid)) - rreturn = [ray.get(jid) for jid in jobs] - for ind, res in rreturn: - if res is not None: - try: - results[ind,:] = res - except: - results[ind,:] = np.nan - else: - results[ind,:] = np.nan - ray.shutdown() - return results if __name__ == '__main__': diff --git a/build_luts.py b/build_luts.py index 9afd125..4d8ec7f 100644 --- a/build_luts.py +++ b/build_luts.py @@ -18,144 +18,356 @@ import numpy as np -from isofit.utils.apply_oe import write_modtran_template, SerialEncoder import os import json -from isofit.radiative_transfer.modtran import ModtranRT -from isofit.radiative_transfer.six_s import SixSRT -from isofit.configs import configs, Config import datetime import ray import argparse -from isofit.core.sunposition import sunpos +from types import SimpleNamespace +import logging +import subprocess + +from scipy.stats import qmc + +from isofit.utils.template_construction import ( + write_modtran_template, + SerialEncoder, + get_lut_subset, +) +from isofit.data import env +from isofit.radiative_transfer.radiative_transfer import confPriority +from isofit.radiative_transfer.engines.modtran import ModtranRT +from isofit.radiative_transfer.engines.six_s import SixSRT +from isofit.configs import configs, Config def main(): - # Parse arguments parser = argparse.ArgumentParser(description="built luts for emulation.") - parser.add_argument('-ip_head', type=str) - parser.add_argument('-redis_password', type=str) + parser.add_argument('-dir', default='', help='Working Dir') parser.add_argument('-n_cores', type=int, default=1) parser.add_argument('-train', type=int, default=1, choices=[0,1]) parser.add_argument('-cleanup', type=int, default=0, choices=[0,1]) - + parser.add_argument('-ip_head', type=str) + parser.add_argument('-redis_password', type=str) + parser.add_argument( + '--configure_and_exit', + default=False, + action=argparse.BooleanOptionalAction + ) + parser.add_argument( + '--sobel', + default=False, + action=argparse.BooleanOptionalAction + ) + parser.add_argument( + '-sixs_path', + default=None, + help='SIXS Installation DIR' + ) + parser.add_argument( + '-modtran_path', + default=None, + help='MODTRAN Installation DIR' + ) + parser.add_argument( + '--coarse', + default=None, + type=str + ) args = parser.parse_args() - args.train = args.train == 1 + args.training = args.train == 1 args.cleanup = args.cleanup == 1 dayofyear = 200 - if args.train: - to_solar_zenith_lut = [0, 12.5, 25, 37.5, 50] - to_solar_azimuth_lut = [180] - to_sensor_azimuth_lut = [180] - to_sensor_zenith_lut = [140, 160, 180] - altitude_km_lut = [2, 4, 7, 10, 15, 25] - elevation_km_lut = [0, 0.75, 1.5, 2.25, 4.5] - h2o_lut_grid = np.round(np.linspace(0.1,5,num=5),3).tolist() + [0.6125] - h2o_lut_grid.sort() - aerfrac_2_lut_grid = np.round(np.linspace(0.01,1,num=5),3).tolist() + [0.5] - aerfrac_2_lut_grid.sort() + paths = Paths( + args, + args.training + ) + + if args.sobel: + if args.train: + # Get bounds from regular grid + bounds = { + 'to_solar_zenith_bnds': [0, 50], + 'to_sensor_zenith_bnds': [140, 180], + 'altitude_km_bnds': [2, 25], + 'elevation_km_bnds': [0.01, 4.5], + 'h2o_bnds': [0.1, 5], + 'aerfrac_2_bnds': [0.01, 1], + } + # consts won't be sampled from the grid + consts = { + 'to_solar_azimuth_bnds': 180, + 'to_sensor_azimuth_bnds': 180, + } + # n is the number of samples to take from the grid. + grid = sobel_lut(bounds, consts, n=3) + + # Add in relative azimuth + relative_azimuth = np.abs( + np.array(grid['to_solar_azimuth_lut']) + - np.array(grid['to_sensor_azimuth_lut']) + ) + grid['relative_azimuth_lut'] = [float(i) for i in np.minimum( + relative_azimuth, + 360 - relative_azimuth + )] + else: + # Get bounds from regular grid + bounds = { + 'to_solar_zenith_bnds': [0, 50], + 'to_sensor_zenith_bnds': [140, 180], + 'altitude_km_bnds': [2, 25], + 'elevation_km_bnds': [0.01, 4.5], + 'h2o_bnds': [0.1, 5], + 'aerfrac_2_bnds': [0.01, 1], + } + consts = { + 'to_solar_azimuth_bnds': 180, + 'to_sensor_azimuth_bnds': 180, + } + grid = sobel_lut(bounds, consts, n=10) + + # Add in relative azimuth + relative_azimuth = np.abs( + np.array(grid['to_solar_azimuth_lut']) + - np.array(grid['to_sensor_azimuth_lut']) + ) + grid['relative_azimuth_lut'] = [float(i) for i in np.minimum( + relative_azimuth, + 360 - relative_azimuth + )] else: - # HOLDOUT SET - to_solar_zenith_lut = [6, 18, 30,45] - to_solar_azimuth_lut = [60, 300] - to_sensor_azimuth_lut = [180] - to_sensor_zenith_lut = [145, 155, 165, 175] - altitude_km_lut = [3, 5.5, 8.5, 12.5, 17.5] - elevation_km_lut = [0.325, 1.025, 1.875, 2.575, 4.2] - h2o_lut_grid = np.round(np.linspace(0.5,4.5,num=4),3) - aerfrac_2_lut_grid = np.round(np.linspace(0.125,0.9,num=4),3) - - - - n_lut_build = np.product([len(to_solar_zenith_lut), - len(to_solar_azimuth_lut), - len(to_sensor_zenith_lut), - len(to_sensor_azimuth_lut), - len(altitude_km_lut), - len(elevation_km_lut), - len(h2o_lut_grid), - len(aerfrac_2_lut_grid)]) + if args.train: + grid = { + #'to_solar_zenith_lut': [0, 12.5, 25, 37.5, 50], + 'to_solar_zenith_lut': [0, 30, 60], + #'to_solar_azimuth_lut': [0, 60, 120, 180], + 'to_solar_azimuth_lut': [0, 90, 180], + 'to_sensor_azimuth_lut': [0], + 'to_sensor_zenith_lut': [140, 160, 180], + 'altitude_km_lut': [2, 4, 7, 10, 15, 25, 99], + #'elevation_km_lut': [0, 0.75, 1.5, 2.25, 4.5, 6], + 'elevation_km_lut': [0, 1.5, 4.5, 6], + 'h2o_lut': list(np.sort( + np.round(np.linspace(0.1, 5, num=5), 3).tolist() + )), #+ [0.7125, 1.9375, 3.162, 4.3875], + 'aerfrac_2_lut': list(np.sort( + np.round(np.linspace(0.01, 1, num=5), 3).tolist() + )) + #'to_solar_zenith_lut': [0, 50], + #'to_solar_azimuth_lut': [0, 180], + #'to_sensor_azimuth_lut': [0], + #'to_sensor_zenith_lut': [140, 180], + #'altitude_km_lut': [2, 99], + #'elevation_km_lut': [0, 6], + #'h2o_lut': list(np.sort( + # np.round(np.linspace(0.1, 5, num=2), 3).tolist() + [0.6125] + #)), + #'aerfrac_2_lut': list(np.sort( + # np.round(np.linspace(0.01, 1, num=5), 3).tolist() + [0.5] + #)) + } + relative_azimuth_lut = np.abs( + np.array(grid['to_solar_azimuth_lut']) + - np.array(grid['to_sensor_azimuth_lut']) + ) + grid['relative_azimuth_lut'] = list(np.minimum( + relative_azimuth_lut, + 360 - relative_azimuth_lut + )) + + else: + # HOLDOUT SET + grid = { + 'to_solar_zenith_lut': [6, 18, 30,45], + 'to_solar_azimuth_lut': [135], + 'to_sensor_azimuth_lut': [0], + 'to_sensor_zenith_lut': [145, 155, 165, 175], + 'altitude_km_lut': [3, 5.5, 8.5, 12.5, 17.5], + 'elevation_km_lut': [0.325, 1.025, 1.875, 2.575, 4.2], + 'h2o_lut': np.round(np.linspace(0.5, 4.5, num=4), 3), + 'aerfrac_2_lut': np.round(np.linspace(0.125, 0.9, num=4), 3), + } + relative_azimuth_lut = np.abs( + np.array(grid['to_solar_azimuth_lut']) + - np.array(grid['to_sensor_azimuth_lut']) + ) + grid['relative_azimuth_lut'] = np.minimum( + relative_azimuth_lut, + 360 - relative_azimuth_lut + ) + + n_lut_build = np.prod([ + len(grid['to_solar_zenith_lut']), + len(grid['to_sensor_zenith_lut']), + len(grid['relative_azimuth_lut']), + len(grid['altitude_km_lut']), + len(grid['elevation_km_lut']), + len(grid['h2o_lut']), + len(grid['aerfrac_2_lut']) + ]) print('Num LUTs to build: {}'.format(n_lut_build)) print('Expected MODTRAN runtime: {} hrs'.format(n_lut_build*1.5)) print('Expected MODTRAN runtime: {} days'.format(n_lut_build*1.5/24)) print('Expected MODTRAN runtime per (40-core) node: {} days'.format(n_lut_build*1.5/24/40)) + # Create wavelength file - wl = np.arange(0.350, 2.550, 0.0005) - wl_file_contents = np.zeros((len(wl),3)) - wl_file_contents[:,0] = np.arange(len(wl),dtype=int) - wl_file_contents[:,1] = wl - wl_file_contents[:,2] = 0.0005 - np.savetxt('support/hifidelity_wavelengths.txt',wl_file_contents,fmt='%.5f') - - # Initialize ray for parallel execution - rayargs = {'address': args.ip_head, - 'redis_password': args.redis_password, - 'local_mode': args.n_cores == 1} - - if args.n_cores < 40: - rayargs['num_cpus'] = args.n_cores - ray.init(**rayargs) - print(ray.cluster_resources()) - - template_dir = 'templates' - if os.path.isdir(template_dir) is False: - os.mkdir(template_dir) - - modtran_template_file = os.path.join(template_dir,'modtran_template.json') - - if args.training: - isofit_config_file = os.path.join(template_dir,'isofit_template_v2.json') - else: - isofit_config_file = os.path.join(template_dir,'isofit_template_holdout.json') - - write_modtran_template(atmosphere_type='ATM_MIDLAT_SUMMER', fid=os.path.splitext(modtran_template_file)[0], - altitude_km=altitude_km_lut[0], - dayofyear=dayofyear, latitude=to_solar_azimuth_lut[0], longitude=to_solar_zenith_lut[0], - to_sensor_azimuth=to_sensor_azimuth_lut[0], to_sensor_zenith=to_sensor_zenith_lut[0], - gmtime=0, elevation_km=elevation_km_lut[0], output_file=modtran_template_file) - + if args.coarse is None: + wl = np.arange(0.350, 2.550, 0.0001) + wl_file_contents = np.zeros((len(wl),3)) + wl_file_contents[:,0] = np.arange(len(wl),dtype=int) + wl_file_contents[:,1] = wl + wl_file_contents[:,2] = 0.0001 + + np.savetxt( + paths.wavelength_file, + wl_file_contents,fmt='%.5f' + ) + + write_modtran_template( + atmosphere_type='ATM_MIDLAT_SUMMER', + fid=os.path.splitext(paths.modtran_template_file)[0], + altitude_km=grid['altitude_km_lut'][0], + dayofyear=dayofyear, + to_sensor_azimuth=grid['to_sensor_azimuth_lut'][0], + to_sensor_zenith=grid['to_sensor_zenith_lut'][0], + to_sun_zenith=grid['to_solar_zenith_lut'][0], + relative_azimuth=grid['relative_azimuth_lut'][0], + gmtime=0, + elevation_km=grid['elevation_km_lut'][0], + output_file=paths.modtran_template_file, + ihaze_type="AER_NONE", + ) # Make sure H2O grid is fully valid - with open(modtran_template_file, 'r') as f: + with open(paths.modtran_template_file, 'r') as f: modtran_config = json.load(f) - modtran_config['MODTRAN'][0]['MODTRANINPUT']['GEOMETRY']['IPARM'] = 12 - modtran_config['MODTRAN'][0]['MODTRANINPUT']['ATMOSPHERE']['H2OOPT'] = '+' - modtran_config['MODTRAN'][0]['MODTRANINPUT']['AEROSOLS']['VIS'] = 100 - with open(modtran_template_file, 'w') as fout: - fout.write(json.dumps(modtran_config, cls=SerialEncoder, indent=4, sort_keys=True)) - - paths = Paths(os.path.join('.',os.path.basename(modtran_template_file)), args.training) - - - build_main_config(paths, isofit_config_file, to_solar_azimuth_lut, to_solar_zenith_lut, - aerfrac_2_lut_grid, h2o_lut_grid, elevation_km_lut, altitude_km_lut, to_sensor_azimuth_lut, - to_sensor_zenith_lut, n_cores=args.n_cores) - config = configs.create_new_config(isofit_config_file) + for n in range(len(modtran_config['MODTRAN'])): + modtran_config['MODTRAN'][n]['MODTRANINPUT']['GEOMETRY']['IPARM'] = 12 + modtran_config['MODTRAN'][n]['MODTRANINPUT']['ATMOSPHERE']['H2OOPT'] = '+' + modtran_config['MODTRAN'][n]['MODTRANINPUT']['AEROSOLS']['VIS'] = 100 + if args.coarse is not None: + modtran_config['MODTRAN'][n]['MODTRANINPUT']['SPECTRAL']['BMNAME'] = '05_2013' + modtran_config['MODTRAN'][n]['MODTRANINPUT']['SPECTRAL']['DV'] = 5 + modtran_config['MODTRAN'][n]['MODTRANINPUT']['SPECTRAL']['FWHM'] = 5 + + with open(paths.modtran_template_file, 'w') as fout: + fout.write(json.dumps( + modtran_config, + cls=SerialEncoder, + indent=4, + sort_keys=True + )) + + configure_and_exit = args.configure_and_exit + build_main_config( + paths, + grid['h2o_lut'], + grid['elevation_km_lut'], + grid['altitude_km_lut'], + grid['to_sensor_zenith_lut'], + grid['to_solar_zenith_lut'], + grid['to_sensor_azimuth_lut'], + grid['to_solar_azimuth_lut'], + grid['relative_azimuth_lut'], + grid['aerfrac_2_lut'], + paths.isofit_config_file, + configure_and_exit=configure_and_exit, + ip_head=args.ip_head, + redis_password=args.redis_password, + n_cores=args.n_cores + ) + + config = configs.create_new_config(paths.isofit_config_file) + + if args.configure_and_exit: + ray.init( + num_cpus=args.n_cores, + _temp_dir=config.implementation.ray_temp_dir , + include_dashboard=False, + local_mode=args.n_cores == 1, + ) + else: + # Initialize ray for parallel execution + rayargs = { + 'address': config.implementation.ip_head, + '_redis_password': config.implementation.redis_password, + 'local_mode': args.n_cores == 1 + } + + if args.n_cores < 40: + rayargs['num_cpus'] = args.n_cores + + ray.init(**rayargs) + print(ray.cluster_resources()) + + # Inits of engines based on radiative_transfer.py + rt_config = config.forward_model.radiative_transfer + instrument_config = config.forward_model.instrument + lut_grid_config = rt_config.lut_grid + + _keys = [ + "interpolator_style", + "overwrite_interpolator", + "lut_grid", + "lut_path", + "wavelength_file", + ] + + modtran_engine_config = rt_config.radiative_transfer_engines[0] + params = { + key: confPriority(key, [ + modtran_engine_config, instrument_config, rt_config + ]) + for key in _keys + } + params["engine_config"] = modtran_engine_config + params["build_interpolators"] = False + params['lut_grid'] = { + key: params['lut_grid'][key] for key in + modtran_engine_config.lut_names.keys() + } + isofit_modtran = ModtranRT(**params) + + sixs_engine_config = rt_config.radiative_transfer_engines[1] + params = { + key: confPriority(key, [ + sixs_engine_config, instrument_config, rt_config + ]) + for key in _keys + } + params["engine_config"] = sixs_engine_config + params["build_interpolators"] = False + params['lut_grid'] = { + key: params['lut_grid'][key] for key in + sixs_engine_config.lut_names.keys() + } + # params['modtran_emulation'] = True - isofit_modtran = ModtranRT(config.forward_model.radiative_transfer.radiative_transfer_engines[0], - config) - - isofit_sixs = SixSRT(config.forward_model.radiative_transfer.radiative_transfer_engines[1], - config) + # always overwrite 6S wavelength with full raw 6S range + sixs_wl = np.arange(350, 2500 + 2.5, 2.5) + sixs_fwhm = np.full(sixs_wl.size, 2.0) + params['wl'] = sixs_wl + params['fwhm'] = sixs_fwhm + isofit_sixs = SixSRT(**params) # cleanup if args.cleanup: for to_rm in ['*r_k', '*t_k', '*tp7', '*wrn', '*psc', '*plt', '*7sc', '*acd']: - cmd = 'find {os.path.join(paths.lut_modtran_directory)} -name "{to_rm}"') + cmd = 'find {os.path.join(paths.lut_modtran_directory)} -name "{to_rm}"' print(cmd) os.system(cmd) - def build_modtran_configs(isofit_config: Config, template_file: str): with open(template_file, 'r') as f: modtran_config = json.load(f) @@ -163,60 +375,159 @@ def build_modtran_configs(isofit_config: Config, template_file: str): class Paths(): - def __init__(self, modtran_tpl, training=True): - self.aerosol_tpl_path = '../support/aerosol_template.json' - self.aerosol_model_path = '../support/aerosol_model.txt' - self.wavelenth_file = '../support/hifidelity_wavelengths.txt' - self.earth_sun_distance_file = '../support/earth_sun_distance.txt' - self.irradiance_file = '../support/prism_optimized_irr.dat' + def __init__(self, args, training=True): - if args.training: - self.lut_modtran_directory = '../modtran_lut' - self.lut_sixs_directory = '../sixs_lut' - else: - self.lut_modtran_directory = '../modtran_lut_holdout_az' - self.lut_sixs_directory = '../sixs_lut_holdout_az' + working_dir = args.dir + + self.support_dir = os.path.join(working_dir, 'support') + self.template_dir = os.path.join(working_dir, 'templates') - self.modtran_template_path = modtran_tpl + # Make dirs + os.makedirs(self.support_dir, exist_ok=True) + os.makedirs(self.template_dir, exist_ok=True) + self.modtran_template_file = os.path.join( + self.template_dir, 'modtran_template.json' + ) + + if training: + self.isofit_config_file = os.path.join( + self.template_dir, 'isofit_template_v2.json' + ) + else: + self.isofit_config_file = os.path.join( + self.template_dir, 'isofit_template_holdout.json' + ) + + self.aerosol_tpl_path = os.path.join( + self.support_dir, + 'aerosol_template.json' + ) + self.aerosol_model_path = os.path.join( + self.support_dir, + 'aerosol_model.txt' + ) + + if args.coarse: + self.wavelength_file = os.path.join( + self.support_dir, + 'coarse_wavelengths.txt' + ) + if os.path.isfile(args.coarse) is False: + print(f'No wavelength file found for args.coarse at {args.coarse}') + exit() + subprocess.call(f'cp {args.coarse} {self.wavelength_file}',shell=True) + else: + self.wavelength_file = os.path.join( + self.support_dir, + 'hifidelity_wavelengths.txt' + ) + + # If we want to format this with the data repo + # self.earth_sun_distance_file = str(env.path("data", "emit_noise.txt") + self.earth_sun_distance_file = os.path.join( + self.support_dir, + 'earth_sun_distance.txt' + ) + self.irradiance_file = os.path.join( + self.support_dir, + 'kurudz_0.1nm.dat' + ) + + if training: + self.lut_modtran_directory = os.path.join( + working_dir, + 'modtran_lut' + ) + self.lut_sixs_directory = os.path.join( + working_dir, + 'sixs_lut' + ) + else: + self.lut_modtran_directory = os.path.join( + working_dir, + 'modtran_lut_holdout_az' + ) + self.lut_sixs_directory = os.path.join( + working_dir, + 'sixs_lut_holdout_az' + ) + + os.makedirs(self.lut_modtran_directory, exist_ok=True) + os.makedirs(self.lut_sixs_directory, exist_ok=True) + + # RTE Installation + if args.modtran_path: + self.modtran_path = args.modtran_path + else: + self.modtran_path = os.getenv("MODTRAN_DIR", env.modtran) -def build_main_config(paths, config_output_path, to_solar_azimuth_lut_grid: np.array, to_solar_zenith_lut_grid: np.array, - aerfrac_2_lut_grid: np.array, h2o_lut_grid: np.array = None, - elevation_lut_grid: np.array = None, altitude_lut_grid: np.array = None, to_sensor_azimuth_lut_grid: np.array = None, - to_sensor_zenith_lut_grid: np.array = None, n_cores: int = 1): + if args.sixs_path: + self.sixs_path = args.sixs_path + else: + self.sixs_path = os.getenv("SIXS_DIR", env.sixs) + + +def build_main_config( + paths, + h2o_lut_grid: np.array, + elevation_lut_grid: np.array, + altitude_lut_grid: np.array, + to_sensor_zenith_lut_grid: np.array, + to_solar_zenith_lut_grid: np.array, + to_sensor_azimuth_lut_grid: np.array, + to_solar_azimuth_lut_grid: np.array, + relative_azimuth_lut_grid: np.array, + aerfrac_2_lut_grid: np.array, + config_output_path, + configure_and_exit: bool = True, + ip_head: str = None, + redis_password: str = None, + n_cores: int = 1, +): """ Write an isofit dummy config file, so we can pass in for luts. Args: - paths: object with relevant path information attatched - config_output_path: path to write config to - to_solar_azimuth_lut_grid: the to-solar azimuth angle look up table grid to build - to_solar_zenith_lut_grid: the to-solar zenith angle look up table grid to build - aerfrac_2_lut_grid: the aerosol 2 look up table grid isofit should use for this solve - h2o_lut_grid: the water vapor look up table grid isofit should use for this solve - elevation_lut_grid: the ground elevation look up table grid isofit should use for this solve - altitude_lut_grid: the acquisition altitude look up table grid isofit should use for this solve - to_sensor_azimuth_lut_grid: the to-sensor azimuth angle look up table grid isofit should use for this solve - to_sensor_zenith_lut_grid: the to-sensor zenith angle look up table grid isofit should use for this solve - n_cores: the number of cores to use during processing + paths: Object containing references to all relevant file locations + lut_params: Configuration parameters for the lut grid + h2o_lut_grid: The water vapor look up table grid isofit should use for this solve + elevation_lut_grid: The ground elevation look up table grid isofit should use for this solve + altitude_lut_grid: The altitude elevation (km) look up table grid isofit should use for this solve + to_sensor_zenith_lut_grid: The to-sensor zenith angle look up table grid isofit should use for this solve + to_solar_zenith_lut_grid: The to-sun zenith angle look up table grid isofit should use for this solve + relative_azimuth_lut_grid: The relative to-sun azimuth angle look up table grid isofit should use for + config_output_path: Path to write config to + n_cores: The number of cores to use during processing """ - + # Initialize the RT config with the engines portion. radiative_transfer_config = { - "radiative_transfer_engines": { "modtran": { "engine_name": 'modtran', - "lut_path": paths.lut_modtran_directory, - "template_file": paths.modtran_template_path, - "wavelength_range": [350,2500], + "sim_path": paths.lut_modtran_directory, + "lut_path": os.path.join( + paths.lut_modtran_directory, + 'lut.nc' + ), + "rt_mode": 'rdn', + "irradiance_file": paths.irradiance_file, + "multipart_transmittance": True, + "template_file": paths.modtran_template_file, + "rte_configure_and_exit": configure_and_exit, + "engine_base_dir": paths.modtran_path, #lut_names - populated below #statevector_names - populated below }, "sixs": { "engine_name": '6s', - "lut_path": paths.lut_sixs_directory, - "wavelength_range": [350, 2500], + "rt_mode": 'rdn', + "sim_path": paths.lut_sixs_directory, + "lut_path": os.path.join( + paths.lut_sixs_directory, + 'lut.nc' + ), "irradiance_file": paths.irradiance_file, "earth_sun_distance_file": paths.earth_sun_distance_file, "month": 6, # irrelevant to readouts we care about @@ -226,7 +537,11 @@ def build_main_config(paths, config_output_path, to_solar_azimuth_lut_grid: np.a "viewaz": to_sensor_azimuth_lut_grid[0], "viewzen": 180 - to_sensor_zenith_lut_grid[0], "solaz": to_solar_azimuth_lut_grid[0], - "solzen": to_solar_zenith_lut_grid[0] + "solzen": to_solar_zenith_lut_grid[0], + "multipart_transmittance": True, + "template_file": paths.modtran_template_file, + "rte_configure_and_exit": configure_and_exit, + "engine_base_dir": paths.sixs_path, # lut_names - populated below # statevector_names - populated below } @@ -236,62 +551,237 @@ def build_main_config(paths, config_output_path, to_solar_azimuth_lut_grid: np.a "H2O_ABSCO": 0.0 } } + + # Question to figure out: Do I need the second key on a lot of these? + # H2O LUT Grid if h2o_lut_grid is not None and len(h2o_lut_grid) > 1: - radiative_transfer_config['lut_grid']['H2OSTR'] = [max(0.0, float(q)) for q in h2o_lut_grid] + radiative_transfer_config['lut_grid']['H2OSTR'] = [ + max(0.0, float(q)) for q in h2o_lut_grid + ] + # Elevation LUT Grid if elevation_lut_grid is not None and len(elevation_lut_grid) > 1: - radiative_transfer_config['lut_grid']['GNDALT'] = [max(0.0, float(q)) for q in elevation_lut_grid] - radiative_transfer_config['lut_grid']['elev'] = [float(q) for q in elevation_lut_grid] + radiative_transfer_config['lut_grid']['surface_elevation_km'] = [ + max(0.0, float(q)) for q in elevation_lut_grid + ] + # Altitude LUT Grid if altitude_lut_grid is not None and len(altitude_lut_grid) > 1: - radiative_transfer_config['lut_grid']['H1ALT'] = [max(0.0, float(q)) for q in altitude_lut_grid] - radiative_transfer_config['lut_grid']['alt'] = [float(q) for q in altitude_lut_grid] - - if to_sensor_azimuth_lut_grid is not None and len(to_sensor_azimuth_lut_grid) > 1: - radiative_transfer_config['lut_grid']['TRUEAZ'] = [float(q) for q in to_sensor_azimuth_lut_grid] - radiative_transfer_config['lut_grid']['viewaz'] = [float(q) for q in to_sensor_azimuth_lut_grid] - - if to_sensor_zenith_lut_grid is not None and len(to_sensor_zenith_lut_grid) > 1: - radiative_transfer_config['lut_grid']['OBSZEN'] = [float(q) for q in to_sensor_zenith_lut_grid] # modtran convension - radiative_transfer_config['lut_grid']['viewzen'] = [180 - float(q) for q in to_sensor_zenith_lut_grid] # sixs convension - - if to_solar_zenith_lut_grid is not None and len(to_solar_zenith_lut_grid) > 1: - radiative_transfer_config['lut_grid']['solzen'] = [float(q) for q in to_solar_zenith_lut_grid] # modtran convension - - if to_solar_azimuth_lut_grid is not None and len(to_solar_azimuth_lut_grid) > 1: - radiative_transfer_config['lut_grid']['solaz'] = [float(q) for q in to_solar_azimuth_lut_grid] # modtran convension - - # add aerosol elements from climatology + radiative_transfer_config['lut_grid']['observer_altitude_km'] = [ + max(0.0, float(q)) for q in altitude_lut_grid + ] + # radiative_transfer_config['lut_grid']['alt'] = [ + # float(q) for q in altitude_lut_grid + # ] + + # Observer Zenith LUT Grid + if ( + to_sensor_zenith_lut_grid is not None + and len(to_sensor_zenith_lut_grid) > 1 + ): + # Does Sixs? automatically convert value from Modtran convention + # Don't think so. Happens in template_construction.py + # Do I have to carry these as two keys? + # radiative_transfer_config['lut_grid']['observer_zenith'] = [ + # float(q) for q in to_sensor_zenith_lut_grid + # ] + # Sixs convension + # radiative_transfer_config['lut_grid']['viewzen'] = [ + radiative_transfer_config['lut_grid']['observer_zenith'] = [ + 180 - float(q) for q in to_sensor_zenith_lut_grid + ] + + # Solar Zenith LUT Grid + if ( + to_solar_zenith_lut_grid is not None + and len(to_solar_zenith_lut_grid) > 1 + ): + # Modtran convension + radiative_transfer_config['lut_grid']['solar_zenith'] = [ + float(q) for q in to_solar_zenith_lut_grid + ] + + # Azimuth LUT Grid + if ( + relative_azimuth_lut_grid is not None + and len(relative_azimuth_lut_grid) > 1 + ): + radiative_transfer_config["lut_grid"]["relative_azimuth"] = [ + float(q) for q in relative_azimuth_lut_grid + ] + + # Aerosol LUT Grid + # Should be able to use AERFRAC_2 for both if len(aerfrac_2_lut_grid) > 1: - radiative_transfer_config['lut_grid']['AERFRAC_2'] = [float(q) for q in aerfrac_2_lut_grid] - radiative_transfer_config['lut_grid']['AOT550'] = [float(q) for q in aerfrac_2_lut_grid] + radiative_transfer_config['lut_grid']['AERFRAC_2'] = [ + float(q) for q in aerfrac_2_lut_grid + ] + radiative_transfer_config['lut_grid']['AOT550'] = [ + float(q) for q in aerfrac_2_lut_grid + ] if paths.aerosol_model_path is not None: - radiative_transfer_config['radiative_transfer_engines']['modtran']['aerosol_model_file'] = paths.aerosol_model_path + radiative_transfer_config[ + 'radiative_transfer_engines' + ]['modtran']['aerosol_model_file'] = paths.aerosol_model_path + if paths.aerosol_tpl_path is not None: - radiative_transfer_config['radiative_transfer_engines']['modtran']["aerosol_template_file"] = paths.aerosol_tpl_path + radiative_transfer_config[ + 'radiative_transfer_engines' + ]['modtran']["aerosol_template_file"] = paths.aerosol_tpl_path # MODTRAN should know about our whole LUT grid and all of our statevectors, so copy them in - radiative_transfer_config['radiative_transfer_engines']['modtran']['lut_names'] = [x for x in ['H2OSTR','AERFRAC_2','GNDALT','H1ALT','TRUEAZ','OBSZEN', 'solzen','solaz'] if x in radiative_transfer_config['lut_grid'].keys()] - radiative_transfer_config['radiative_transfer_engines']['sixs']['lut_names'] = [x for x in ['H2OSTR','AOT550','elev','alt','viewaz','viewzen','solzen','solaz'] if x in radiative_transfer_config['lut_grid'].keys()] + # Populate the lut_names and the statevector_names + modtran_lut_names = { + x: None for x in [ + 'H2OSTR','surface_elevation_km','observer_altitude_km', + 'observer_zenith','solar_zenith', 'relative_azimuth', + 'AERFRAC_2' + ] if x in radiative_transfer_config['lut_grid'].keys() + } + radiative_transfer_config[ + 'radiative_transfer_engines' + ]['modtran']['lut_names'] = modtran_lut_names + radiative_transfer_config["radiative_transfer_engines"]["modtran"][ + "statevector_names"] = list(modtran_lut_names.keys()) + + sixs_lut_names = { + x: None for x in [ + 'H2OSTR','surface_elevation_km','observer_altitude_km', + # 'viewzen','solar_zenith', 'relative_azimuth', + 'observer_zenith','solar_zenith', 'relative_azimuth', + 'AOT550' + ] if x in radiative_transfer_config['lut_grid'].keys() + } + radiative_transfer_config[ + 'radiative_transfer_engines' + ]['sixs']['lut_names'] = sixs_lut_names + radiative_transfer_config["radiative_transfer_engines"]["sixs"][ + "statevector_names" + ] = list(sixs_lut_names.keys()) - # make isofit configuration - isofit_config_modtran = {'input': {}, - 'output': {}, - 'forward_model': { - "radiative_transfer": radiative_transfer_config, - "instrument": {"wavelength_file": paths.wavelenth_file} - }, - "implementation": { - "inversion": {"windows": [[1,2],[3,4]]}, - "n_cores": n_cores} - } + # Inversion windows - Not sure what to use here + # inversion_windows = [[1,2],[3,4]] + inversion_windows = [[350.0, 1360.0], [1410, 1800.0], [1970.0, 2500.0]] + # make isofit configuration + isofit_config_modtran = { + 'input': {}, + 'output': {}, + 'forward_model': { + "radiative_transfer": radiative_transfer_config, + "instrument": {"wavelength_file": paths.wavelength_file} + }, + "implementation": { + "inversion": {"windows": inversion_windows}, + "n_cores": n_cores, + "ip_head": ip_head, + "redis_password": redis_password + } + } isofit_config_modtran['implementation']["rte_configure_and_exit"] = True # write modtran_template with open(config_output_path, 'w') as fout: - fout.write(json.dumps(isofit_config_modtran, cls=SerialEncoder, indent=4, sort_keys=True)) + fout.write(json.dumps( + isofit_config_modtran, + cls=SerialEncoder, + indent=4, + sort_keys=True + )) + + +def sobel_lut(bounds, consts={}, n=3, log=False): + """ + Create a look using a Sobel sequence + + Args: + bounds (dict): keys are the lut element names. Values (list) are the bounds (idx 0 and 1) of the sequence. idx 2 is the number of samples to take + ns (list): Number of samples to take from sequence for each bounds element. The length of bounds + should math the length of n + consts (dict): keys are the lut element names. Values are the constant values to use + log (bool): Flag for log transform before sampling + Returns: + grid (dict): lut grid. + """ + + # Handle cos transform + cosd = lambda a : np.cos(np.deg2rad(a)) + cos_keys = [ + 'to_solar_zenith_bnds', + 'to_sensor_zenith_bnds', + 'to_solar_azimuth_bnds', + 'to_sensor_azimuth_bnds', + ] + for key, bnd in bounds.items(): + if key in cos_keys: + bnd = [cosd(bnd[0]), cosd(bnd[1])] + bounds[key] = bnd + + # Handle input log transform + log_keys = [ + 'altitude_km_bnds', + 'elevation_km_bnds', + 'h2o_bnds', + 'aerfrac_2_bnds', + ] + if log: + for key, bnd in bounds.items(): + if key in log_keys: + bnd = [np.log(bnd[0]), np.log(bnd[1])] + bounds[key] = bnd + + # Sort bounds + for key, bnd in bounds.items(): + bounds[key] = sorted(bnd) + + l_bound = [bnd[0] for bnd in bounds.values()] + r_bound = [bnd[1] for bnd in bounds.values()] + + # Make the sobel sequence + # Not 100% sure if we need to scamble the sequence gen + sobol = qmc.Sobol(d=len(bounds), scramble=False) + seq = sobol.random(n=n) + seq = qmc.scale( + seq, + l_bounds=l_bound, + u_bounds=r_bound + ) + + # Transform sample out of cos + arccosd = lambda a : np.rad2deg(np.arccos(a)) + cos_i = [i for i, key in enumerate(bounds.keys()) if key in cos_keys] + seq[:, cos_i] = arccosd(seq[:, cos_i]) + + # Transform out of log space + if log: + log_i = [i for i, key in enumerate(bounds.keys()) if key in log_keys] + seq[:, log_i] = np.exp(seq[:, log_i]) + + # Construct the grid-dict + grid = {} + for i, key in enumerate(bounds.keys()): + key_lut = '_'.join(key.split('_')[:-1] + ['lut']) + grid[key_lut] = [float(i) for i in seq[:, i]] + + # Add the constants onto the grid + for key, value in consts.items(): + key_lut = '_'.join(key.split('_')[:-1] + ['lut']) + grid[key_lut] = [float(value)] + + # Trim out non-physical components + # (e.g., elevation > observation altitude) + # The best way I can do this with uneven clolumns is to remove + # alt that are below max elev + grid['altitude_km_lut'] = np.array(grid['altitude_km_lut'])[ + grid['altitude_km_lut'] > np.max(grid['elevation_km_lut']) + ] + + # Do some cleaning a sorting + for key, value in grid.items(): + grid[key] = np.round(sorted(value), 4) + + return grid if __name__ == '__main__': diff --git a/convolve.py b/convolve.py new file mode 100644 index 0000000..fc708a2 --- /dev/null +++ b/convolve.py @@ -0,0 +1,214 @@ + + + + + + + +import numpy as np +import argparse +import h5py +import subprocess +import os +import ray +import multiprocessing as mp + + + +def spectral_response_function(response_range: np.array, mu: float, sigma: float): + """Calculate the spectral response function. + + Args: + response_range: signal range to calculate over + mu: mean signal value + sigma: signal variation + + Returns: + np.array: spectral response function + + """ + + u = (response_range - mu) / abs(sigma) + y = (1.0 / (np.sqrt(2.0 * np.pi) * abs(sigma))) * np.exp(-u * u / 2.0) + srf = y / y.sum() + return srf + + +def resample_spectrum( + x: np.array, wl: np.array, wl2: np.array, fwhm2: np.array, fill: bool = False, H: np.array = None +) -> np.array: + """Resample a spectrum to a new wavelength / FWHM. + Assumes Gaussian SRFs. + + Args: + x: radiance vector + wl: sample starting wavelengths + wl2: wavelengths to resample to + fwhm2: full-width-half-max at resample resolution + fill: boolean indicating whether to fill in extrapolated regions + + Returns: + np.array: interpolated radiance vector + + """ + if H is None: + H = np.array( + [ + spectral_response_function(wl, wi, fwhmi / 2.355) + for wi, fwhmi in zip(wl2, fwhm2) + ] + ) + H[np.isnan(H)] = 0 + + dims = len(x.shape) + if fill: + if dims > 1: + raise Exception("resample_spectrum(fill=True) only works with vectors") + + x = x.reshape(-1, 1) + xnew = np.dot(H, x).ravel() + good = np.isfinite(xnew) + for i, xi in enumerate(xnew): + if not good[i]: + nearest_good_ind = np.argmin(abs(wl2[good] - wl2[i])) + xnew[i] = xnew[nearest_good_ind] + return xnew + else: + # Replace NaNs with zeros + x[np.isnan(x)] = 0 + + # Matrix + if dims > 1: + return np.dot(H, x.T).T + + # Vector + else: + x = x.reshape(-1, 1) + return np.dot(H, x).ravel() + + +def main(): + parser = argparse.ArgumentParser(description="built luts for emulation.") + parser.add_argument('input_file', type=str) + parser.add_argument('output_file', type=str) + parser.add_argument('--output_sixs_wl_file', type=str, default=None) + parser.add_argument('--output_modtran_wl_file', type=str, default=None) + parser.add_argument('--irr_file', type=str, default=None) + args = parser.parse_args() + + output_sixs = None + output_sixs_wl = None + npzf = np.load(args.input_file) + + + # SIXS + input_sixs_wl = npzf['sixs_wavelengths'] + + + input_sixs = npzf['sixs_results'].copy() + input_solar_irr = npzf['sol_irr'] + points = npzf['points'] + point_names = npzf['point_names'].tolist() + + if args.output_sixs_wl_file is not None: + if args.irr_file is not None: + irr_wl, irr = np.loadtxt(args.irr_file, comments="#").T + irr = irr / 10 # convert to uW cm-2 sr-1 nm-1 + input_sixs_irr = np.array(resample_spectrum(irr, irr_wl, input_sixs_wl, np.ones(len(irr_wl))*(input_sixs_wl[1]-input_sixs_wl[0])),dtype=np.float32) + + n_bands = int(input_sixs.shape[1] / len(npzf['keys'])) + for key in range(len(npzf['keys'])): + input_sixs[:, n_bands * key : n_bands * (key + 1)] = input_sixs[:, n_bands * key : n_bands * (key + 1)] * input_sixs_irr * np.pi / np.cos(np.deg2rad(points[:, point_names.index('solar_zenith')]))[:,np.newaxis] + + output_sixs_wl = np.genfromtxt(args.output_sixs_wl_file)[:, 1] + output_sixs_fwhm = np.genfromtxt(args.output_sixs_wl_file)[:, 2] + if np.all(output_sixs_wl < 1000): + output_sixs_wl *= 1000 + output_sixs_fwhm *= 1000 + + output_sixs_irr = np.array(resample_spectrum(irr, irr_wl, output_sixs_wl, output_sixs_fwhm),dtype=np.float32) + H = np.array( [ spectral_response_function(input_sixs_wl, wi, fwhmi) for wi, fwhmi in zip(output_sixs_wl, output_sixs_fwhm) ] ) + H[np.isnan(H)] = 0 + output_sixs = np.zeros((input_sixs.shape[0], len(output_sixs_wl)*len(npzf['keys']))) + n_bands_o = len(output_sixs_wl) + for key in range(len(npzf['keys'])): + output_sixs[:, n_bands_o * key : n_bands_o * (key + 1)] = np.dot(input_sixs[:, n_bands * key : n_bands * (key + 1)], H.T) / output_sixs_irr / np.pi * np.cos(np.deg2rad(points[:, point_names.index('solar_zenith')]))[:,np.newaxis] + print(output_sixs.shape) + else: + output_sixs = input_sixs + output_sixs_wl = npzf['sixs_wavelengths'] + + + # MODTRAN + output_modtran = None + output_solar_irr = None + output_modtran_wl = None + input_modtran_wl = npzf['modtran_wavelengths'] + + + input_modtran = npzf['modtran_results'].copy() + input_solar_irr = npzf['sol_irr'] + points = npzf['points'] + point_names = npzf['point_names'].tolist() + + if args.output_modtran_wl_file is not None: + if args.irr_file is not None: + irr_wl, irr = np.loadtxt(args.irr_file, comments="#").T + irr = irr / 10 + input_modtran_irr = np.array(resample_spectrum(irr, irr_wl, input_modtran_wl, np.ones(len(irr_wl))*(input_modtran_wl[1] - input_modtran_wl[0])),dtype=np.float32) + + n_bands = int(input_modtran.shape[1] / len(npzf['keys'])) + #for key in range(len(npzf['keys'])): + # input_modtran[:, n_bands * key : n_bands * (key + 1)] = input_modtran[:, n_bands * key : n_bands * (key + 1)] * input_modtran_irr * np.pi / np.cos(np.deg2rad(points[:, point_names.index('solar_zenith')]))[:,np.newaxis] + + + output_modtran_wl = np.genfromtxt(args.output_modtran_wl_file)[:, 1] + output_modtran_fwhm = np.genfromtxt(args.output_modtran_wl_file)[:, 2] + if np.all(output_modtran_wl < 1000): + output_modtran_wl *= 1000 + output_modtran_fwhm *= 1000 + + output_modtran_irr = np.array(resample_spectrum(irr, irr_wl, output_modtran_wl, output_modtran_fwhm),dtype=np.float32) + H = np.array( [ spectral_response_function(input_modtran_wl, wi, fwhmi) for wi, fwhmi in zip(output_modtran_wl, output_modtran_fwhm) ] ) + H[np.isnan(H)] = 0 + output_modtran = np.zeros((input_modtran.shape[0], len(output_modtran_wl)*len(npzf['keys']))) + n_bands_o = len(output_modtran_wl) + input_modtran[np.isfinite(input_modtran) == False] = 0 + for key in range(len(npzf['keys'])): + print(n_bands_o * key , n_bands_o * (key + 1), n_bands * key , n_bands * (key + 1)) + output_modtran[:, n_bands_o * key : n_bands_o * (key + 1)] = np.dot(H, input_modtran[:, n_bands * key : n_bands * (key + 1)].T).T #/ output_modtran_irr / np.pi * np.cos(np.deg2rad(points[:, point_names.index('solar_zenith')]))[:,np.newaxis] + #import ipdb; ipdb.set_trace() + #test=resample_spectrum(input_modtran[0, n_bands * key : n_bands * (key + 1)], input_modtran_wl, output_modtran_wl, output_modtran_fwhm) + + print(output_modtran.shape) + + + else: + output_modtran = input_modtran + output_modtran_irr = npzf['sol_irr'] + output_modtran_wl = npzf['modtran_wavelengths'] + + + np.savez(args.output_file, modtran_results=output_modtran, + sixs_results=output_sixs, + points=npzf['points'], + sol_irr=output_modtran_irr, + sixs_wavelengths=output_sixs_wl, + modtran_wavelengths=output_modtran_wl, + point_names=npzf['point_names'], + keys=npzf['keys'] + ) + + + + + +if __name__ == "__main__": + main() + + + + + + + diff --git a/plot_munged.py b/plot_munged.py new file mode 100644 index 0000000..56c5966 --- /dev/null +++ b/plot_munged.py @@ -0,0 +1,166 @@ + + + + +import numpy as np +import matplotlib.pyplot as plt +import matplotlib.gridspec as gridspec +import os +import argparse +from matplotlib.lines import Line2D +from scipy import interpolate +import xarray as xr + +def d2_subset(data,ranges): + a = data.copy() + a = a[ranges[0],:] + a = a[:,ranges[1]] + return a + +def load( + file: str, + **kwargs, +) -> xr.Dataset: + + ds = xr.open_dataset(file, mode="r", lock=False, **kwargs) + dims = ds.drop_dims("wl").dims + + # Create the point dimension + ds = ds.stack(point=dims).transpose("point", "wl") + ds.load() + + return ds + +def main(): + + # Parse arguments + parser = argparse.ArgumentParser(description="built luts for emulation.") + parser.add_argument('input_netcdf', type=str) + parser.add_argument('--comparison_netcdf', type=str, default=None) + parser.add_argument('--input_netcdf_name', type=str, default=None) + parser.add_argument('--comparison_netcdf_name', type=str, default=None) + parser.add_argument('--variables', type=str, nargs='+', default=['rhoatm','sphalb','transm_down_dir','transm_down_dif', 'transm_up_dir','transm_up_dif']) + parser.add_argument('-fig_dir', type=str, default='figs') + + args = parser.parse_args() + + np.random.seed(13) + + rtm = load(args.input_netcdf) + rtm_points = np.vstack([x.data.tolist() for x in rtm.point]) + if args.comparison_netcdf is not None: + rtm_comp = load(args.comparison_netcdf) + rtm_comp_points = np.vstack([x.data.tolist() for x in rtm.point]) + if np.all(rtm_points == rtm_comp_points) is False: + print('Points do not match between input and comparison netcdf, terminating') + exit() + + for name in list(args.variables): + if name not in list(rtm.variables): + print(f'Could not find {name} in primary .nc, terminating') + exit() + if args.comparison_netcdf and name not in list(rtm_comp.variables): + print(f'Could not find {name} in secondary .nc, terminating') + exit() + + + + point_names = list(rtm.point.to_index().names) + + bad_points = np.zeros(rtm_points.shape[0],dtype=bool) + if 'surface_elevation_km' in point_names and 'observer_altitude_km' in point_names: + bad_points = rtm_points[:, point_names.index('surface_elevation_km')] >= rtm_points[:, point_names.index('observer_altitude_km')] -3 + bad_points[np.any(rtm_comp['transm_down_dif'].data[:,:] > 10,axis=1)] = True + bad_points[np.any(rtm['transm_down_dif'].data[:,:] > 10,axis=1)] = True + good_points = np.logical_not(bad_points) + bad_ind = np.any(rtm_comp['transm_down_dif'].data[:,:] > 10,axis=1) + + rtm['transm_up_dif'].data[rtm['transm_up_dif'].data[:,:] > 1] = 0 + rtm['transm_up_dif'].data[rtm['transm_up_dif'].data[:,:] < 0] = 0 + rtm_comp['transm_up_dif'].data[rtm_comp['transm_up_dif'].data[:,:] > 1] = 0 + rtm_comp['transm_up_dif'].data[rtm_comp['transm_up_dif'].data[:,:] < 0] = 0 + + rtm['sphalb'].data[:,rtm.wl > 1200][rtm['sphalb'].data[:,rtm.wl > 1200] > 0.1] = 0 + rtm['sphalb'].data[rtm['sphalb'].data[:,:] < 0] = 0 + + print(np.sum(rtm_comp['sphalb'].data[:,rtm_comp.wl > 1200] > 0.1)) + subs = rtm_comp['sphalb'].data[:,rtm_comp.wl > 1200].copy() + subs[subs > 0.1] = 0 + rtm_comp['sphalb'].data[:,rtm_comp.wl > 1200] = subs + print(np.sum(rtm_comp['sphalb'].data[:,rtm_comp.wl > 1200] > 0.1)) + + rtm_comp['sphalb'].data[rtm_comp['sphalb'].data[:,:] < 0] = 0 + + lims_main = [[0, 1], [0, 0.25], [0, 0.35]] + lims_diff = [[0, 0.25], [0, 0.1], [0, 0.1]] + + cmap = plt.get_cmap('coolwarm') + for dim in range(len(point_names)): + + fig = plt.figure(figsize=(30, 10)) + gs = gridspec.GridSpec(ncols=len(args.variables), nrows=2,wspace=0.2, hspace=0.2) + plt.suptitle(point_names[dim]) + + #slice = np.take(points,np.arange(0,points.shape[0]),axis=dim) + slice = rtm_points[:,dim] + un_vals = np.unique(slice) + + for key_ind, key in enumerate(args.variables): + ax = fig.add_subplot(gs[0,key_ind]) + plt.title(f'{args.input_netcdf_name}: {args.variables[key_ind]}') + + leg_lines = [] + leg_names = [] + for _val, val in enumerate(un_vals): + rtm_slice = np.nanmean(rtm[key].data[np.logical_and(slice == val, good_points), :], axis=0) + plt.plot(rtm.wl, rtm_slice, c=cmap(float(_val)/len(un_vals)), linewidth=1) + + leg_lines.append(Line2D([0], [0], color=cmap(float(_val)/len(un_vals)), lw=2)) + leg_names.append(str(round(val,2))) + plt.grid() + + #plt.ylim(lims_main[key_ind]) + plt.xlabel('Wavelength [nm]') + + if key_ind == 1: + + if point_names[dim] == 'H2OSTR': + pn = 'Water Vapor\n[g cm$^{-1}$]' + elif point_names[dim] == 'AOT550' or point_names[dim] == 'AERFRAC_2': + pn = 'Aerosol Optical\nDepth' + elif point_names[dim] == 'observer_altitude_km': + pn = 'Observer\nAltitude [km]' + elif point_names[dim] == 'surface_elevation_km': + pn = 'Surface\nElevation [km]' + elif point_names[dim] == 'solar_zenith': + pn = 'Solar\nZenith Angle [deg]' + else: + pn = point_names[dim] + + plt.legend(leg_lines, leg_names, title=pn) + + #pointstr = '{}: {} - {}'.format(point_names[dim], un_vals[0],un_vals[-1]) + #plt.text(50, 0.9 * (lims_main[key_ind][1] - lims_main[key_ind][0]) + lims_main[key_ind][0], pointstr, + # verticalalignment='top') + #elif key_ind == 2 and args.comparison_netcdf: + # + # leg_lines = [Line2D([0], [0], color='black', lw=2), + # Line2D([0], [0], color='black', lw=2, ls='--') + # ] + # plt.legend(leg_lines, [args.input_netcdf_name, args.comparison_netcdf_name], title='RTM') + + if args.comparison_netcdf: + for key_ind, key in enumerate(args.variables): + ax = fig.add_subplot(gs[1,key_ind]) + plt.title(f'{args.comparison_netcdf_name}: {key}') + for _val, val in enumerate(un_vals): + rtm_comp_slice = np.nanmean(rtm_comp[key].data[np.logical_and(slice == val, good_points), :], axis=0) + plt.plot(rtm_comp.wl, rtm_comp_slice, c=cmap(float(_val)/len(un_vals)), linewidth=1) + plt.grid() + + plt.savefig(f'{args.fig_dir}/dim_{point_names[dim]}.png', dpi=200, bbox_inches='tight') + plt.clf() + + +if __name__ == '__main__': + main() diff --git a/train_combined_emulator.py b/train_combined_emulator.py index 037fa58..ed4971d 100644 --- a/train_combined_emulator.py +++ b/train_combined_emulator.py @@ -6,8 +6,8 @@ import matplotlib.pyplot as plt import matplotlib.gridspec as gridspec import os -from isofit.radiative_transfer.modtran import ModtranRT -from isofit.radiative_transfer.six_s import SixSRT +#from isofit.radiative_transfer.modtran import ModtranRT +#from isofit.radiative_transfer.six_s import SixSRT from isofit.configs import configs import argparse from sklearn import linear_model @@ -22,25 +22,79 @@ from scipy import interpolate import pickle +#os.environ["CUDA_VISIBLE_DEVICES"] = "-1" + + +def rho_to_rdn(rho,solar_irr,coszen): + return rho / np.pi * solar_irr[np.newaxis,:] * coszen + def beckman_rdn(simulated_modtran, wl, n_bands=424): - refl_file = "../isofit/examples/20171108_Pasadena/insitu/BeckmanLawn.txt" - solar_irr = np.load('../isofit/examples/20171108_Pasadena/solar_irr.npy')[:-1] - coszen = 0.6155647578988601 - rfl = np.genfromtxt(refl_file) - rfl = resample_spectrum(rfl[:,1],rfl[:,0],wl,rfl[:,2]*1000) - rho = simulated_modtran[:,n_bands:2*n_bands] - rdn_atm = rho / np.pi*(solar_irr * coszen) + refl_file = "/store/brodrick/repos/isofit-tutorials/20171108_Pasadena/insitu/BeckmanLawn.txt" + solar_irr = np.genfromtxt('/store/brodrick/repos/isofit-data/kurudz_0.1nm.dat')[:-1] - transm = simulated_modtran[:,:n_bands] - rdn_down = (solar_irr * coszen) / np.pi * transm - sphalb = simulated_modtran[:,n_bands*2:n_bands*3] + coszen = 0.6155647578988601 + rfl = np.genfromtxt(refl_file) + fwhm = rfl[:,2]*1000 + print(fwhm[0]) + rfl = resample_spectrum(rfl[:,1],rfl[:,0],wl,fwhm) + + irr = np.loadtxt('/store/brodrick/repos/isofit-data/kurudz_0.1nm.dat', comments="#") + iwl, irr = irr.T + irr = irr / 10.0 # convert, uW/nm/cm2 + #irr = irr / self.irr_factor**2 # consider solar distance + solar_irr = resample_spectrum(irr, iwl, wl, fwhm) + + n_bands = len(wl) + + #keys=['rhoatm','sphalb','transm_down_dir','transm_down_dif', 'transm_up_dir','transm_up_dif'] + rho = simulated_modtran[:,:n_bands] + sphalb = simulated_modtran[:,n_bands*1:n_bands*2] + transm_down_dir = simulated_modtran[:,n_bands*2:n_bands*3] + transm_down_dif = simulated_modtran[:,n_bands*3:n_bands*4] + transm_up_dir = simulated_modtran[:,n_bands*4:n_bands*5] + transm_up_dif = simulated_modtran[:,n_bands*5:n_bands*6] + + + # bi-directional (downward direct * upward direct) + # hemispherical-directional (downward diffuse * upward direct) + # directional-hemispherical (downward direct * upward diffuse) + # bi-hemispherical (downward diffuse * upward diffuse) + ## at-sensor radiance model, including topography, adjacency effects, and glint + #ret = ( + # L_atm + # + ( + # L_bi_direct * rfl_dir # bi-directional radiance + # + L_hemi_direct * rfl_dif # hemispherical-directional radiance + # + L_direct_hemi * bg_dir # directional-hemispherical radiance + # + L_bi_hemi * bg_dif # bi-hemispherical radiance + # ) + # / (1.0 - s_alb * bg_dif) + # + L_up + #) + + + L_bi_direct = rho_to_rdn(transm_down_dir * transm_up_dir, solar_irr, coszen) + L_hemi_direct = rho_to_rdn(transm_down_dif * transm_up_dir, solar_irr, coszen) + L_direct_hemi = rho_to_rdn(transm_down_dir * transm_up_dif, solar_irr, coszen) + L_bi_hemi = rho_to_rdn(transm_down_dif * transm_up_dif, solar_irr, coszen) + L_atm = rho_to_rdn(rho, solar_irr, coszen) + + rdn = L_atm +\ + (+\ + L_bi_direct * rfl +\ + L_hemi_direct * rfl +\ + L_direct_hemi * rfl +\ + L_bi_hemi * rfl\ + ) / (1.0 - sphalb * rfl) + + + rdn_atm = rho_to_rdn(rho, solar_irr, coszen) - rdn = rdn_atm + rdn_down * rfl / (1.0 - sphalb * rfl) return rdn, rdn_atm @@ -60,12 +114,42 @@ def nn_model(in_data_shape, out_data_shape, num_layers=5): output_layer = keras.layers.Dense(units=out_data_shape[-1], activation='linear')(output_layer) model = keras.models.Model(inputs=[inlayer], outputs=[output_layer]) optimizer=keras.optimizers.Adam(learning_rate=0.0001) + #optimizer=keras.optimizers.Adam() model.compile(loss='mse', optimizer=optimizer) return model +def nn_model_ind(in_data_shape, out_data_shape, num_keys, num_layers=1): + + layer_depths = np.linspace(int(in_data_shape[-1]/num_keys),int(out_data_shape[-1]/num_keys),num=num_layers+1,dtype=int) + + inlayer = keras.layers.Input(shape=(in_data_shape[-1],)) + + instack = [] + outstack = [] + input_width = int(in_data_shape[-1]/num_keys) + for nk in range(num_keys): + #output_layer = keras.layers.Cropping1D(cropping=(0,(num_keys-1)*(in_data_shape[-1]/num_keys)))(inlayer) + output_layer = inlayer[:, nk*input_width:(nk+1)*input_width] + + for _l in range(num_layers): + output_layer = keras.layers.Dense(units=layer_depths[_l])(output_layer) + output_layer = keras.layers.LeakyReLU(alpha=0.4)(output_layer) + + output_layer = keras.layers.Dense(units=int(out_data_shape[-1]/num_keys), activation='linear')(output_layer) + outstack.append(output_layer) + + output_merge = keras.layers.Concatenate()(outstack) + + model = keras.models.Model(inputs=[inlayer], outputs=[output_merge]) + optimizer=keras.optimizers.Adam(learning_rate=0.001) + #optimizer=keras.optimizers.Adam() + model.compile(loss='mse', optimizer=optimizer) + + return model + class SplitModel(): def __init__(self,in_data_shape, out_data_shape, n_splits): @@ -117,11 +201,22 @@ def main(): np.random.seed(13) + #es = keras.callbacks.EarlyStopping(monitor='val_loss', mode='min', verbose=1, patience=20, restore_best_weights=True) + #train_sixs = np.ones((13905,5166)) + #train_modtran = np.ones((13905, 129006)) + #keys = np.ones(6) + #model = nn_model_ind(train_sixs.shape, train_modtran.shape, len(keys)) + #print(model.summary()) + #model.fit(train_sixs, train_modtran, batch_size=10, epochs=400, + # validation_data=(train_sixs, train_modtran),callbacks=[es]) + #exit() - npzf = np.load(args.munged_file) - modtran_results = npzf['modtran_results'] + + npzf = np.load(args.munged_file, allow_pickle=True) + sixs_results = npzf['sixs_results'] + modtran_results = npzf['modtran_results'] points = npzf['points'] keys = npzf['keys'] point_names = npzf['point_names'] @@ -134,6 +229,37 @@ def main(): n_bands_modtran = int(modtran_results.shape[-1]/len(keys)) n_bands_sixs = int(sixs_results.shape[-1]/len(keys)) + modtran_results[np.isnan(modtran_results)] = 0 + sixs_results[np.isnan(sixs_results)] = 0 + modtran_results[np.isfinite(modtran_results) == False] = 0 + sixs_results[np.isfinite(sixs_results) == False] = 0 + + max_wl = np.min([np.max(simulator_wavelengths),np.max(emulator_wavelengths)]) + min_wl = np.max([np.min(simulator_wavelengths),np.min(emulator_wavelengths)]) + emulator_idx = np.where(np.logical_and(emulator_wavelengths >= min_wl, emulator_wavelengths <= max_wl))[0] + + modtran_results_clipped = np.zeros((modtran_results.shape[0], len(emulator_idx)*len(keys))) + for n in range(len(keys)): + modtran_results_clipped[:,n*len(emulator_idx):(n+1)*len(emulator_idx)] = modtran_results[:, emulator_idx + n_bands_modtran*n] + + print(modtran_results.shape) + modtran_results = modtran_results_clipped + del modtran_results_clipped + n_bands_modtran = int(modtran_results.shape[-1]/len(keys)) + emulator_wavelengths = emulator_wavelengths[emulator_idx] + + print(modtran_results.shape) + + #fig = plt.figure(figsize=(20,5)) + #gs = gridspec.GridSpec(ncols=1, nrows=2, wspace=0.3, hspace=0.4) + #ax = fig.add_subplot(gs[0, 0]) + #plt.plot(np.mean(modtran_results,axis=0)) + #ax = fig.add_subplot(gs[1, 0]) + #plt.plot(np.mean(sixs_results,axis=0)) + #plt.savefig('figs/h2o_comp/inputs.png',dpi=200,bbox_inches='tight') + #exit() + + if args.holdout_dim == -1: perm = np.random.permutation(points.shape[0]) @@ -174,6 +300,8 @@ def main(): train_modtran = modtran_results-sixs_results_match_modtran + train_modtran[np.isfinite(train_modtran) == False] = 0 + train_sixs[np.isfinite(train_sixs) == False] = 0 print(train_modtran.shape) base_save_name = os.path.join(args.save_dir,'emulator') @@ -187,10 +315,14 @@ def main(): monitor='val_loss' es = keras.callbacks.EarlyStopping(monitor=monitor, mode='min', verbose=1, patience=20, restore_best_weights=True) - model = nn_model(train_sixs.shape, modtran_results.shape) + #model = nn_model_ind(train_sixs.shape, train_modtran.shape, len(keys)) + model = nn_model(train_sixs.shape, train_modtran.shape) + print(model.summary()) + simple_response_scaler = np.ones(train_modtran.shape[1])*100 train_modtran *= simple_response_scaler + #import ipdb; ipdb.set_trace() model.fit(train_sixs[train,:], train_modtran[train,:], batch_size=1000, epochs=400, validation_data=(train_sixs[test,:], train_modtran[test,:]),callbacks=[es]) train_modtran /= simple_response_scaler @@ -200,13 +332,11 @@ def main(): pred = full_pred[test, :] - model.save(base_save_name) + model.save(base_save_name + '.h5') np.savez(base_save_name + '_aux.npz', lut_names=point_names, rt_quantities=keys, - feature_scaler_mean=feature_scaler.mean_,response_scaler_mean=response_scaler.mean_, - feature_scaler_var=feature_scaler.var_,response_scaler_var=response_scaler.var_, - feature_scaler_scale=feature_scaler.scale_,response_scaler_scale=response_scaler.scale_, - solar_irr=solar_irr, emulator_wavelengths=emulator_wavelengths, + solar_irr=solar_irr, + emulator_wavelengths=emulator_wavelengths, response_scaler=simple_response_scaler, simulator_wavelengths=simulator_wavelengths) @@ -221,16 +351,27 @@ def main(): - fig = plt.figure(figsize=(10, 10 * 2 / (0.5 + 2))) + fig = plt.figure(figsize=(20, 20 * 2 / (0.5 + 2))) gs = gridspec.GridSpec(ncols=len(keys), nrows=2, wspace=0.3, hspace=0.4) - ref_rdn = np.genfromtxt('../isofit/examples/20171108_Pasadena/remote/ang20171108t184227_rdn_v2p11_BeckmanLawn_424.txt')[:,1] + ref_rdn = np.genfromtxt('/store/brodrick/repos/isofit-tutorials/20171108_Pasadena/remote/ang20171108t184227_rdn_v2p11_BeckmanLawn.txt') rdn_modtran, rdn_modtran_atm = beckman_rdn(modtran_results,emulator_wavelengths) + ref_rdn = resample_spectrum(ref_rdn[:,1],ref_rdn[:,0],emulator_wavelengths, np.ones(len(emulator_wavelengths))*7.5) + cf = 100 - varset = np.where(np.logical_and.reduce((points[:,3] == 2.55, points[:,-1] == 0, points[:,-2] == 180, points[:,1] == 0, points[:,2] == 4, points[:,0] == 0.01)))[0] + + varset = np.where(np.logical_and.reduce((points[:,0] == np.median(points[:,0]), + points[:,1] == np.median(points[:,1]), + points[:,2] == np.median(points[:,2]), + points[:,3] == np.median(points[:,3]), + points[:,4] == np.median(points[:,4]), + points[:,5] == np.median(points[:,5]), + points[:,6] == np.median(points[:,6]), + )))[0] + best_modtran = np.argmin(np.sum(np.power(rdn_modtran[:,:cf] - ref_rdn[:cf],2),axis=1)) best_emu = np.argmin(np.sum(np.power(rdn[:,:cf] - ref_rdn[:cf],2),axis=1)) plt.plot(emulator_wavelengths, ref_rdn, c='black', linewidth=0.8) @@ -252,7 +393,7 @@ def main(): plt.text(1000, 10 , pointstr_modtran, verticalalignment='top') plt.text(2000, 10 , pointstr_emu, verticalalignment='top') - plt.savefig('rdn_plots/best_matches.png',dpi=200,bbox_inches='tight') + plt.savefig(os.path.join(args.fig_dir, 'rdn_plots_best_matches.png'),dpi=200,bbox_inches='tight') @@ -278,7 +419,7 @@ def main(): pointstr = '{}: {} - {}'.format(point_names[dim], un_vals[0],un_vals[-1]) plt.text(2000, 8 , pointstr, verticalalignment='top') - plt.savefig('{}/dim_{}.png'.format('rdn_plots', dim), dpi=200, bbox_inches='tight') + plt.savefig('{}/dim_{}.png'.format(args.fig_dir, dim), dpi=200, bbox_inches='tight') plt.clf() @@ -302,15 +443,18 @@ def main(): pointstr = '{}: {} - {}'.format(point_names[dim], un_vals[0],un_vals[-1]) plt.text(2000, 8 , pointstr, verticalalignment='top') - plt.savefig('{}/modtran_dim_{}.png'.format('rdn_plots', dim), dpi=200, bbox_inches='tight') + plt.savefig('{}/modtran_dim_{}.png'.format(args.fig_dir, dim), dpi=200, bbox_inches='tight') plt.clf() - fig = plt.figure(figsize=(10, 10 * 3 / (0.5 + 3))) - gs = gridspec.GridSpec(ncols=len(keys), nrows=3, wspace=0.3, hspace=0.4) - print_keys = ['Total\nTransmittance', 'Atmospheric Path\nReflectance', 'Spherical Albedo'] + #print_keys = ['Total\nTransmittance', 'Atmospheric Path\nReflectance', 'Spherical Albedo'] + print_keys=['rhoatm','sphalb','transm_down_dir','transm_down_dif', 'transm_up_dir','transm_up_dif'] n_bands = int(modtran_results.shape[-1]/len(keys)) + + fig = plt.figure(figsize=(20, 20 * 2 / (0.5 + len(print_keys)))) + gs = gridspec.GridSpec(ncols=len(keys), nrows=3, wspace=0.3, hspace=0.4) + for key_ind, key in enumerate(keys): ax = fig.add_subplot(gs[0, key_ind]) @@ -385,7 +529,8 @@ def main(): plt.savefig('{}/mean_test_set.png'.format(args.fig_dir), bbox_inches='tight') plt.clf() - fig = plt.figure(figsize=(10, 10 * 2 / (0.5 + 2))) + #fig = plt.figure(figsize=(10, 10 * 2 / (0.5 + 2))) + fig = plt.figure(figsize=(20, 20 * 2 / (0.5 + len(print_keys)))) gs = gridspec.GridSpec(ncols=len(keys), nrows=2, wspace=0.3, hspace=0.4) error_mean = np.zeros(len(test)) @@ -397,8 +542,8 @@ def main(): order = np.argsort(error_mean) #order = np.argsort(points[test,0]) - lims_main = [[0, 1], [0, 0.25], [0, 0.35],[0,1]] - lims_diff = [[0, 0.25], [0, 0.1], [0, 0.1],[0,1]] + lims_main = [[0, 1], [0, 1], [0, 1], [0, 1], [0, 1], [0, 1]] + lims_diff = [[0, 0.25], [0, 0.25], [0, 0.25], [0, 0.25], [0, 0.25], [0, 0.25]] #for row in range(np.sum(test)): for row_ind, row in enumerate(order[np.linspace(0,len(order)-1,20,dtype=int)].tolist()):