#!/usr/bin/env python # This was developped by Mathieu Compiegne at HYGEOS import sys sys.path.append("../f2py_utils") import numpy as np import pyartdeco_runlib as pyartdeco from f2py_utils import run_artdeco def artdeco_calib_glitter(dir_data, instrument, band_ind, aer_type, scene_var, ascii_kdis=True, verbose=False): """ -dir_data directory where to find input files -instrument instrument to simulate, it must correspond a the kdis filename -band_ind index of the band to model (BE CAREFUL it starts at 1) -aer_type either: none if no aerosols opac_arctic opac_cont_avg opac_cont_clean opac_cont_polluted opac_desert opac_maritime_clean opac_maritime_polluted opac_maritime_tropical opac_urban -scene_var dictionnary containing scene variable on which to loop e.g. scene_var = { "sza" : np.array([35.0]), degre "vza" : np.array([35.0]), degre "raa" : np.array([0.0]), relative azimuth in degre, 0 = forward scattering "windspeed" : np.array([5.0]), m/s "aot550" : np.array([0.0]), "watercolumn" : np.array([0.0]), g / cm2 "ozone" : np.array([0.0]), DU "chloro" : np.array([0.0]), pigment concentration (in mg/m^3) } -ascii_kdis either ascii or h5 Output: -reflectance [nscene_var, 3] : 3 because I,Q,U """ ########################################### # artdeco_in technical parameters structure if band_ind<1: print("") print(" ERROR ") print(" band_ind starts at 1") print("") return keywords = ['nowarn', 'od_no_check'] mode = 'kdis '+instrument filters = ['none'] wavel = [-1, -1] rt_model = 'doad' thermal = False corint = True nmat = 3 trunc_method = 'dm' # for the entire hemisphere nstreams_init = 8 artdeco_in = pyartdeco.artdeco_in(keywords, mode, filters, wavel, \ trunc_method, nstreams_init, rt_model, \ corint, thermal, nmat) #################### # read kdis coeff print("load kdis...") if ascii_kdis: kdis = pyartdeco.kdis_coeff(artdeco_in, dir_data+"kdis/"+instrument+"/", 'ascii', channel_list=np.array([band_ind])) else: kdis = pyartdeco.kdis_coeff(artdeco_in, dir_data+"kdis/"+instrument+"/", 'h5', channel_list=np.array([band_ind])) print("kdis loaded...") ###################### # load solar TOA flux solrad = pyartdeco.solar_irradiance(artdeco_in, file_path="none", file_format="ascii", channel_list=np.array([band_ind]), solar_spec="none") #print(solrad.__dict__) #################### # load atmosphere gas = ["o2", "co2", "o3", "h2o"] ppmv = [208493.3747, 400.3639, -1, -1] atm = 'atm_afgl_us62_red.dat' wldepol = np.array([0.1, 200.0]) # microns depol = np.array([0.0, 0.0]) atmos = pyartdeco.atmosphere(artdeco_in, dir_data+"atm/", atm, gas, ppmv, 'ascii', wldepol, depol, interp=True, kdis=kdis) ############################### # set ptcle structure wlref = 0.550 wlptcle = [np.min(kdis.wvlband), np.max(kdis.wvlband)] if aer_type == "none": ptcle_opt_path=[] else: ptcle_opt_path=[{"name":aer_type, "file_path":dir_data+"opt/opt_opac.h5"}] print("read aer optical properties") ptcle_opt = pyartdeco.ptcle_optical_properties(ptcle_opt_path, artdeco_in.nstreams, wlptcle, wlref, read_betal=True, opt_interp=True) print("read aer optical properties ok") # scale heigh of OPAC aerosols z_aer={"opac_antarctic":8.0, "opac_arctic":99.0, "opac_cont_avg":8.0, "opac_cont_clean":8.0, "opac_cont_polluted":8.0, "opac_desert":2.0, "opac_maritime_clean":1.0, "opac_maritime_polluted":1.0, "opac_maritime_tropical":1.0, "opac_urban":8.0, } if aer_type == "none": ptcle_def = [] else: ptcle_def = [ {"name":aer_type, "tau":np.nan, "alt_distrib":"sh", "z":z_aer[aer_type], "humidity":80.0} ] if len(ptcle_def) > 1: print("") print(" ERROR ") print(" Only one ptcle maximum admitted") print("") return ptcles = pyartdeco.particle(artdeco_in, ptcle_def, ptcle_opt) # print(ptcles.__dict__) # after that, there is no more I/O # start loop over the scene_var nscenes = len(scene_var["sza"]) refl = np.full( (nscenes,3), np.nan ) for isc in range(nscenes): ######################## # Set glitter surface = pyartdeco.surface("surface", "brdf", "ocean", ocean_whitecaps = True, shadow=False, xsal=34.3, pcl=scene_var["chloro"][isc], wdspd=scene_var["windspeed"][isc] ) ######################## # Geometry sza = np.array([ scene_var["sza"][isc] ]) vza = np.array([ scene_var["vza"][isc] ]) raa = np.array([ scene_var["raa"][isc] ]) geom = pyartdeco.geometry(sza, vza, raa) ######################## # Set aer AOT if len(ptcle_def) == 1: ptcles.tau_reflamb[0] = scene_var["aot550"][isc] ######################## # set ozone and water vapor atmos.apply_scale_O3(scene_var["ozone"][isc]) atmos.apply_scale_h2o(scene_var["watercolumn"][isc]) ############### # run ARTDECO lamb, rad, rad_levels, flux, alt = run_artdeco(artdeco_in, atmos, surface, solrad, ptcles, geom, kdis=kdis, verbose=verbose) # flux = np.zeros((artdeco.mcommon.nsza, 3, artdeco.mcommon.nalt_atm, artdeco.mcommon.nlambda)) # alt_flux = np.zeros(artdeco.mcommon.nalt_atm) # lamb = np.zeros(artdeco.mcommon.nlambda) # radiance = np.zeros((artdeco.mcommon.nsza, artdeco.mcommon.nvza, artdeco.mcommon.nvaa, artdeco.mcommon.nmat, artdeco.mcommon.nlambda)) # radiance_levels = np.zeros((artdeco.mcommon.nsza, artdeco.mcommon.nvza, artdeco.mcommon.nvaa, artdeco.mcommon.nmat, artdeco.mcommon.nalt_atm, artdeco.mcommon.nlambda)) if verbose: print("fluxes") print(" alt total down up direct down") for ialt in range(len(alt)): print(" %12.6f "%alt[ialt]+" %16.9f"%flux[0,0,ialt,0] +" %16.9f"%flux[0,1,ialt,0] +" %16.9f"%flux[0,2,ialt,0]) print('radiances') print(lamb[0], vza[0], rad[0,0,0,0:3,0]) refl[isc,:] = np.pi * rad[0,0,0,0:3,0] / np.cos(sza[0] * np.pi /180.0) return lamb[0], refl if __name__=='__main__': print("Library file")