From 5ca2363c4fc6595129e88e5b4c4c65060f4fe3be Mon Sep 17 00:00:00 2001 From: Evan Greenberg Date: Sat, 7 Dec 2024 11:09:30 -0800 Subject: [PATCH 01/14] Configure_and_exit working. Sobol in progress --- build_luts.py | 617 ++++++++++++++++++++++++++++++++++++++------------ 1 file changed, 469 insertions(+), 148 deletions(-) diff --git a/build_luts.py b/build_luts.py index 9afd125..30c6797 100644 --- a/build_luts.py +++ b/build_luts.py @@ -18,35 +18,65 @@ 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 + +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( + '-sixs_path', + default=None, + help='SIXS Installation DIR' + ) + parser.add_argument( + '-modtran_path', + default=None, + help='MODTRAN Installation DIR' + ) args = parser.parse_args() - args.train = args.train == 1 + args.training = args.train == 1 args.cleanup = args.cleanup == 1 + dayofyear = 200 + paths = Paths( + args, + args.training + ) + if args.train: to_solar_zenith_lut = [0, 12.5, 25, 37.5, 50] to_solar_azimuth_lut = [180] @@ -54,10 +84,22 @@ def main(): 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 = ( + 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 = ( + np.round(np.linspace(0.01, 1, num=5), 3).tolist() + [0.5] + ) aerfrac_2_lut_grid.sort() + relative_azimuth_lut = np.abs( + np.array(to_solar_azimuth_lut) + - np.array(to_sensor_azimuth_lut) + ) + relative_azimuth_lut = np.minimum( + relative_azimuth_lut, + 360 - relative_azimuth_lut + ) else: # HOLDOUT SET @@ -67,19 +109,26 @@ def main(): 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)]) + 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) + relative_azimuth_lut = np.abs( + np.array(to_solar_azimuth_lut) + - np.array(to_sensor_azimuth_lut) + ) + relative_azimuth_lut = np.minimum( + relative_azimuth_lut, + 360 - relative_azimuth_lut + ) + + n_lut_build = np.prod([ + len(to_solar_zenith_lut), + len(to_sensor_zenith_lut), + len(relative_azimuth_lut), + len(altitude_km_lut), + len(elevation_km_lut), + len(h2o_lut_grid), + len(aerfrac_2_lut_grid) + ]) print('Num LUTs to build: {}'.format(n_lut_build)) print('Expected MODTRAN runtime: {} hrs'.format(n_lut_build*1.5)) @@ -92,70 +141,135 @@ def main(): 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) + 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=altitude_km_lut[0], + dayofyear=dayofyear, + to_sensor_azimuth=to_sensor_azimuth_lut[0], + to_sensor_zenith=to_sensor_zenith_lut[0], + to_sun_zenith=to_solar_azimuth_lut[0], + relative_azimuth=relative_azimuth_lut[0], + gmtime=0, + elevation_km=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) - - - isofit_modtran = ModtranRT(config.forward_model.radiative_transfer.radiative_transfer_engines[0], - config) + 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, + h2o_lut_grid, + elevation_km_lut, + altitude_km_lut, + to_sensor_zenith_lut, + to_solar_zenith_lut, + to_sensor_azimuth_lut, + to_solar_azimuth_lut, + relative_azimuth_lut, + aerfrac_2_lut_grid, + 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, + ) - isofit_sixs = SixSRT(config.forward_model.radiative_transfer.radiative_transfer_engines[1], - config) + 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 = 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['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['lut_grid'] = { + key: params['lut_grid'][key] for key in + sixs_engine_config.lut_names.keys() + } + # params['modtran_emulation'] = True + 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 +277,145 @@ 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.modtran_template_path = modtran_tpl + self.support_dir = os.path.join(working_dir, 'support') + self.template_dir = os.path.join(working_dir, 'templates') + # Make dirs + os.makedirs(self.support_dir, exist_ok=True) + os.makedirs(self.template_dir, exist_ok=True) -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): + 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' + ) + 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, + 'prism_optimized_irr.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) + + 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' + ), + "multipart_transmittance": False, + "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], + "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 +425,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": False, + "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 +439,180 @@ 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', + 'AERFRAC_2' + ] 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 sobol_lut(): + 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() + + # Create a Sobol sequence generator + sobol = qmc.Sobol(d=2, scramble=False) # d is the number of dimensions + + # Generate 100 Sobol points + sample = sobol.random(n=100) + + # Scale to the desired range (e.g., [0, 1]) + sample = qmc.scale(sample, l_bounds=[0, 0], u_bounds=[1, 1]) + + + relative_azimuth_lut = np.abs( + np.array(to_solar_azimuth_lut) + - np.array(to_sensor_azimuth_lut) + ) + relative_azimuth_lut = np.minimum( + relative_azimuth_lut, + 360 - relative_azimuth_lut + ) if __name__ == '__main__': From 3ea1b9e3edcbf78ba3e98d5311e7a390c9c6ba23 Mon Sep 17 00:00:00 2001 From: Evan Greenberg Date: Sat, 7 Dec 2024 14:37:11 -0800 Subject: [PATCH 02/14] Sobel almost working. Worker issue. --- build_luts.py | 289 ++++++++++++++++++++++++++++++++++---------------- 1 file changed, 195 insertions(+), 94 deletions(-) diff --git a/build_luts.py b/build_luts.py index 30c6797..7b899b2 100644 --- a/build_luts.py +++ b/build_luts.py @@ -54,6 +54,11 @@ def main(): default=False, action=argparse.BooleanOptionalAction ) + parser.add_argument( + '--sobel', + default=False, + action=argparse.BooleanOptionalAction + ) parser.add_argument( '-sixs_path', default=None, @@ -69,7 +74,6 @@ def main(): args.training = args.train == 1 args.cleanup = args.cleanup == 1 - dayofyear = 200 paths = Paths( @@ -77,57 +81,91 @@ def main(): args.training ) - 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() - relative_azimuth_lut = np.abs( - np.array(to_solar_azimuth_lut) - - np.array(to_sensor_azimuth_lut) - ) - relative_azimuth_lut = np.minimum( - relative_azimuth_lut, - 360 - relative_azimuth_lut - ) - + if args.sobel: + if args.train: + # Get bounds from regular grid + bounds = { + 'to_solar_zenith_bnds': [0, 50, 2], + 'to_sensor_zenith_bnds': [140, 180, 2], + 'altitude_km_bnds': [2, 25, 2], + 'elevation_km_bnds': [0.01, 4.5, 2], + 'h2o_bnds': [0.1, 5, 3], + 'aerfrac_2_bnds': [0.01, 1, 3], + } + consts = { + 'to_solar_azimuth_bnds': 180, + 'to_sensor_azimuth_bnds': 180, + } + grid = sobel_lut(bounds, consts) + else: + # Get bounds from regular grid + bounds = { + 'to_solar_zenith_bnds': [0, 50, 3], + 'to_sensor_zenith_bnds': [140, 180, 3], + 'altitude_km_bnds': [2, 25, 3], + 'elevation_km_bnds': [0.01, 4.5, 3], + 'h2o_bnds': [0.1, 5, 5], + 'aerfrac_2_bnds': [0.01, 1, 5], + } + consts = { + 'to_solar_azimuth_bnds': 180, + 'to_sensor_azimuth_bnds': 180, + } + rid = sobel_lut(bounds, consts) 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) - relative_azimuth_lut = np.abs( - np.array(to_solar_azimuth_lut) - - np.array(to_sensor_azimuth_lut) - ) - relative_azimuth_lut = np.minimum( - relative_azimuth_lut, - 360 - relative_azimuth_lut - ) + if args.train: + grid = { + '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': list(np.sort( + np.round(np.linspace(0.1, 5, num=5), 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': [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': 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(to_solar_zenith_lut), - len(to_sensor_zenith_lut), - len(relative_azimuth_lut), - len(altitude_km_lut), - len(elevation_km_lut), - len(h2o_lut_grid), - len(aerfrac_2_lut_grid) + 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)) @@ -149,14 +187,14 @@ def main(): write_modtran_template( atmosphere_type='ATM_MIDLAT_SUMMER', fid=os.path.splitext(paths.modtran_template_file)[0], - altitude_km=altitude_km_lut[0], + altitude_km=grid['altitude_km_lut'][0], dayofyear=dayofyear, - to_sensor_azimuth=to_sensor_azimuth_lut[0], - to_sensor_zenith=to_sensor_zenith_lut[0], - to_sun_zenith=to_solar_azimuth_lut[0], - relative_azimuth=relative_azimuth_lut[0], + to_sensor_azimuth=grid['to_sensor_azimuth_lut'][0], + to_sensor_zenith=grid['to_sensor_zenith_lut'][0], + to_sun_zenith=grid['to_solar_azimuth_lut'][0], + relative_azimuth=grid['relative_azimuth_lut'][0], gmtime=0, - elevation_km=elevation_km_lut[0], + elevation_km=grid['elevation_km_lut'][0], output_file=paths.modtran_template_file, ihaze_type="AER_NONE", ) @@ -179,15 +217,15 @@ def main(): configure_and_exit = args.configure_and_exit build_main_config( paths, - h2o_lut_grid, - elevation_km_lut, - altitude_km_lut, - to_sensor_zenith_lut, - to_solar_zenith_lut, - to_sensor_azimuth_lut, - to_solar_azimuth_lut, - relative_azimuth_lut, - aerfrac_2_lut_grid, + 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, @@ -222,8 +260,7 @@ def main(): # Inits of engines based on radiative_transfer.py rt_config = config.forward_model.radiative_transfer instrument_config = config.forward_model.instrument - - lut_grid = rt_config.lut_grid + lut_grid_config = rt_config.lut_grid _keys = [ "interpolator_style", @@ -579,40 +616,104 @@ def build_main_config( )) -def sobol_lut(): - 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() +def sobel_lut(bounds, consts={}, log=False): + """ + Create a look using a Sobel sequence - # Create a Sobol sequence generator - sobol = qmc.Sobol(d=2, scramble=False) # d is the number of dimensions + 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. + """ - # Generate 100 Sobol points - sample = sobol.random(n=100) + # 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]), bnd[2]] + 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]), bnd[2]] + bounds[key] = bnd + + # Sort bounds + for key, bnd in bounds.items(): + bounds[key] = sorted(bnd[:-1]) + [bnd[2]] + + # Make the sobel sequence + grid = {} + for key, value in bounds.items(): + l_bound, r_bound = value[:-1] + + # Not 100% sure if we need to scamble the sequence gen + sobol = qmc.Sobol(d=1, scramble=False) + seq = sobol.random(n=value[-1]) + seq = qmc.scale( + seq, + l_bounds=l_bound, + u_bounds=r_bound + ) - # Scale to the desired range (e.g., [0, 1]) - sample = qmc.scale(sample, l_bounds=[0, 0], u_bounds=[1, 1]) + # Transform sample out of cos + if key in cos_keys: + seq = np.rad2deg(np.arccos(seq)) + # Transform out of log space + if log: + if key in log_keys: + seq = np.exp(seq) - relative_azimuth_lut = np.abs( - np.array(to_solar_azimuth_lut) - - np.array(to_sensor_azimuth_lut) - ) - relative_azimuth_lut = np.minimum( - relative_azimuth_lut, - 360 - relative_azimuth_lut + key_lut = '_'.join(key.split('_')[:-1] + ['lut']) + grid[key_lut] = [float(i) for i in seq] + + # Add the constants onto the grid + for key, value in consts.items(): + key_lut = '_'.join(key.split('_')[:-1] + ['lut']) + grid[key_lut] = [float(value)] + + # 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 + )] + + # 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__': From 81e27f1d7a152c2d5d6d51377288530211e9577f Mon Sep 17 00:00:00 2001 From: brodrick Date: Wed, 11 Dec 2024 11:18:08 -0800 Subject: [PATCH 03/14] updates for functional testing --- build_luts.py | 84 +++++++++++++++++++++++++++++++++++---------------- 1 file changed, 58 insertions(+), 26 deletions(-) diff --git a/build_luts.py b/build_luts.py index 7b899b2..2b62217 100644 --- a/build_luts.py +++ b/build_luts.py @@ -25,6 +25,7 @@ import argparse from types import SimpleNamespace import logging +import subprocess from scipy.stats import qmc @@ -69,6 +70,11 @@ def main(): default=None, help='MODTRAN Installation DIR' ) + parser.add_argument( + '--coarse', + default=None, + type=str + ) args = parser.parse_args() args.training = args.train == 1 @@ -115,14 +121,26 @@ def main(): else: if args.train: grid = { - '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], + #'to_solar_zenith_lut': [0, 12.5, 25, 37.5, 50], + #'to_solar_azimuth_lut': [0, 60, 120, 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], + #'h2o_lut': list(np.sort( + # np.round(np.linspace(0.1, 5, num=5), 3).tolist() + [0.6125] + #)), + #'aerfrac_2_lut': list(np.sort( + # np.round(np.linspace(0.01, 1, num=5), 3).tolist() + [0.5] + #)) + '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=5), 3).tolist() + [0.6125] + 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] @@ -141,8 +159,8 @@ def main(): # HOLDOUT SET grid = { 'to_solar_zenith_lut': [6, 18, 30,45], - 'to_solar_azimuth_lut': [60, 300], - 'to_sensor_azimuth_lut': [180], + '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], @@ -174,16 +192,19 @@ def main(): 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( - paths.wavelength_file, - wl_file_contents,fmt='%.5f' - ) + if args.coarse is None: + 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( + 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], @@ -346,10 +367,21 @@ def __init__(self, args, training=True): self.support_dir, 'aerosol_model.txt' ) - self.wavelength_file = os.path.join( - self.support_dir, - 'hifidelity_wavelengths.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") @@ -359,7 +391,7 @@ def __init__(self, args, training=True): ) self.irradiance_file = os.path.join( self.support_dir, - 'prism_optimized_irr.dat' + 'kurudz_0.1nm.dat' ) if training: @@ -439,7 +471,7 @@ def build_main_config( paths.lut_modtran_directory, 'lut.nc' ), - "multipart_transmittance": False, + "multipart_transmittance": True, "template_file": paths.modtran_template_file, "rte_configure_and_exit": configure_and_exit, "engine_base_dir": paths.modtran_path, @@ -463,7 +495,7 @@ def build_main_config( "viewzen": 180 - to_sensor_zenith_lut_grid[0], "solaz": to_solar_azimuth_lut_grid[0], "solzen": to_solar_zenith_lut_grid[0], - "multipart_transmittance": False, + "multipart_transmittance": True, "template_file": paths.modtran_template_file, "rte_configure_and_exit": configure_and_exit, "engine_base_dir": paths.sixs_path, From 6bc508221cdeffe8558977b31c281b794bb721ba Mon Sep 17 00:00:00 2001 From: brodrick Date: Fri, 13 Dec 2024 14:58:42 -0800 Subject: [PATCH 04/14] training setup, read rework --- aggregated_combined_data.py | 98 +++++++++++++++++++++++++++++-------- build_luts.py | 40 +++++++++------ 2 files changed, 101 insertions(+), 37 deletions(-) diff --git a/aggregated_combined_data.py b/aggregated_combined_data.py index bc14810..f7aa0a5 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 @@ -45,7 +46,7 @@ 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') args = parser.parse_args() @@ -59,10 +60,65 @@ def main(): 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) + + rt_config = config.forward_model.radiative_transfer + instrument_config = config.forward_model.instrument + + + 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]) + 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()} + params['engine_config'].rte_configure_and_exit = False + params['engine_config'].rt_mode = 'rdn' + params['wavelength_file'] = confPriority('wavelength_file', [params['engine_config'], instrument_config, rt_config]) + isofit_sixs = SixSRT(**params) + + + point_dims = list(isofit_modtran.lut_grid.keys()) + rtm_key = 'transm_down_dir' + for sp,mp,std_dir,mtd_dir, std_dif, mtd_dif, stu_dir, mtu_dir, stu_dif, mtu_dif, s_r, m_r in zip(isofit_sixs.points, isofit_modtran.points, + isofit_sixs.lut['transm_down_dir'], isofit_modtran.lut['transm_down_dir'], + isofit_sixs.lut['transm_down_dif'], isofit_modtran.lut['transm_down_dif'], + isofit_sixs.lut['transm_up_dir'], isofit_modtran.lut['transm_up_dir'], + isofit_sixs.lut['transm_up_dif'], isofit_modtran.lut['transm_up_dif'], + isofit_sixs.lut['rhoatm'], isofit_modtran.lut['rhoatm'], + ): + name='_'.join([f'{point_dims[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'figs/{name.replace("\n","_")}.png',dpi=100) + plt.clf() + + exit() + #import ipdb; ipdb.set_trace() + #sixs_results = 1 + #modtran_results = 1 + sixs_results = get_obj_res(isofit_sixs, key, resample=False) modtran_results = get_obj_res(isofit_modtran, key) @@ -108,25 +164,25 @@ def main(): 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' + #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')] + #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])])) @@ -176,7 +232,7 @@ def read_data_piece(ind, maxind, point, fn, key, resample, obj): def get_obj_res(obj, key, resample=True): # We don't want the VectorInterpolator, but rather the raw inputs - ray.init(temp_dir='/tmp/ray/brodrick/') + ray.init(_temp_dir='/local/ray/brodrick/') if hasattr(obj,'sixs_ngrid_init'): results = np.zeros((obj.points.shape[0],obj.sixs_ngrid_init), dtype=float) else: diff --git a/build_luts.py b/build_luts.py index 2b62217..924ed92 100644 --- a/build_luts.py +++ b/build_luts.py @@ -122,29 +122,32 @@ def main(): 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_sensor_azimuth_lut': [0], - #'to_sensor_zenith_lut': [140, 160, 180], - #'altitude_km_lut': [2, 4, 7, 10, 15, 25, 99], + '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() + )), + '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=5), 3).tolist() + [0.6125] + # 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] #)) - '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']) @@ -299,6 +302,7 @@ def main(): for key in _keys } params["engine_config"] = modtran_engine_config + params["build_interpolator"] = False params['lut_grid'] = { key: params['lut_grid'][key] for key in modtran_engine_config.lut_names.keys() @@ -313,6 +317,7 @@ def main(): for key in _keys } params["engine_config"] = sixs_engine_config + params["build_interpolator"] = False params['lut_grid'] = { key: params['lut_grid'][key] for key in sixs_engine_config.lut_names.keys() @@ -471,6 +476,8 @@ def build_main_config( 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, @@ -480,6 +487,7 @@ def build_main_config( }, "sixs": { "engine_name": '6s', + "rt_mode": 'rdn', "sim_path": paths.lut_sixs_directory, "lut_path": os.path.join( paths.lut_sixs_directory, From d8c33adb84b73742882ed6c76ef4c3f8facd77e7 Mon Sep 17 00:00:00 2001 From: brodrick Date: Sun, 15 Dec 2024 17:52:27 -0800 Subject: [PATCH 05/14] updates, sims and data aggregation working --- aggregated_combined_data.py | 222 +++++++++--------------------------- build_luts.py | 24 ++-- 2 files changed, 70 insertions(+), 176 deletions(-) diff --git a/aggregated_combined_data.py b/aggregated_combined_data.py index f7aa0a5..e5ee847 100644 --- a/aggregated_combined_data.py +++ b/aggregated_combined_data.py @@ -48,13 +48,15 @@ def main(): parser.add_argument('-config_file', type=str, default='templates/isofit_template.json') 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) args = parser.parse_args() np.random.seed(13) + outdict_sixs, outdict_modtran = {}, {} + munge_file = os.path.join(args.munge_dir, 'combined_data.npz') 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) @@ -84,174 +86,60 @@ def main(): isofit_sixs = SixSRT(**params) - point_dims = list(isofit_modtran.lut_grid.keys()) - rtm_key = 'transm_down_dir' - for sp,mp,std_dir,mtd_dir, std_dif, mtd_dif, stu_dir, mtu_dir, stu_dif, mtu_dif, s_r, m_r in zip(isofit_sixs.points, isofit_modtran.points, - isofit_sixs.lut['transm_down_dir'], isofit_modtran.lut['transm_down_dir'], - isofit_sixs.lut['transm_down_dif'], isofit_modtran.lut['transm_down_dif'], - isofit_sixs.lut['transm_up_dir'], isofit_modtran.lut['transm_up_dir'], - isofit_sixs.lut['transm_up_dif'], isofit_modtran.lut['transm_up_dif'], - isofit_sixs.lut['rhoatm'], isofit_modtran.lut['rhoatm'], - ): - name='_'.join([f'{point_dims[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'figs/{name.replace("\n","_")}.png',dpi=100) - plt.clf() - - exit() - #import ipdb; ipdb.set_trace() - #sixs_results = 1 - #modtran_results = 1 - - - 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,:] - - - good_data = np.all(np.isnan(modtran_results) == False,axis=1) - good_data[np.any(np.isnan(sixs_results),axis=1)] = False - - modtran_results = modtran_results[good_data,:] - sixs_results = sixs_results[good_data,:] - points = points[good_data,...] - points_sixs = points_sixs[good_data,...] - - 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 - - + 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[key][:, point_names.index('surface_elevation_km')] >= isofit_sixs.lut[key][:, point_names.index('observer_altitude_km')] -2 + 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,:] -def get_obj_res(obj, key, resample=True): + np.savez(munge_file, modtran_results=outdict_modtran, + sixs_results=outdict_sixs, + sol_irr=isofit_sixs.lut.solar_irr) - # We don't want the VectorInterpolator, but rather the raw inputs - ray.init(_temp_dir='/local/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 924ed92..04e39f9 100644 --- a/build_luts.py +++ b/build_luts.py @@ -227,9 +227,15 @@ def main(): 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 + 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, @@ -302,7 +308,7 @@ def main(): for key in _keys } params["engine_config"] = modtran_engine_config - params["build_interpolator"] = False + params["build_interpolators"] = False params['lut_grid'] = { key: params['lut_grid'][key] for key in modtran_engine_config.lut_names.keys() @@ -317,7 +323,7 @@ def main(): for key in _keys } params["engine_config"] = sixs_engine_config - params["build_interpolator"] = False + params["build_interpolators"] = False params['lut_grid'] = { key: params['lut_grid'][key] for key in sixs_engine_config.lut_names.keys() @@ -581,9 +587,9 @@ def build_main_config( 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']['AOT550'] = [ + float(q) for q in aerfrac_2_lut_grid + ] if paths.aerosol_model_path is not None: radiative_transfer_config[ @@ -615,7 +621,7 @@ def build_main_config( 'H2OSTR','surface_elevation_km','observer_altitude_km', # 'viewzen','solar_zenith', 'relative_azimuth', 'observer_zenith','solar_zenith', 'relative_azimuth', - 'AERFRAC_2' + 'AOT550' ] if x in radiative_transfer_config['lut_grid'].keys() } radiative_transfer_config[ From 19e8cb7e1a90cad147a33137ac6208454d6ff7fe Mon Sep 17 00:00:00 2001 From: Phil Brodrick Date: Sun, 15 Dec 2024 18:13:05 -0800 Subject: [PATCH 06/14] add basic plotting - probably to remove and place in isofit-plots --- plot_munged.py | 147 +++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 147 insertions(+) create mode 100644 plot_munged.py diff --git a/plot_munged.py b/plot_munged.py new file mode 100644 index 0000000..b64a085 --- /dev/null +++ b/plot_munged.py @@ -0,0 +1,147 @@ + + + + +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')] -2 + good_points = np.logical_not(bad_points) + bad_ind = np.any(rtm_comp['transm_down_dif'].data[:,:] > 10,axis=1) + import ipdb; ipdb.set_trace() + + 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.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, linestyle='--') + + plt.savefig(f'{args.fig_dir}/dim_{point_names[dim]}.png', dpi=200, bbox_inches='tight') + plt.clf() + + +if __name__ == '__main__': + main() From f1c7682090113d79b1bf77e55424e13f420eb5a9 Mon Sep 17 00:00:00 2001 From: Evan Greenberg Date: Thu, 19 Dec 2024 10:37:21 -0800 Subject: [PATCH 07/14] add working sobel grid --- build_luts.py | 116 ++++++++++++++++++++++++++++---------------------- 1 file changed, 65 insertions(+), 51 deletions(-) diff --git a/build_luts.py b/build_luts.py index 04e39f9..28f5681 100644 --- a/build_luts.py +++ b/build_luts.py @@ -91,33 +91,56 @@ def main(): if args.train: # Get bounds from regular grid bounds = { - 'to_solar_zenith_bnds': [0, 50, 2], - 'to_sensor_zenith_bnds': [140, 180, 2], - 'altitude_km_bnds': [2, 25, 2], - 'elevation_km_bnds': [0.01, 4.5, 2], - 'h2o_bnds': [0.1, 5, 3], - 'aerfrac_2_bnds': [0.01, 1, 3], + '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, } - grid = sobel_lut(bounds, consts) + # 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, 3], - 'to_sensor_zenith_bnds': [140, 180, 3], - 'altitude_km_bnds': [2, 25, 3], - 'elevation_km_bnds': [0.01, 4.5, 3], - 'h2o_bnds': [0.1, 5, 5], - 'aerfrac_2_bnds': [0.01, 1, 5], + '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, } - rid = sobel_lut(bounds, consts) + 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: if args.train: grid = { @@ -125,7 +148,6 @@ def main(): '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], @@ -207,7 +229,6 @@ def main(): wl_file_contents,fmt='%.5f' ) - write_modtran_template( atmosphere_type='ATM_MIDLAT_SUMMER', fid=os.path.splitext(paths.modtran_template_file)[0], @@ -215,7 +236,7 @@ def main(): 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_azimuth_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], @@ -662,7 +683,7 @@ def build_main_config( )) -def sobel_lut(bounds, consts={}, log=False): +def sobel_lut(bounds, consts={}, n=3, log=False): """ Create a look using a Sobel sequence @@ -686,7 +707,7 @@ def sobel_lut(bounds, consts={}, log=False): ] for key, bnd in bounds.items(): if key in cos_keys: - bnd = [cosd(bnd[0]), cosd(bnd[1]), bnd[2]] + bnd = [cosd(bnd[0]), cosd(bnd[1])] bounds[key] = bnd # Handle input log transform @@ -699,54 +720,47 @@ def sobel_lut(bounds, consts={}, log=False): if log: for key, bnd in bounds.items(): if key in log_keys: - bnd = [np.log(bnd[0]), np.log(bnd[1]), bnd[2]] + bnd = [np.log(bnd[0]), np.log(bnd[1])] bounds[key] = bnd # Sort bounds for key, bnd in bounds.items(): - bounds[key] = sorted(bnd[:-1]) + [bnd[2]] + 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 - grid = {} - for key, value in bounds.items(): - l_bound, r_bound = value[:-1] - - # Not 100% sure if we need to scamble the sequence gen - sobol = qmc.Sobol(d=1, scramble=False) - seq = sobol.random(n=value[-1]) - seq = qmc.scale( - seq, - l_bounds=l_bound, - u_bounds=r_bound - ) + # 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 - if key in cos_keys: - seq = np.rad2deg(np.arccos(seq)) + # 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: - if key in log_keys: - seq = np.exp(seq) + # 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] + 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)] - # 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 - )] - # Trim out non-physical components # (e.g., elevation > observation altitude) # The best way I can do this with uneven clolumns is to remove From e226795e166b2d5e32b5d95b1b272fd1ddfa0c80 Mon Sep 17 00:00:00 2001 From: Evan Greenberg Date: Thu, 19 Dec 2024 10:39:08 -0800 Subject: [PATCH 08/14] missing lut element in regular training --- build_luts.py | 1 + 1 file changed, 1 insertion(+) diff --git a/build_luts.py b/build_luts.py index 28f5681..3248845 100644 --- a/build_luts.py +++ b/build_luts.py @@ -148,6 +148,7 @@ def main(): '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], From 4772d2b3f3e10414d98dcb514db4560597baa123 Mon Sep 17 00:00:00 2001 From: brodrick Date: Fri, 27 Dec 2024 03:44:50 -0800 Subject: [PATCH 09/14] training and aggregation updates - functional if not optimal --- aggregated_combined_data.py | 87 ++++++++++++++-------- train_combined_emulator.py | 141 ++++++++++++++++++++++++++++-------- 2 files changed, 167 insertions(+), 61 deletions(-) diff --git a/aggregated_combined_data.py b/aggregated_combined_data.py index e5ee847..1478a40 100644 --- a/aggregated_combined_data.py +++ b/aggregated_combined_data.py @@ -54,38 +54,45 @@ def main(): np.random.seed(13) - outdict_sixs, outdict_modtran = {}, {} - munge_file = os.path.join(args.munge_dir, 'combined_data.npz') - for key_ind, key in enumerate(args.keys): - 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 - rt_config = config.forward_model.radiative_transfer - instrument_config = config.forward_model.instrument + config = configs.create_new_config(args.config_file) - - 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]) - 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()} - params['engine_config'].rte_configure_and_exit = False - params['engine_config'].rt_mode = 'rdn' - params['wavelength_file'] = confPriority('wavelength_file', [params['engine_config'], instrument_config, rt_config]) - isofit_sixs = SixSRT(**params) + # 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 + + + 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]) + 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()} + params['engine_config'].rte_configure_and_exit = False + params['engine_config'].rt_mode = 'rdn' + params['wavelength_file'] = confPriority('wavelength_file', [params['engine_config'], instrument_config, rt_config]) + isofit_sixs = SixSRT(**params) + + + + + 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): + if args.figs_dir is not None: point_dims = list(isofit_modtran.lut_grid.keys()) point_dims_s = list(isofit_sixs.lut_grid.keys()) @@ -130,15 +137,35 @@ def main(): 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[key][:, point_names.index('surface_elevation_km')] >= isofit_sixs.lut[key][:, point_names.index('observer_altitude_km')] -2 + 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,:] - np.savez(munge_file, modtran_results=outdict_modtran, - sixs_results=outdict_sixs, - sol_irr=isofit_sixs.lut.solar_irr) + #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 + ) diff --git a/train_combined_emulator.py b/train_combined_emulator.py index 037fa58..69d3291 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 @@ -23,24 +23,76 @@ import pickle +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,6 +112,7 @@ 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 @@ -118,7 +171,7 @@ def main(): np.random.seed(13) - npzf = np.load(args.munged_file) + npzf = np.load(args.munged_file, allow_pickle=True) modtran_results = npzf['modtran_results'] sixs_results = npzf['sixs_results'] @@ -134,6 +187,14 @@ 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 + + print(modtran_results.shape) + + if args.holdout_dim == -1: perm = np.random.permutation(points.shape[0]) @@ -174,6 +235,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') @@ -189,8 +252,11 @@ def main(): 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) + simple_response_scaler = np.ones(train_modtran.shape[1])*100 + simple_response_scaler = 1/np.mean(train_modtran,axis=0) 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 +266,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 +285,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 +327,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 +353,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 +377,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 +463,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 +476,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()): From 2ad8f14b377a7f42db11dd1bc85d90f09f2a4304 Mon Sep 17 00:00:00 2001 From: brodrick Date: Sat, 28 Dec 2024 04:42:55 -0800 Subject: [PATCH 10/14] expand WV range for training RG --- build_luts.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/build_luts.py b/build_luts.py index 04e39f9..9a8fc0f 100644 --- a/build_luts.py +++ b/build_luts.py @@ -132,7 +132,7 @@ def main(): '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() )) From b6c215ed70a51da6670ba637f17ad303a317d3ac Mon Sep 17 00:00:00 2001 From: brodrick Date: Tue, 31 Dec 2024 11:55:11 -0800 Subject: [PATCH 11/14] integrate hires changes, training still under dev --- aggregated_combined_data.py | 72 ++++++++++++++++++++++++++++++++++--- build_luts.py | 13 +++++-- train_combined_emulator.py | 64 ++++++++++++++++++++++++++++++--- 3 files changed, 137 insertions(+), 12 deletions(-) diff --git a/aggregated_combined_data.py b/aggregated_combined_data.py index 1478a40..ba2f167 100644 --- a/aggregated_combined_data.py +++ b/aggregated_combined_data.py @@ -40,6 +40,36 @@ 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(): @@ -49,6 +79,7 @@ def main(): 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() @@ -70,17 +101,29 @@ def main(): 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]) - params['engine_config'].rte_configure_and_exit = False - params['engine_config'].rt_mode = 'rdn' + 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()} - params['engine_config'].rte_configure_and_exit = False - params['engine_config'].rt_mode = 'rdn' - 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' + + # 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) @@ -145,6 +188,25 @@ def main(): 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) diff --git a/build_luts.py b/build_luts.py index 9a8fc0f..c1671ba 100644 --- a/build_luts.py +++ b/build_luts.py @@ -132,7 +132,7 @@ def main(): '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], + )), #+ [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() )) @@ -193,14 +193,15 @@ def main(): 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 if args.coarse is None: - wl = np.arange(0.350, 2.550, 0.0005) + 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.0005 + wl_file_contents[:,2] = 0.0001 np.savetxt( paths.wavelength_file, @@ -329,6 +330,12 @@ def main(): sixs_engine_config.lut_names.keys() } # params['modtran_emulation'] = True + + # 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 diff --git a/train_combined_emulator.py b/train_combined_emulator.py index 69d3291..db3ed8c 100644 --- a/train_combined_emulator.py +++ b/train_combined_emulator.py @@ -22,6 +22,8 @@ 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 @@ -119,6 +121,35 @@ def nn_model(in_data_shape, out_data_shape, num_layers=5): +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): @@ -170,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, allow_pickle=True) - modtran_results = npzf['modtran_results'] sixs_results = npzf['sixs_results'] + modtran_results = npzf['modtran_results'] points = npzf['points'] keys = npzf['keys'] point_names = npzf['point_names'] @@ -192,6 +234,20 @@ def main(): 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) @@ -250,14 +306,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)) + print(model.summary()) simple_response_scaler = np.ones(train_modtran.shape[1])*100 - simple_response_scaler = 1/np.mean(train_modtran,axis=0) train_modtran *= simple_response_scaler #import ipdb; ipdb.set_trace() - model.fit(train_sixs[train,:], train_modtran[train,:], batch_size=1000, epochs=400, + model.fit(train_sixs[train,:], train_modtran[train,:], batch_size=10, epochs=400, validation_data=(train_sixs[test,:], train_modtran[test,:]),callbacks=[es]) train_modtran /= simple_response_scaler From 8df99870f0e01f57d01377b9c69ee4ec39a11121 Mon Sep 17 00:00:00 2001 From: Phil Brodrick Date: Tue, 31 Dec 2024 11:56:46 -0800 Subject: [PATCH 12/14] small plotting updates --- plot_munged.py | 27 +++++++++++++++++++++++---- 1 file changed, 23 insertions(+), 4 deletions(-) diff --git a/plot_munged.py b/plot_munged.py index b64a085..56c5966 100644 --- a/plot_munged.py +++ b/plot_munged.py @@ -69,10 +69,27 @@ def main(): 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')] -2 - good_points = np.logical_not(bad_points) + 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) - import ipdb; ipdb.set_trace() + + 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]] @@ -100,6 +117,7 @@ def main(): 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]') @@ -137,7 +155,8 @@ def main(): 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, linestyle='--') + 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() From 36469fbef6448b9ac4ef54b325816a32b7eefcc1 Mon Sep 17 00:00:00 2001 From: Phil Brodrick Date: Tue, 31 Dec 2024 11:57:27 -0800 Subject: [PATCH 13/14] add basic convolve script, untested --- convolve.py | 443 ++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 443 insertions(+) create mode 100644 convolve.py diff --git a/convolve.py b/convolve.py new file mode 100644 index 0000000..cf7e9b1 --- /dev/null +++ b/convolve.py @@ -0,0 +1,443 @@ + + + + + + + +import numpy as np +import argparse +import h5py +import subprocess +import os +import ray +import multiprocessing as mp + + +class resample(): + """ Resampling class where you can initialize and then use. It has a saftey logic against non finite + elements in the data to interpoloate. + + Args: + wvs_a: High resolution data's spectral grid + wvs_b: Instrument's resolution spectra grid + fwhm_b: Instrument's FWHM + """ + def __init__(self, wvs_a, wvs_b, fwhm_b): + self.wvs_a = wvs_a + self.wvs_b = wvs_b + self.fwhm_b = fwhm_b + self.get_transform_matrix() + + + def get_transform_matrix(self): + base_wl = np.array(self.wvs_a) + self.base_wl = base_wl + target_wl = np.array(self.wvs_b) + self.target_wl = target_wl + target_fwhm = np.array(self.fwhm_b) + + doTheResample = lambda id: self.spectrumResample(id, base_wl, target_wl, target_fwhm) + ww = np.array([doTheResample(W) for W in range(len(target_wl))]) + self.transform_matrix = ww + + + def process_data(self, ys, n_samples, num_processes=4): + """ + Process data in parallel. + Args: + ys: The data to be processed. + n_samples: The number of samples expected in ys. + num_processes: The number of parallel processes to use. + Returns: + Processed data. + """ + num_processes = mp.cpu_count() + # Check dimensions and align correctly + ys = self._reshape_data(ys, n_samples) + + # Initialize Ray + #if not ray.is_initialized(): + # ray.init(num_cpus=num_processes) + + # Dispatch tasks + #result_ids = [self.process_single_sample.remote(self, ys[i, :]) for i in range(ys.shape[0])] + #results = ray.get(result_ids) + results = [self.process_single_sample(ys[i, :],i) for i in range(ys.shape[0])] + + return np.array(results) + + #@ray.remote(num_cpus=1) + def process_single_sample(self, y, i): + """ + Process a single sample. This method will be called in parallel. + Args: + y: A single sample from ys. + Returns: + The processed sample. + """ + # Processing logic for a single sample + # For example, this could be a call to self.__call__(y) or other processing + if i % 100: + print(i) + return self.__call__(y) # or other processing logic + + + def _reshape_data(self, ys, n_samples): + """ + Reshape data to the correct dimensions for convolution. + """ + counter = 0 + max_attempts = 5 + while ys.shape[0] != n_samples: + if counter >= max_attempts: + raise ValueError("Unable to reshape data to the correct dimensions for convolution") + ys = ys.T + counter += 1 + return ys + + def __call__(self, y): + # Convert input to 2D array and transpose if necessary + + # Convert input to 2D array + spectrum = np.atleast_2d(y) + + # Check if transpose is necessary + transpose_needed = spectrum.shape[0] == 1 + + # Initialize an output array + resampled_spectrum = np.zeros(self.transform_matrix.shape[0]) + + # Identify valid (non-NaN, non-inf, non--inf) elements in the spectrum + valid_indices = np.isfinite(spectrum[0]) if transpose_needed else np.isfinite(spectrum[:, 0]) + + # Optimize the loop with vectorized operations + if transpose_needed: + spectrum = spectrum.T + + resampled_spectrum = np.dot(self.transform_matrix[:, valid_indices], spectrum[valid_indices]) + + return np.squeeze(resampled_spectrum) + + + def srf(self, x, mu, sigma): + """Spectral Response Function """ + u = (x-mu)/abs(sigma) + y = (1.0/(np.sqrt(2.0*np.pi)*abs(sigma)))*np.exp(-u*u/2.0) + if y.sum()==0: + return y + else: + return y/y.sum() + + + def spectrumResample(self, idx, wl, wl2, fwhm2=10, fill=False): + """Resample a spectrum to a new wavelength / FWHM. + I assume Gaussian SRFs""" + + resampled = np.array(self.srf(wl, wl2[idx], fwhm2[idx]/2.35482)) + + + return resampled + + +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 + if args.output_sixs_wl_file is not None: + + npzf = np.load(args.input_file) + input_sixs_wl = npzf['sixs_wavelengths'] + + 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 + + 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 + + input_sixs = npzf['sixs_results'].copy() + input_solar_irr = npzf['solar_irradiance'] + points = npzf['points'] + point_names = npzf['point_names'] + + 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, 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')])) + + output_sixs = np.dot(input_sixs, H.T) + print(output_sixs.shape) + else: + output_sixs = npzf['sixs_results'] + output_sixs_wl = npzf['sixs_wavelengths'] + + output_modtran = None + output_solar_irr = None + output_modtran_wl = None + if args.output_modtran_wl_file is not None: + npzf = np.load(args.input_file) + input_modtran_wl = npzf['modtran_wavelengths'] + + 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 + + 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 + + input_modtran = npzf['modtran_results'].copy() + input_solar_irr = npzf['solar_irradiance'] + points = npzf['points'] + + 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, input_modtran_wl[1] - input_modtran_wl[0]),dtype=np.float32) + output_modtran_irr = np.array(resample_spectrum(irr, irr_wl, output_modtran_wl, output_modtran_fwhm),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')])) + + output_modtran = np.dot(input_modtran, H.T) + print(output_modtran.shape) + + + else: + output_modtran = npzf['modtran_results'] + 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() +#input_file = "/Users/brodrick/Downloads/emit_emulator_v1.hdf5" +##out_file = "/Users/brodrick/Downloads/emit_emulator_v1_conv01.hdf5" +#out_file = "/Users/brodrick/Downloads/emit_emulator_v1_conv_EMIT_Wavelengths_20220817_clipped.hdf5" +#out_file = "/Users/brodrick/repos/kf/convolved/neon_emulator_v1_conv.hdf5" +#subprocess.call(f"cp {input_file} {out_file}", shell=True) + +input_file = "data_202411121/emit_emulator_v2_0.hdf5" +out_file = "data_202411121/emit_emulator_v2_conv_EMIT_Wavelengths_20220817_clipped.hdf5" + +input = h5py.File(input_file, "r") +output = h5py.File(out_file, "w") + +def copy_groups(name, obj): + if isinstance(obj, h5py.Group): + output.create_group(name) +input.visititems(copy_groups) + +def copy_attributes(name, obj): + if isinstance(obj, h5py.Group) or isinstance(obj, h5py.Dataset): + for key, value in obj.attrs.items(): + output[name].attrs[key] = value +input.visititems(copy_attributes) + + +wl = input["MISCELLANEOUS"]["Wavelengths"][:] +X = input["sample_space"]["sample space"][:,:] +Y_all = input["mod_output/modtran output"][:,:,:] + + +output_wl = np.genfromtxt('wl/EMIT_Wavelengths_20220817.txt')[:,1:][::-1,:]*1000 +subset = np.where((output_wl[:,0] >= wl[0]) & (output_wl[:,0] <= wl[-1]))[0] +output_wl = output_wl[subset,:] +subset = np.where((output_wl[:,0] >= 380) & (output_wl[:,0] <= 2493))[0] +output_wl = output_wl[subset,:] +#output_wl = np.genfromtxt('wl/neon_wl.txt')[:,1:] + + +print('prepping') +for key, value in input["sample_space"].items(): + output.create_dataset(f"sample_space/{key}", data=value) +for key, value in input["MISCELLANEOUS"].items(): + output.create_dataset(f"MISCELLANEOUS/{key}", data=value) + +# overwrite wavelengths +del output["MISCELLANEOUS/Wavelengths"] +output["MISCELLANEOUS/Wavelengths"] = output_wl[:,0] + +yn = input['mod_output'].attrs['Products'].tolist() +del input +print('convolving') +##res = resample(wl, np.arange(wl[0], wl[-1], 0.1), 0.1*np.ones(len(np.arange(wl[0], wl[-1], 0.1)))) +##res = resample(wl, np.arange(wl[0], wl[-1], 5.0), 5.0*np.ones(len(np.arange(wl[0], wl[-1], 5.0)))) +print(output_wl[:,0].shape) +print(Y_all.shape) +print(wl.shape) +res = resample(wl, output_wl[:,0], output_wl[:,1]) + +yn_out = [] +Y_out = [] +yn_out.append('path_radiance') +Y_out.append(res.process_data(Y_all[:,yn.index('path_radiance'),:], Y_all.shape[0])) + +#yn_out.append('t_down_dir_t_up_dir') +yn_out.append('bi-direct') +Y_out.append(res.process_data(Y_all[:,yn.index('t_down_dir'),:]*\ + Y_all[:,yn.index('t_up_dir'),:]*\ + Y_all[:,yn.index('ToA_irrad'),:], Y_all.shape[0])) + +#yn_out.append('t_down_dif_t_up_dir') +yn_out.append('hemi-direct') +Y_out.append(res.process_data(Y_all[:,yn.index('t_down_dif'),:]*\ + Y_all[:,yn.index('t_up_dir'),:]*\ + Y_all[:,yn.index('ToA_irrad'),:], Y_all.shape[0])) + +#yn_out.append('t_down_dir_t_up_dif') +yn_out.append('direct-hemi') +Y_out.append(res.process_data(Y_all[:,yn.index('t_down_dir'),:]*\ + Y_all[:,yn.index('t_up_dif'),:]*\ + Y_all[:,yn.index('ToA_irrad'),:], Y_all.shape[0])) + +#yn_out.append('t_down_dif_t_up_dif') +yn_out.append('bi-hemi') +Y_out.append(res.process_data(Y_all[:,yn.index('t_down_dif'),:]*\ + Y_all[:,yn.index('t_up_dif'),:]*\ + Y_all[:,yn.index('ToA_irrad'),:], Y_all.shape[0])) + +yn_out.append('sphalb') +Y_out.append(res.process_data(Y_all[:,yn.index('sphalb_num'),:]/\ + Y_all[:,yn.index('sphalb_denom'),:], Y_all.shape[0])) + +Y_out = np.stack(Y_out, axis=1) + + +print('writing') +output["mod_output/modtran output"] = Y_out +output['mod_output'].attrs['Products'] = yn_out +del output + +print('checking') +input = h5py.File(out_file, "r") +Y_all = input["mod_output/modtran output"][:,:,:] +print(Y_all.shape) +del input + + +#Y_out = [] +#for n in range(Y_all.shape[1]): +# Y_out.append(res.process_data(Y_all[:,n,:], Y_all.shape[0])) +# +#Y_out = np.stack(Y_out, axis=1) +# +#del input["mod_output/modtran output"] +#input["mod_output/modtran output"] = Y_out +#del input["MISCELLANEOUS/Wavelengths"] +#input["MISCELLANEOUS/Wavelengths"] = output_wl[:,0] +#del input +#input = h5py.File(out_file, "r") +#Y_all = input["mod_output/modtran output"][:,:,:] +#print(Y_all.shape) + + + + + + + + + From 984212dd580a1512db2dca6e02675ca3999ff657 Mon Sep 17 00:00:00 2001 From: brodrick Date: Fri, 3 Jan 2025 03:47:06 -0800 Subject: [PATCH 14/14] convolve script fixed and trimmed, model training reset to use convolved (temporary) --- convolve.py | 337 ++++++------------------------------- train_combined_emulator.py | 14 +- 2 files changed, 66 insertions(+), 285 deletions(-) diff --git a/convolve.py b/convolve.py index cf7e9b1..fc708a2 100644 --- a/convolve.py +++ b/convolve.py @@ -14,131 +14,6 @@ import multiprocessing as mp -class resample(): - """ Resampling class where you can initialize and then use. It has a saftey logic against non finite - elements in the data to interpoloate. - - Args: - wvs_a: High resolution data's spectral grid - wvs_b: Instrument's resolution spectra grid - fwhm_b: Instrument's FWHM - """ - def __init__(self, wvs_a, wvs_b, fwhm_b): - self.wvs_a = wvs_a - self.wvs_b = wvs_b - self.fwhm_b = fwhm_b - self.get_transform_matrix() - - - def get_transform_matrix(self): - base_wl = np.array(self.wvs_a) - self.base_wl = base_wl - target_wl = np.array(self.wvs_b) - self.target_wl = target_wl - target_fwhm = np.array(self.fwhm_b) - - doTheResample = lambda id: self.spectrumResample(id, base_wl, target_wl, target_fwhm) - ww = np.array([doTheResample(W) for W in range(len(target_wl))]) - self.transform_matrix = ww - - - def process_data(self, ys, n_samples, num_processes=4): - """ - Process data in parallel. - Args: - ys: The data to be processed. - n_samples: The number of samples expected in ys. - num_processes: The number of parallel processes to use. - Returns: - Processed data. - """ - num_processes = mp.cpu_count() - # Check dimensions and align correctly - ys = self._reshape_data(ys, n_samples) - - # Initialize Ray - #if not ray.is_initialized(): - # ray.init(num_cpus=num_processes) - - # Dispatch tasks - #result_ids = [self.process_single_sample.remote(self, ys[i, :]) for i in range(ys.shape[0])] - #results = ray.get(result_ids) - results = [self.process_single_sample(ys[i, :],i) for i in range(ys.shape[0])] - - return np.array(results) - - #@ray.remote(num_cpus=1) - def process_single_sample(self, y, i): - """ - Process a single sample. This method will be called in parallel. - Args: - y: A single sample from ys. - Returns: - The processed sample. - """ - # Processing logic for a single sample - # For example, this could be a call to self.__call__(y) or other processing - if i % 100: - print(i) - return self.__call__(y) # or other processing logic - - - def _reshape_data(self, ys, n_samples): - """ - Reshape data to the correct dimensions for convolution. - """ - counter = 0 - max_attempts = 5 - while ys.shape[0] != n_samples: - if counter >= max_attempts: - raise ValueError("Unable to reshape data to the correct dimensions for convolution") - ys = ys.T - counter += 1 - return ys - - def __call__(self, y): - # Convert input to 2D array and transpose if necessary - - # Convert input to 2D array - spectrum = np.atleast_2d(y) - - # Check if transpose is necessary - transpose_needed = spectrum.shape[0] == 1 - - # Initialize an output array - resampled_spectrum = np.zeros(self.transform_matrix.shape[0]) - - # Identify valid (non-NaN, non-inf, non--inf) elements in the spectrum - valid_indices = np.isfinite(spectrum[0]) if transpose_needed else np.isfinite(spectrum[:, 0]) - - # Optimize the loop with vectorized operations - if transpose_needed: - spectrum = spectrum.T - - resampled_spectrum = np.dot(self.transform_matrix[:, valid_indices], spectrum[valid_indices]) - - return np.squeeze(resampled_spectrum) - - - def srf(self, x, mu, sigma): - """Spectral Response Function """ - u = (x-mu)/abs(sigma) - y = (1.0/(np.sqrt(2.0*np.pi)*abs(sigma)))*np.exp(-u*u/2.0) - if y.sum()==0: - return y - else: - return y/y.sum() - - - def spectrumResample(self, idx, wl, wl2, fwhm2=10, fill=False): - """Resample a spectrum to a new wavelength / FWHM. - I assume Gaussian SRFs""" - - resampled = np.array(self.srf(wl, wl2[idx], fwhm2[idx]/2.35482)) - - - return resampled - def spectral_response_function(response_range: np.array, mu: float, sigma: float): """Calculate the spectral response function. @@ -223,46 +98,69 @@ def main(): output_sixs = None output_sixs_wl = None - if args.output_sixs_wl_file is not None: + npzf = np.load(args.input_file) - npzf = np.load(args.input_file) - input_sixs_wl = npzf['sixs_wavelengths'] - 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 - - 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 + # SIXS + input_sixs_wl = npzf['sixs_wavelengths'] + - input_sixs = npzf['sixs_results'].copy() - input_solar_irr = npzf['solar_irradiance'] - points = npzf['points'] - point_names = npzf['point_names'] + 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, input_sixs_wl[1]-input_sixs_wl[0]),dtype=np.float32) + 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')])) + 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 = np.dot(input_sixs, H.T) + 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 = npzf['sixs_results'] + 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: - npzf = np.load(args.input_file) - input_modtran_wl = npzf['modtran_wavelengths'] + 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] @@ -270,29 +168,23 @@ def main(): 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 - - input_modtran = npzf['modtran_results'].copy() - input_solar_irr = npzf['solar_irradiance'] - points = npzf['points'] - - 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, input_modtran_wl[1] - input_modtran_wl[0]),dtype=np.float32) - output_modtran_irr = np.array(resample_spectrum(irr, irr_wl, output_modtran_wl, output_modtran_fwhm),dtype=np.float32) - - n_bands = int(input_modtran.shape[1] / len(npzf['keys'])) + 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'])): - 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')])) - - output_modtran = np.dot(input_modtran, H.T) + 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 = npzf['modtran_results'] + output_modtran = input_modtran output_modtran_irr = npzf['sol_irr'] output_modtran_wl = npzf['modtran_wavelengths'] @@ -313,127 +205,6 @@ def main(): if __name__ == "__main__": main() -#input_file = "/Users/brodrick/Downloads/emit_emulator_v1.hdf5" -##out_file = "/Users/brodrick/Downloads/emit_emulator_v1_conv01.hdf5" -#out_file = "/Users/brodrick/Downloads/emit_emulator_v1_conv_EMIT_Wavelengths_20220817_clipped.hdf5" -#out_file = "/Users/brodrick/repos/kf/convolved/neon_emulator_v1_conv.hdf5" -#subprocess.call(f"cp {input_file} {out_file}", shell=True) - -input_file = "data_202411121/emit_emulator_v2_0.hdf5" -out_file = "data_202411121/emit_emulator_v2_conv_EMIT_Wavelengths_20220817_clipped.hdf5" - -input = h5py.File(input_file, "r") -output = h5py.File(out_file, "w") - -def copy_groups(name, obj): - if isinstance(obj, h5py.Group): - output.create_group(name) -input.visititems(copy_groups) - -def copy_attributes(name, obj): - if isinstance(obj, h5py.Group) or isinstance(obj, h5py.Dataset): - for key, value in obj.attrs.items(): - output[name].attrs[key] = value -input.visititems(copy_attributes) - - -wl = input["MISCELLANEOUS"]["Wavelengths"][:] -X = input["sample_space"]["sample space"][:,:] -Y_all = input["mod_output/modtran output"][:,:,:] - - -output_wl = np.genfromtxt('wl/EMIT_Wavelengths_20220817.txt')[:,1:][::-1,:]*1000 -subset = np.where((output_wl[:,0] >= wl[0]) & (output_wl[:,0] <= wl[-1]))[0] -output_wl = output_wl[subset,:] -subset = np.where((output_wl[:,0] >= 380) & (output_wl[:,0] <= 2493))[0] -output_wl = output_wl[subset,:] -#output_wl = np.genfromtxt('wl/neon_wl.txt')[:,1:] - - -print('prepping') -for key, value in input["sample_space"].items(): - output.create_dataset(f"sample_space/{key}", data=value) -for key, value in input["MISCELLANEOUS"].items(): - output.create_dataset(f"MISCELLANEOUS/{key}", data=value) - -# overwrite wavelengths -del output["MISCELLANEOUS/Wavelengths"] -output["MISCELLANEOUS/Wavelengths"] = output_wl[:,0] - -yn = input['mod_output'].attrs['Products'].tolist() -del input -print('convolving') -##res = resample(wl, np.arange(wl[0], wl[-1], 0.1), 0.1*np.ones(len(np.arange(wl[0], wl[-1], 0.1)))) -##res = resample(wl, np.arange(wl[0], wl[-1], 5.0), 5.0*np.ones(len(np.arange(wl[0], wl[-1], 5.0)))) -print(output_wl[:,0].shape) -print(Y_all.shape) -print(wl.shape) -res = resample(wl, output_wl[:,0], output_wl[:,1]) - -yn_out = [] -Y_out = [] -yn_out.append('path_radiance') -Y_out.append(res.process_data(Y_all[:,yn.index('path_radiance'),:], Y_all.shape[0])) - -#yn_out.append('t_down_dir_t_up_dir') -yn_out.append('bi-direct') -Y_out.append(res.process_data(Y_all[:,yn.index('t_down_dir'),:]*\ - Y_all[:,yn.index('t_up_dir'),:]*\ - Y_all[:,yn.index('ToA_irrad'),:], Y_all.shape[0])) - -#yn_out.append('t_down_dif_t_up_dir') -yn_out.append('hemi-direct') -Y_out.append(res.process_data(Y_all[:,yn.index('t_down_dif'),:]*\ - Y_all[:,yn.index('t_up_dir'),:]*\ - Y_all[:,yn.index('ToA_irrad'),:], Y_all.shape[0])) - -#yn_out.append('t_down_dir_t_up_dif') -yn_out.append('direct-hemi') -Y_out.append(res.process_data(Y_all[:,yn.index('t_down_dir'),:]*\ - Y_all[:,yn.index('t_up_dif'),:]*\ - Y_all[:,yn.index('ToA_irrad'),:], Y_all.shape[0])) - -#yn_out.append('t_down_dif_t_up_dif') -yn_out.append('bi-hemi') -Y_out.append(res.process_data(Y_all[:,yn.index('t_down_dif'),:]*\ - Y_all[:,yn.index('t_up_dif'),:]*\ - Y_all[:,yn.index('ToA_irrad'),:], Y_all.shape[0])) - -yn_out.append('sphalb') -Y_out.append(res.process_data(Y_all[:,yn.index('sphalb_num'),:]/\ - Y_all[:,yn.index('sphalb_denom'),:], Y_all.shape[0])) - -Y_out = np.stack(Y_out, axis=1) - - -print('writing') -output["mod_output/modtran output"] = Y_out -output['mod_output'].attrs['Products'] = yn_out -del output - -print('checking') -input = h5py.File(out_file, "r") -Y_all = input["mod_output/modtran output"][:,:,:] -print(Y_all.shape) -del input - - -#Y_out = [] -#for n in range(Y_all.shape[1]): -# Y_out.append(res.process_data(Y_all[:,n,:], Y_all.shape[0])) -# -#Y_out = np.stack(Y_out, axis=1) -# -#del input["mod_output/modtran output"] -#input["mod_output/modtran output"] = Y_out -#del input["MISCELLANEOUS/Wavelengths"] -#input["MISCELLANEOUS/Wavelengths"] = output_wl[:,0] -#del input -#input = h5py.File(out_file, "r") -#Y_all = input["mod_output/modtran output"][:,:,:] -#print(Y_all.shape) - - diff --git a/train_combined_emulator.py b/train_combined_emulator.py index db3ed8c..ed4971d 100644 --- a/train_combined_emulator.py +++ b/train_combined_emulator.py @@ -250,6 +250,15 @@ def main(): 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: @@ -306,14 +315,15 @@ def main(): monitor='val_loss' es = keras.callbacks.EarlyStopping(monitor=monitor, mode='min', verbose=1, patience=20, restore_best_weights=True) - model = nn_model_ind(train_sixs.shape, train_modtran.shape, len(keys)) + #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=10, epochs=400, + 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