diff --git a/Readme.rtf b/Readme.rtf new file mode 100644 index 0000000..6945173 --- /dev/null +++ b/Readme.rtf @@ -0,0 +1,17 @@ +{\rtf1\ansi\ansicpg1252\cocoartf2761 +\cocoatextscaling0\cocoaplatform0{\fonttbl\f0\froman\fcharset0 Times-Bold;\f1\froman\fcharset0 Times-Roman;} +{\colortbl;\red255\green255\blue255;\red0\green0\blue0;} +{\*\expandedcolortbl;;\cssrgb\c0\c0\c0;} +\margl1440\margr1440\vieww25460\viewh17460\viewkind0 +\deftab720 +\pard\pardeftab720\sa280\partightenfactor0 + +\f0\b\fs36 \cf0 \expnd0\expndtw0\kerning0 +Steps to Generate the Lookup Table and Create Plots: +\f1\b0 \ +\pard\tx220\tx720\pardeftab720\li720\fi-720\partightenfactor0 +\cf0 1. Run `modify_json.py` to generate the JSON files for various combinations of environmental conditions.\ +2. Run `run_n_modtran.py` to execute Modtran using those JSON files.\ +3. Run `generate_LUT_and_plot.py` to generate the lookup table (stored in all_r.npy) and create plots. The lookup table has five dimensions: sensor zenith angle, ground altitude, water vapor content, and solar zenith angle.\ +\ +} \ No newline at end of file diff --git a/all_r.npy b/all_r.npy new file mode 100644 index 0000000..bb1a78b Binary files /dev/null and b/all_r.npy differ diff --git a/center.npy b/center.npy new file mode 100644 index 0000000..a02f4a0 Binary files /dev/null and b/center.npy differ diff --git a/generate_LUT_and_plot.py b/generate_LUT_and_plot.py new file mode 100644 index 0000000..2241d44 --- /dev/null +++ b/generate_LUT_and_plot.py @@ -0,0 +1,205 @@ +# Five degrees lookup table +import numpy as np +import os +import read_chn_files +import scipy.ndimage +from matplotlib import pyplot as plt +import matplotlib +matplotlib.use('TkAgg') # Or 'MacOSX' + +invertflag=1 # 1: convert zenith angle from 120-180 to 0-60 degrees +useexist=1 # 1: use existing LUT; 0: create new LUT from .chn files +actualpath=0 # 1: calculate the actual path length that light travels through the atmosphere, given the solar zenith angle and the sensor zenith angle + +def get_5deg_lookup_index(zenith, ground, water, solarz,conc): + idx = np.asarray([[zenith],[ground],[water],[solarz],[conc]]) + # idx = np.asarray([zenith, ground, water, solarz, conc]) + return idx + +def spline_5deg_lookup(grid_data, zenith, ground, water, solarz, conc): + order=1 + coords = get_5deg_lookup_index( + zenith=zenith, ground=ground, water=water, solarz=solarz, conc=conc) + coords_fractional_part, coords_whole_part = np.modf(coords) + # print(coords_whole_part) + coords_near_slice = tuple((slice(int(cc.item()), int(cc.item()+2)) for cc in coords_whole_part)) + near_grid_data = grid_data[coords_near_slice] + new_coord = np.concatenate((coords_fractional_part * np.ones((1, near_grid_data.shape[-1])), + np.arange(near_grid_data.shape[-1])[None, :]), axis=0) + if order == 1: + lookup = scipy.ndimage.map_coordinates(near_grid_data, coordinates=new_coord, order=1, mode='nearest') + elif order == 3: + lookup = np.asarray([scipy.ndimage.map_coordinates( + im, coordinates=coords_fractional_part, order=order, mode='nearest') for im in np.moveaxis(near_grid_data, 5, 0)]) + return lookup.squeeze() + +def generate_library(gas_concentration_vals, zenith, ground, water, solarz, grid): + rads = np.empty((len(gas_concentration_vals), grid.shape[-1])) + for i, ppmm in enumerate(gas_concentration_vals): + rads[i, :] = spline_5deg_lookup( + grid, zenith=zenith, ground=ground, water=water, solarz=solarz, conc=ppmm) + return rads + +# Generate LUT using .chn files +def create_4d_matrix(num): + #path = '/Users/xiang/Desktop/data/c21' + current_folder = os.getcwd() + path = os.path.join(current_folder, 'chn_files') + all_r=np.zeros((num, num, num, num, num,285)) + + z_ind=0 + for z in np.linspace(120,180,num): + z_ind=z_ind+1 + g_ind=0 + for g in np.linspace(0,3,num): + g_ind=g_ind+1 + w_ind=0 + for w in np.linspace(0,6,num): + w_ind=w_ind+1 + s_ind=0 + for s in np.linspace(0,60,num): + s_ind=s_ind+1 + c_ind=0 + for c in np.linspace(0,120,num): + c_ind=c_ind+1 + + base_name=f'z{z_ind}_g{g_ind}_w{w_ind}_s{s_ind}_c{c_ind}' + fn=base_name+'.chn' + filename= os.path.join(path, fn) + + if os.path.isfile(filename): # radiance is already resampled in .chn files + t, wl, r, center, fwhm = read_chn_files.load_chn(filename, 1) + #logr = np.log(r, out=np.zeros_like(wl), where=r > 0) + all_r[z_ind-1][g_ind-1][w_ind-1][s_ind-1][c_ind-1] = r + else: + print('Missed file') + return all_r, center + +def find_ind(num_table,value_min,value_max,value): # convert value of parameter to its index in LUT + ind=(value-value_min)/(value_max-value_min)*(num_table-1) + return ind + + +num=11 # number of lines (parameter values) in the plot +num_table=5 # number of values for each parameter in the LUT + +if (os.path.exists('all_r.npy')) & useexist == 1: # If LUT exists, load the data from the file + all_r = np.load('all_r.npy') # LUT + center = np.load('center.npy') # corresponding wavelength +else: # If LUT doesn't exist, generate LUT and save it to the file + print('Generating new lookup table') + all_r, center = create_4d_matrix(num_table) + np.save('all_r.npy', all_r) + np.save('center.npy', center) + +smaller_range=1 +maxz, minz, maxg, ming, maxw, minw, maxs, mins = 180, 120, 3, 0, 6, 0, 60, 0 # range of sensor zenith angle, ground altitude, water vapor and solar zenith angle +SCALING = 1e5 +concentrations_ind=range(1,4,1) +c_list=np.linspace(0,120,num_table) # Attendiont: make it consistant with modify_json.py used to generate LUT +concentrations = c_list[concentrations_ind] +#cl = [c * 0.5 * 1000 for c in concentrations] # convert concentration to concentration length + +# Generate plots for each of the four parameters—zenith, ground altitude, water content, and solar zenith angle. +# Each plot shows how the spectrum varies with one specific parameter, while keeping the remaining parameters fixed at their median value. +for para in [0,1,2,3]: + plt.figure(dpi=150) + + if para==0: # sensor zenith angle + g=1.5 + w=3 + s=30 + g_ind=find_ind(num_table,ming,maxg, g) + w_ind=find_ind(num_table,minw,maxw, w) + s_ind = find_ind(num_table, mins, maxs, s) + para_list=np.linspace(minz,maxz,num) + + if para==1: # ground altitude + z=150 + w=3 + s=30 + z_ind = find_ind(num_table, minz, maxz, z) + w_ind=find_ind(num_table,minw, maxw, w) + s_ind = find_ind(num_table, mins, maxs, s) + para_list=np.linspace(ming, maxg, num) + + if para==2: # w + z=150 + g=1.5 + s=30 + z_ind = find_ind(num_table, minz, maxz, z) + g_ind=find_ind(num_table,ming,maxg, g) + s_ind=find_ind(num_table,mins,maxs, s) + para_list=np.linspace(minw, maxw, num) + + if para==3: # sensor zenith angle + z=150 + g=1.5 + w=3 + z_ind = find_ind(num_table, minz, maxz, z) + g_ind=find_ind(num_table,ming,maxg,g) + w_ind=find_ind(num_table,minw,maxw,w) + para_list=np.linspace(mins,maxs,num) + + spectra_matrix = np.empty((len(para_list), len(center))) + j=0 + + if (para == 0) & (invertflag == 1): + para_list = para_list[::-1] + para_list2 = 180 - para_list + else: + para_list2 = para_list + + + for i in para_list: + if para==0: + z_ind = find_ind(num_table, minz, maxz, i) + z=i + if para==1: + g_ind=find_ind(num_table,ming,maxg,i) + g=i + if para==2: + w_ind=find_ind(num_table,minw,maxw,i) + w=i + if para == 3: + s_ind=find_ind(num_table,mins,maxs,i) + s=i + + rads= generate_library(concentrations_ind, z_ind, g_ind, w_ind, s_ind, all_r) + logr = np.log(rads, out=np.zeros_like(rads), where=rads > 0) + + layer_thickness=0.5*1000 + if actualpath==1: # calculate the actual path length that light travels through the atmosphere, given the solar zenith angle and the sensor zenith angle + path_length=layer_thickness/np.cos(s*np.pi/180)+layer_thickness/np.cos((180-z)*np.pi/180) + else: + path_length=layer_thickness + cl = [c * path_length for c in concentrations] # Multiply each element in the list + + slope, _, _, _ = np.linalg.lstsq(np.stack((np.ones_like(cl), cl)).T, logr, rcond=None) + spectrum = slope[1, :] * SCALING + spectra_matrix[j, :] = spectrum + j=j+1 + + cmap = plt.get_cmap('Blues') + minc=min(para_list2) + minc=minc-0.3 #increase the shade of the lightest line + maxc=max(para_list2) + color = cmap((para_list2[j-1]-minc)/ (maxc-minc)) + plt.plot(center, spectrum, label=f'Line {j}', c=color) + plt.xlabel('Wavelength') + plt.ylabel(r'$dlog(L)/dl_{\mathrm{CH4}}$') + if smaller_range==1: # plot methane absorption region + plt.xlim(2100,2500) + + degree_symbol = "\u00B0" + names = [f"Sensor zenith angle ({degree_symbol})", "Ground altitude (km)", "Water vapor (g)", + f"Solar zenith angle ({degree_symbol})"] + sm = plt.cm.ScalarMappable(cmap=cmap, norm=plt.Normalize(vmin=minc, vmax=maxc)) + sm.set_array([]) + cbar = plt.colorbar(sm, ax=plt.gca(), pad=0.1) + cbar.set_label(f'{names[para]}') + # plt.ylim(-0.1,1e-5) + # plt.title('New Lookup Table') + + plt.tight_layout() + plt.show() \ No newline at end of file diff --git a/modify_json.py b/modify_json.py new file mode 100644 index 0000000..bab0df2 --- /dev/null +++ b/modify_json.py @@ -0,0 +1,72 @@ +import json +import numpy +import os + +def np_encoder(object): + if isinstance(object, numpy.generic): + return object.item() + +output_folder = "json_files" +if not os.path.exists(output_folder): + os.makedirs(output_folder) + +num=5 # the number of values per parameter +z_ind=0 +for z in numpy.linspace(120,180,num): + z_ind=z_ind+1 + g_ind=0 + for g in numpy.linspace(0,3,num): + g_ind=g_ind+1 + w_ind=0 + for w in numpy.linspace(0,6,num): + w_ind=w_ind+1 + c_ind=0 + for c in numpy.linspace(0,120,num): + c_ind=c_ind+1 + s_ind=0 + for s in numpy.linspace(0,60,num): + s_ind=s_ind+1 + + with open("phil_original.json", "r") as jsonFile: + data = json.load(jsonFile) + data["MODTRAN"][0]["MODTRANINPUT"]["GEOMETRY"]["H1ALT"] = 400 + + data["MODTRAN"][0]["MODTRANINPUT"]["GEOMETRY"]["OBSZEN"] = z + data["MODTRAN"][0]["MODTRANINPUT"]["SURFACE"]["GNDALT"] = g + data["MODTRAN"][0]["MODTRANINPUT"]["ATMOSPHERE"]["H2OSTR"] = w + data["MODTRAN"][0]["MODTRANINPUT"]["GEOMETRY"]["PARM2"] = s + + dist=data["MODTRAN"][0]["MODTRANINPUT"]["ATMOSPHERE"]["PROFILES"][0]["PROFILE"] + dist[0]=g + dist[1]=g+0.5 + dist[2]=g+0.5001 # methane is concentrated in the bottom layer + + new_elements = [] + new_inds = [] + + for index, element in enumerate(dist[3:]): + if element > dist[2]: + new_elements.append(element) + new_inds.append(index + 3) + + new_dist = dist[:3] + new_elements + data["MODTRAN"][0]["MODTRANINPUT"]["ATMOSPHERE"]["PROFILES"][0]["PROFILE"]=new_dist + data["MODTRAN"][0]["MODTRANINPUT"]["ATMOSPHERE"]["NLAYERS"]=len(new_dist) + + conc=data["MODTRAN"][0]["MODTRANINPUT"]["ATMOSPHERE"]["PROFILES"][1]["PROFILE"] + indices_of_new_dist= list(range(3))+new_inds + conc[0:2]=[c,c] # methane is concentrated in the bottom layer + new_conc = [conc[index] for index in indices_of_new_dist] + data["MODTRAN"][0]["MODTRANINPUT"]["ATMOSPHERE"]["PROFILES"][1]["PROFILE"]=new_conc + + + base_name = 'z' + str(z_ind) + '_g' + str(g_ind) + '_w' + str(w_ind) + '_s' + str(s_ind) + '_c' + str(c_ind) + fn=base_name+'.json' + filename = os.path.join(output_folder, fn) + data["MODTRAN"][0]["MODTRANINPUT"]["NAME"]=base_name + + with open(filename, "w") as jsonFile: + json.dump(data, jsonFile, default=np_encoder,indent=4) + + + diff --git a/phil_original.json b/phil_original.json new file mode 100644 index 0000000..0590ec7 --- /dev/null +++ b/phil_original.json @@ -0,0 +1,189 @@ +{ + "MODTRAN": [ + { + "MODTRANINPUT": { + "AEROSOLS": { + "IHAZE": "AER_RURAL", + "VIS": -0.01 + }, + "ATMOSPHERE": { + "CO2MX": 410.0, + "H2OOPT": "+", + "H2OSTR": 2.8, + "H2OUNIT": "g", + "HMODEL": "CH4_PL", + "M1": "ATM_MIDLAT_SUMMER", + "M2": "ATM_MIDLAT_SUMMER", + "M3": "ATM_MIDLAT_SUMMER", + "M4": "ATM_MIDLAT_SUMMER", + "M5": "ATM_MIDLAT_SUMMER", + "M6": "ATM_MIDLAT_SUMMER", + "MODEL": "ATM_USER_ALT_PROFILE", + "NLAYERS": 49, + "NPROF": 2, + "O3STR": 0.3, + "O3UNIT": "a", + "PROFILES": [ + { + "PROFILE": [ + 1, + 1.5, + 3.501, + 4.0, + 5.0, + 6.0, + 7.0, + 8.0, + 9.0, + 10.0, + 11.0, + 12.0, + 13.0, + 14.0, + 15.0, + 16.0, + 17.0, + 18.0, + 19.0, + 20.0, + 21.0, + 22.0, + 23.0, + 24.0, + 25.0, + 27.5, + 30.0, + 32.5, + 35.0, + 37.5, + 40.0, + 42.5, + 45.0, + 47.5, + 50.0, + 55.0, + 60.0, + 65.0, + 70.0, + 75.0, + 80.0, + 85.0, + 90.0, + 95.0, + 100.0, + 105.0, + 110.0, + 115.0, + 120.0 + ], + "TYPE": "PROF_ALTITUDE", + "UNITS": "UNT_KILOMETERS" + }, + { + "PROFILE": [ + 257.85, + 257.85, + 1.85, + 1.84674, + 1.83585, + 1.81953, + 1.7945, + 1.77274, + 1.7575, + 1.71832, + 1.67806, + 1.64106, + 1.6095, + 1.57903, + 1.54747, + 1.51265, + 1.47565, + 1.43974, + 1.39403, + 1.332, + 1.25582, + 1.16006, + 1.05885, + 0.957647, + 0.8584, + 0.766771, + 0.687221, + 0.608541, + 0.544988, + 0.484591, + 0.426153, + 0.368803, + 0.31265, + 0.259435, + 0.211553, + 0.171288, + 0.163235, + 0.163235, + 0.163235, + 0.163235, + 0.163235, + 0.163235, + 0.152353, + 0.141471, + 0.130588, + 0.119706, + 0.103382, + 0.0652941, + 0.0326471 + ], + "TYPE": "PROF_CH4", + "UNITS": "UNT_DPPMV" + } + ] + }, + "CASE": 0, + "DESCRIPTION": "", + "FILEOPTIONS": { + "CKPRNT": true, + "NOPRNT": 2 + }, + "GEOMETRY": { + "GMTIME": 0, + "H1ALT": 400, + "IDAY": 200, + "IPARM": 12, + "ITYPE": 3, + "OBSZEN": 150, + "PARM1": 0, + "PARM2": 60, + "TRUEAZ": 0 + }, + "NAME": "IND_7881_solzen-60_solaz-0_GNDALT-1_H2OSTR-2.8_AERFRAC2-0.01_ch4enh-128_senzen-150_senaz-0_altitude-400", + "RTOPTIONS": { + "DISALB": true, + "IEMSCT": "RT_SOLAR_AND_THERMAL", + "IMULT": "RT_DISORT", + "LYMOLC": false, + "MODTRN": "RT_CORRK_FAST", + "NSTR": 8, + "SOLCON": 0.0, + "T_BEST": false + }, + "SPECTRAL": { + "BMNAME": "p1_2013", + "DV": 0.1, + "FILTNM": "/beegfs/scratch/brodrick/emit/reflectance_analyses/modtran_comp/emit20220830t052228/lut_full/wavelengths_modtran_381.0055927_2492.9238809.flt", + "FLAGS": "NT A ", + "FWHM": 0.1, + "V1": 340.0, + "V2": 2520.0, + "XFLAG": "N", + "YFLAG": "R" + }, + "SURFACE": { + "GNDALT": 1, + "NSURF": 1, + "SURFP": { + "CSALB": "LAMB_CONST_5_PCT" + }, + "SURFTYPE": "REFL_LAMBER_MODEL" + } + } + } + ] +} \ No newline at end of file diff --git a/read_chn_files.py b/read_chn_files.py new file mode 100644 index 0000000..a17095c --- /dev/null +++ b/read_chn_files.py @@ -0,0 +1,114 @@ +from sys import platform +import os +import logging +from os.path import split +import re +import json +from copy import deepcopy +import scipy as s +from scipy.stats import norm as normal +from scipy.interpolate import interp1d +import numpy as np +import pdb + +def load_chn(infile, coszen, return_specrad = False): + """Load a '.chn' output file and parse critical coefficient vectors. + + These are: + wl - wavelength vector + sol_irr - solar irradiance + sphalb - spherical sky albedo at surface + transm - diffuse and direct irradiance along the + sun-ground-sensor path + transup - transmission along the ground-sensor path only + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + Be careful with these! They are to be used only by the modtran_tir functions because + MODTRAN must be run with a reflectivity of 1 for them to be used in the RTM defined + in radiative_transfer.py. + thermal_upwelling - atmospheric path radiance + thermal_downwelling - sky-integrated thermal path radiance reflected off the ground + and back into the sensor. + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + We parse them one wavelength at a time.""" + with open(infile) as f: + sols, transms, sphalbs, wls, rhoatms, transups, spec_rads = [], [], [], [], [], [], [] + rdnatms = [] + centers, fwhms=[], [] + thermal_upwellings, thermal_downwellings = [], [] + lines = f.readlines() + for i, line in enumerate(lines): + if i < 5: + continue + toks = line.strip().split(' ') + toks = re.findall(r"[\S]+", line.strip()) + wl, wid = float(toks[0]), float(toks[8]) # nm + solar_irr = float(toks[18]) * 1e6 * \ + np.pi / wid / coszen # uW/nm/sr/cm2 + rdnatm = float(toks[4]) * 1e6 # uW/nm/sr/cm2 + rhoatm = rdnatm * np.pi / (solar_irr * coszen) + #sphalb = float(toks[23]) + sphalb = 1 + transm = float(toks[22]) + float(toks[21]) + transup = float(toks[24]) + + spec_rad = float(toks[4]) + + # Be careful with these! See note in function comments above + thermal_emission = float(toks[11]) + thermal_scatter = float(toks[12]) + thermal_upwelling = (thermal_emission + thermal_scatter) / wid * 1e6 # uW/nm/sr/cm2 + + # Be careful with these! See note in function comments above + grnd_rflt = float(toks[16]) + thermal_downwelling = grnd_rflt / wid * 1e6 # uW/nm/sr/cm2 + center=float(toks[-5]) + fwhm=float(toks[-2]) + + + sols.append(solar_irr) + transms.append(transm) + sphalbs.append(sphalb) + rhoatms.append(rhoatm) + transups.append(transup) + rdnatms.append(rdnatm) + centers.append(center) + fwhms.append(fwhm) + + spec_rads.append(spec_rad) + + thermal_upwellings.append(thermal_upwelling) + thermal_downwellings.append(thermal_downwelling) + + wls.append(wl) + + params = [np.array(i) for i in [wls, sols, rhoatms, transms, sphalbs, transups, thermal_upwellings, thermal_downwellings, spec_rads]] + + return tuple(params), np.array(wls),np.array(rdnatms), np.array(centers), np.array(fwhms) + +def change_temperature_profile(input_modtran_json_filename, output_modtran_json_filename, deltaT_K, tropopause_km = 17.): + with open(input_modtran_json_filename) as f: + js_in = json.load(f) + + nprof = len(js_in['MODTRAN'][0]['MODTRANINPUT']['ATMOSPHERE']['PROFILES']) + for i in range(nprof): + if js_in['MODTRAN'][0]['MODTRANINPUT']['ATMOSPHERE']['PROFILES'][i]['TYPE'] == 'PROF_TEMPERATURE': + prof_temperature = np.array(js_in['MODTRAN'][0]['MODTRANINPUT']['ATMOSPHERE']['PROFILES'][i]['PROFILE']) + + if js_in['MODTRAN'][0]['MODTRANINPUT']['ATMOSPHERE']['PROFILES'][i]['TYPE'] == 'PROF_ALTITUDE': + prof_altitude = np.array(js_in['MODTRAN'][0]['MODTRANINPUT']['ATMOSPHERE']['PROFILES'][i]['PROFILE']) + + new_prof_temperature = np.where(prof_altitude <= tropopause_km, prof_temperature + deltaT_K, prof_temperature) + + for i in range(nprof): + if js_in['MODTRAN'][0]['MODTRANINPUT']['ATMOSPHERE']['PROFILES'][i]['TYPE'] == 'PROF_TEMPERATURE': + js_in['MODTRAN'][0]['MODTRANINPUT']['ATMOSPHERE']['PROFILES'][i]['PROFILE'] = new_prof_temperature.tolist() + + + with open(output_modtran_json_filename, 'w') as f: + json.dump(js_in, f, sort_keys = True, indent = 4) + + + + \ No newline at end of file diff --git a/read_tp7.py b/read_tp7.py new file mode 100644 index 0000000..04e07ac --- /dev/null +++ b/read_tp7.py @@ -0,0 +1,39 @@ +from sys import platform +import os +import logging +from os.path import split +import re +import json +from copy import deepcopy +import scipy as s +from scipy.stats import norm as normal +from scipy.interpolate import interp1d +import numpy as np +import pdb + +def load_data(infile, return_specrad = False): + + with open(infile) as f: + rads, freqs = [], [] + + lines = f.readlines() + for i, line in enumerate(lines): + if i < 11: + continue + toks = line.strip().split(' ') + + # print(len(toks)) + if len(toks) >= 10: + + toks = re.findall(r"[\S]+", line.strip()) + freq, rad = float(toks[0]), float(toks[9]) # nm + + freqs.append(freq) + rads.append(rad) + + return s.array(rads),s.array(freqs) + + + + + \ No newline at end of file diff --git a/run_n_modtran.py b/run_n_modtran.py new file mode 100644 index 0000000..6e13e25 --- /dev/null +++ b/run_n_modtran.py @@ -0,0 +1,98 @@ + + + + +import subprocess +import glob +import os, sys +import time +import random +import argparse +import multiprocessing +import json +import numpy as np + +parser = argparse.ArgumentParser() +parser.add_argument('lut_dir',type=str) +parser.add_argument('-n_cpu',type=int, default=-1) +parser.add_argument('-n_jobs',type=int, default=-1) +parser.add_argument('-modtran_exe',type=str,default='/beegfs/store/shared/MODTRAN6/MODTRAN6.0.0/bin/linux/mod6c_cons') +args = parser.parse_args() + + + + +def run_single_modtran(modtran_exe,json_file,check_files,lock_file): + completed = True + for _f in check_files: + if (os.path.isfile(_f) is False): + completed = False + if completed: + print('already completed: {}'.format(json_file)) + return + + underway = os.path.isfile(lock_file) + if underway: + print('already underway: {}'.format(lock_file)) + return + + print('starting: '+json_file) + + cmd_str = 'touch {}'.format(lock_file) + subprocess.call(cmd_str,shell=True) + time.sleep(random.randint(1, 10) / 10.) + + cmd_str = '{} {}'.format(modtran_exe, json_file) + subprocess.call(cmd_str,shell=True) + cmd_str = 'rm {}'.format(lock_file) + subprocess.call(cmd_str,shell=True) + for name in ['r_k', + "t_k", + "wrn", + "psc", + "plt", + "7sc", + "acd"]: + subprocess.call('rm {}.{}'.format(os.path.splitext(check_files[0])[0], name), shell=True) + + + return + +subprocess.call('rm -f /dev/shm/*',shell=True) +#subprocess.call('cp run_one_modtran.csh {}'.format(args.lut_dir),shell=True) +os.chdir(args.lut_dir) +subprocess.call('mkdir logs',shell=True) + +run_files = glob.glob('./*.json', recursive=True) + +lock_files = [os.path.splitext(x)[0].replace('LUT_','') + '.lock' for x in run_files] +tp6_files = [os.path.splitext(x)[0].replace('LUT_','') + '.tp6' for x in run_files] +chn_files = [os.path.splitext(x)[0].replace('LUT_','') + '.chn' for x in run_files] + +if (args.n_cpu == -1): + args.n_cpu = multiprocessing.cpu_count() + +pool = multiprocessing.Pool(processes=args.n_cpu) +results = [] +if args.n_jobs == -1: + max_runs = 100000 +else: + max_runs = args.n_jobs +run=0 +for _f in range(len(run_files)): + if np.all([os.path.isfile(cf) is False for cf in [lock_files[_f], tp6_files[_f], chn_files[_f]]]): + results.append(pool.apply_async(run_single_modtran, args=(args.modtran_exe,run_files[_f],[tp6_files[_f],chn_files[_f]],lock_files[_f],))) + run+=1 + + if run > max_runs: + break + + +#results = [p.get() for p in results] +pool.close() +pool.join() + + + + +