#!/usr/bin/env python from __future__ import print_function, division, absolute_import from builtins import range import sys import os import string import numpy as np from scipy.interpolate import interp1d import h5py from luts import Idx, read_mlut_hdf5 import matplotlib ## Force matplotlib to not use any Xwindows backend. #matplotlib.use('Agg') import matplotlib.pyplot as plt fsize = 12 matplotlib.rcParams.update({'font.size': fsize}) size_limit_loop = 1e8 ######################################################################################################### def skipcomment(f): pos = f.tell() tmp = f.readline() while tmp.startswith('#'): pos = f.tell() tmp = f.readline() f.seek(pos) def find_nearest(array, value): array = np.asarray(array) idx = (np.abs(array - value)).argmin() return idx def get_solrad(sol_spec): dir_pykdis = "/home/mathieu/work/RTM/pykdis/" if sol_spec=="kurudz92": f = open(dir_pykdis+"data/solrad_kurudz92_full.dat") skipcomment(f) tmp = f.readline() F0 = float(tmp) skipcomment(f) tmp = f.readline() nFlamb = int(tmp) Flamb = np.zeros(nFlamb) lamb = np.zeros(nFlamb) skipcomment(f) for i in range(nFlamb): tmp = f.readline() lamb[i] = float(tmp.split()[0]) # microns Flamb[i] = float(tmp.split()[1]) # W m-2 microns-1 f.close() sorted_lamb = np.argsort(lamb) n = 40 # print(lamb.shape) # print(lamb.shape[0]%n) cut = lamb.shape[0]%n lamb = lamb[sorted_lamb][:-cut] Flamb = Flamb[sorted_lamb][:-cut] # print(lamb.shape) # print(lamb.shape[0]%n) Flamb = np.average(Flamb.reshape(-1, n), axis=1) lamb = np.average(lamb.reshape(-1, n), axis=1) # print(lamb.shape) sol_const = {"wn":1e4/lamb, "lamb":lamb, "F0":F0, "Flamb":Flamb} elif sol_spec=="meftah2018": tmp_sol = np.genfromtxt(dir_pykdis+"data/solrad_meftah.dat", comments='#', names=['wvl', 'F'], usecols=(0,1)) lamb = tmp_sol['wvl'] Flamb = tmp_sol['F'] lamb = lamb[np.isfinite(Flamb)] Flamb = Flamb[np.isfinite(Flamb)] sorted_lamb = np.argsort(lamb) F0 = 1372.3 # see Meftah, M.; Damé, L.; Bolsée, D.; Hauchecorne, A.; Pereira, N.; Sluse, D.; Cessateur, G.; Irbah, A.; Bureau, J.; Weber, M.; Bramstedt, K.; Hilbig, T.; Thiéblemont, R.; Marchand, M.; Lefèvre, F.; Sarkissian, A. & Bekki, S. SOLAR-ISS: A new reference spectrum based on SOLAR/SOLSPEC observations A&A, 2018, 611, A1 sol_const = {"wn":1e7/lamb[sorted_lamb], "lamb":lamb[sorted_lamb]*1e-3, "F0":F0, "Flamb":Flamb[sorted_lamb]*1e3} elif sol_spec=="thuillier2003": tmp_sol = np.genfromtxt(dir_pykdis+"data/solrad_thuillier2003.dat", comments='#', names=['wvl', 'F'], usecols=(0,1)) lamb = tmp_sol['wvl'] Flamb = tmp_sol['F'] lamb = lamb[np.isfinite(Flamb)] Flamb = Flamb[np.isfinite(Flamb)] res_sol = 0.1 #nm #print(np.max(lamb), np.min(lamb)) num = int((np.max(lamb)-np.min(lamb)) / res_sol) newlamb = np.linspace(np.min(lamb), np.max(lamb), num=num) Flamb = interp1d(lamb,Flamb)(newlamb) sorted_lamb = np.argsort(newlamb) sol_const = {"wn":1e7/newlamb[sorted_lamb], "lamb":newlamb[sorted_lamb]*1e-3, "F0":-1, "Flamb":Flamb[sorted_lamb]*10} sol_const['Fwn'] = sol_const['Flamb'] * sol_const['lamb'] / sol_const['wn'] # add F solar in W m-2 cm return sol_const ###################################################################################################### def doit(lut_spectral_file, res_file, kdis_name, kdis_dir, wvl_axis_name="wavelengths", use_solspec=True, excluded_dataset=[]): if os.path.isfile(res_file): print("") print("") print("Res file already present") print(" "+res_file) print("") print("") return kdis_file = kdis_dir+"kdis_"+kdis_name+".h5" if not os.path.isfile(kdis_file): print("") print("") print(" Missing kdis file") print(" "+kdis_file) print("") print("") return fkdis = h5py.File(kdis_file, "r") wvl = np.copy(fkdis["def"]["central_wvl"]) fwhm = fkdis["def"]["max_wvl"][:] - fkdis["def"]["min_wvl"][:] if "isrf" in list(fkdis.keys()): isrfs = [] ind_isrf = np.array([],dtype=int) isrf_list = fkdis['isrf'].attrs.get("name_list").decode().split(",") for iwvl, name in enumerate(isrf_list): #print(wvl[iwvl], name) isrfs.append( {"wvl_c":np.copy(fkdis['isrf'][name]["wvl_c"])*1e-3, "trans":np.copy(fkdis['isrf'][name]["trans"]), "wvl":np.copy(fkdis['isrf'][name]["wvl"])*1e-3} ) ind_isrf = np.append(ind_isrf, iwvl) else: ind_isrf = np.full_like(wvl, -1, dtype=int) if use_solspec: solspec = fkdis['solrad']['solrad'].attrs.get("solspec").decode() fkdis.close() if use_solspec: solrad = get_solrad(solspec) if not 0.55 in wvl: wvl = np.append(0.55, wvl) fwhm = np.append(0., fwhm) ind_isrf = np.append(-1, ind_isrf) if not 0.65 in wvl: wvl = np.append(0.65, wvl) fwhm = np.append(0., fwhm) ind_isrf = np.append(-1, ind_isrf) if not 0.865 in wvl: wvl = np.append(0.865, wvl) fwhm = np.append(0., fwhm) ind_isrf = np.append(-1, ind_isrf) nwvl = len(wvl) fwhm = fwhm[np.argsort(wvl)] ind_isrf = ind_isrf[np.argsort(wvl)] wvl = np.sort(wvl) is_isrf = np.full_like(ind_isrf, False, dtype=bool) is_isrf[ind_isrf != -1] = True hh = h5py.File(lut_spectral_file, "r") groups = list(hh.keys()) hh.close() for group in groups: m_spectral = read_mlut_hdf5(lut_spectral_file, group=group, lazy=True) if not wvl_axis_name in m_spectral.axes.keys(): print("") print("") print(" No -"+wvl_axis_name+"- axis on LUT ") print(" Following axes are present: " ) print(" "+m_spectral.axes.keys() ) print("") print("") return f = h5py.File(res_file, 'a') grp = f.require_group(group) print("Create axis in res file...") subgrp = grp.require_group("axis") dset = subgrp.create_dataset(wvl_axis_name, data=wvl, dtype='float64') dset = subgrp.create_dataset("FWHM", data=fwhm, dtype='float64') dset = subgrp.create_dataset("is_isrf", data=is_isrf, dtype=bool) for ax in m_spectral.axes.keys(): if ax != wvl_axis_name: dset = subgrp.create_dataset(ax, data=m_spectral.axes[ax]) subgrp = grp.require_group("data") print("Create dataset in res file...") hh = h5py.File(lut_spectral_file, "r") data_list = [] for i, (name, dataset, axes, attrs) in enumerate(m_spectral.data): if not name in excluded_dataset: data_list.append(name) #print("") #print(name) dim = () attr_ax = "" for ia, ax in enumerate(axes): #print(ia,ax) if ax == wvl_axis_name: dim = dim + (len(wvl),) else: dim = dim + (len(m_spectral.axes[ax]),) attr_ax = attr_ax + ax + "," #print(dim) attr_ax= attr_ax[:-1] #print(dim) #print(attr_ax) dset = subgrp.require_dataset(name, dim, dtype='float64') #dset.attrs.create("dimensions", np.string_(attr_ax)) for attrsname, attrsvalue in hh[group]['data'][name].attrs.items(): dset.attrs.create(attrsname, attrsvalue) hh.close() f.close() print("Compute dataset ...") for name in data_list: if not name in excluded_dataset: print(" ", name, group, kdis_name) del m_spectral m_spectral = read_mlut_hdf5(lut_spectral_file, group=group, datasets=[name] ) #m_spectral.print_info() if not wvl_axis_name in m_spectral.axes.keys(): f = h5py.File(res_file, 'a') grp = f.require_group(group) subgrp = grp.require_group("data") subgrp[name][...] = m_spectral[name].data[...] f.close() continue print(" %d bytes" % (m_spectral[name].data.size * m_spectral[name].data.itemsize)) loop_iax = 0 loop_ax_name = m_spectral[name].names[loop_iax] if loop_ax_name == wvl_axis_name: loop = False else: loop = True if (m_spectral[name].data.size * m_spectral[name].data.itemsize) < size_limit_loop: loop = False if loop: ntmp = m_spectral[name].data.shape[loop_iax] else: ntmp= 1 if ((m_spectral[name].data.size * m_spectral[name].data.itemsize) / ntmp) > size_limit_loop: print("") print("ERROR") print("You will be in trouble with memory usage") print("") return if loop: print(" Loop") for iwvl in range(len(wvl)): #print("\n treat wvl=%.3f"%wvl[iwvl], group) if is_isrf[iwvl]: if np.round_(isrfs[ind_isrf[iwvl]]['wvl_c'], decimals=8) != np.round_(wvl[iwvl], decimals=8): print("") print("") print(" Pb with the ISRF") print("") print(isrfs[ind_isrf[iwvl]]['wvl_c'], wvl[iwvl]) print("") exit() if use_solspec and (wvl[iwvl]<4.0): iwvl_min = find_nearest(solrad['lamb'], np.min(isrfs[ind_isrf[iwvl]]['wvl']) ) if iwvl_min!= 0: iwvl_min = iwvl_min-1 iwvl_max = find_nearest(solrad['lamb'], np.max(isrfs[ind_isrf[iwvl]]['wvl']) ) if iwvl_max != len(solrad['lamb']-1) : iwvl_max = iwvl_max+1 solspec = solrad['Flamb'][ iwvl_min:iwvl_max+1 ] wvl_solspec = solrad['lamb'][ iwvl_min:iwvl_max+1 ] # plt.figure() # plt.plot(isrfs[ind_isrf[iwvl]]['wvl'], isrfs[ind_isrf[iwvl]]['trans']) # plt.plot(wvl_solspec, solspec/np.max(solspec)) # plt.savefig(kdis_name+"_%i"%(iwvl+1)+".png") fisrf = interp1d(isrfs[ind_isrf[iwvl]]['wvl'], isrfs[ind_isrf[iwvl]]['trans'], bounds_error = False, fill_value = 0.0) srf = solspec * fisrf(wvl_solspec) wvl_srf = wvl_solspec else: srf = isrfs[ind_isrf[iwvl]]['trans'] wvl_srf = isrfs[ind_isrf[iwvl]]['wvl'] elif fwhm[iwvl] > 0.0: wvl_srf = np.linspace(wvl[iwvl]-fwhm[iwvl]/2.0, wvl[iwvl]+fwhm[iwvl]/2.0, num=100) if use_solspec and (wvl[iwvl]<4.0): iwvl_min = find_nearest(solrad['lamb'], wvl[iwvl]-fwhm[iwvl]/2.0 ) if iwvl_min!= 0: iwvl_min = iwvl_min-1 iwvl_max = find_nearest(solrad['lamb'], wvl[iwvl]+fwhm[iwvl]/2.0 ) if iwvl_max != len(solrad['lamb']-1) : iwvl_max = iwvl_max+1 solspec = solrad['Flamb'][ iwvl_min:iwvl_max+1 ] wvl_solspec = solrad['lamb'][ iwvl_min:iwvl_max+1 ] fsolspec = interp1d(wvl_solspec,solspec) srf = fsolspec(wvl_srf) else: srf = np.ones_like(srf) else: wvl_srf = wvl[iwvl] if is_isrf[iwvl] or fwhm[iwvl] > 0.0: norm_srf = np.trapz(srf, x=wvl_srf) if loop: # print(loop_iax) # print(loop_ax_name) for ii in range(m_spectral[name].data.shape[loop_iax]): sub_lut = m_spectral.sub({loop_ax_name:ii}).sub({wvl_axis_name:Idx(wvl_srf)}) # sub_lut = m_spectral.sub( {loop_ax_name:ii, wvl_axis_name:Idx(wvl_srf)} ) if is_isrf[iwvl] or fwhm[iwvl] > 0.0: iaxis_wvl_tmp = 0 for ia, ax in enumerate(sub_lut[name].names): if ax == wvl_axis_name: iaxis_wvl_tmp = ia res_tmp = np.trapz(sub_lut[name].data * srf, x=wvl_srf, axis=iaxis_wvl_tmp) / norm_srf else: res_tmp = sub_lut[name].data[...] # print(res_tmp.shape) # print(indice) indice = () for ia, ax in enumerate(m_spectral[name].names): #print(ia,ax) if ax == wvl_axis_name: indice = indice + (iwvl,) elif ax == loop_ax_name: indice = indice + (ii,) else: indice = indice + (slice(len(m_spectral.axes[ax])),) f = h5py.File(res_file, 'a') grp = f.require_group(group) subgrp = grp.require_group("data") subgrp[name][indice] = res_tmp f.close() del sub_lut else: indice = () iaxis_wvl = 0 for ia, ax in enumerate(m_spectral[name].names): #print(ia,ax) if ax == wvl_axis_name: iaxis_wvl = ia indice = indice + (iwvl,) else: indice = indice + (slice(len(m_spectral.axes[ax])),) sub_lut = m_spectral.sub({wvl_axis_name:Idx(wvl_srf)}) if is_isrf[iwvl] or fwhm[iwvl] > 0.0: res_tmp = np.trapz(sub_lut[name].data * srf, x=wvl_srf, axis=iaxis_wvl) / norm_srf else: res_tmp = sub_lut[name].data[...] f = h5py.File(res_file, 'a') grp = f.require_group(group) subgrp = grp.require_group("data") subgrp[name][indice] = res_tmp f.close() del sub_lut if __name__=='__main__': # kdis_name = "ckdmip_s3a_olci_+1p5nm" # kdis_dir = "/rfs/proj/pykdis/output/"+kdis_name+"/" # lut_spectral_file = "/rfs/data/mopsmap/aer_opt_lut/ahmad2010_aer_spectral_enhanced_fromrhindiv.h5" # file_res = "/rfs/data/mopsmap/aer_opt_lut/opt_ahmad2010_aer_enhanced_fromrhindiv_"+kdis_name+".h5" # doit(lut_spectral_file, file_res, kdis_name, kdis_dir, excluded_dataset=["normed_ext_coeff"]) # kdis_name = "ckdmip_s3a_olci_-1p5nm" # kdis_dir = "/rfs/proj/pykdis/output/"+kdis_name+"/" # lut_spectral_file = "/rfs/data/mopsmap/aer_opt_lut/ahmad2010_aer_spectral_enhanced_fromrhindiv.h5" # file_res = "/rfs/data/mopsmap/aer_opt_lut/opt_ahmad2010_aer_enhanced_fromrhindiv_"+kdis_name+".h5" # doit(lut_spectral_file, file_res, kdis_name, kdis_dir, excluded_dataset=["normed_ext_coeff"]) # kdis_name = "ckdmip_map" # kdis_dir = "/rfs/proj/map-clim-sds/input/artdeco_static/"+kdis_name+"/" # ptcles = "baum_ice_asc" # lut_spectral_file = "/rfs/proj/map-clim-sds/input/artdeco_static/opt_"+ptcles+".h5" # file_res = "/scratch/map-clim-sds/opt_"+ptcles+"_"+kdis_name+".h5" # doit(lut_spectral_file, file_res, kdis_name, kdis_dir, excluded_dataset=["normed_ext_coeff"]) # ptcles = "baum_ice_ghm" # lut_spectral_file = "/rfs/proj/map-clim-sds/input/artdeco_static/opt_"+ptcles+".h5" # file_res = "/scratch/map-clim-sds/opt_"+ptcles+"_"+kdis_name+".h5" # doit(lut_spectral_file, file_res, kdis_name, kdis_dir, excluded_dataset=["normed_ext_coeff"]) # ptcles = "libradtran_liquid" # lut_spectral_file = "/rfs/proj/map-clim-sds/input/artdeco_static/opt_"+ptcles+".h5" # file_res = "/scratch/map-clim-sds/opt_"+ptcles+"_"+kdis_name+".h5" # doit(lut_spectral_file, file_res, kdis_name, kdis_dir, excluded_dataset=["normed_ext_coeff"]) # ptcles = "cams_aer" # lut_spectral_file = "/rfs/proj/map-clim-sds/input/artdeco_static/"+ptcles+".h5" # file_res = "/scratch/map-clim-sds/opt_"+ptcles+"_"+kdis_name+".h5" # doit(lut_spectral_file, file_res, kdis_name, kdis_dir, excluded_dataset=["normed_ext_coeff"]) # kdis_name = "ckdmip_msg3_seviri_visswir" # kdis_dir = "/rfs/proj/pykdis/output/"+kdis_name+"/" # ptcles = "baum_ice_asc" # lut_spectral_file = "/rfs/proj/geosds/input/artdeco_static/opt_"+ptcles+".h5" # file_res = "/scratch/geosds/opt_"+ptcles+"_msg3_seviri_visswir.h5" # doit(lut_spectral_file, file_res, kdis_name, kdis_dir, excluded_dataset=["normed_ext_coeff"]) # ptcles = "baum_ice_ghm" # lut_spectral_file = "/rfs/proj/geosds/input/artdeco_static/opt_"+ptcles+".h5" # file_res = "/scratch/geosds/opt_"+ptcles+"_msg3_seviri_visswir.h5" # doit(lut_spectral_file, file_res, kdis_name, kdis_dir, excluded_dataset=["normed_ext_coeff"]) # ptcles = "libradtran_liquid" # lut_spectral_file = "/rfs/proj/geosds/input/artdeco_static/opt_"+ptcles+".h5" # file_res = "/scratch/geosds/opt_"+ptcles+"_msg3_seviri_visswir.h5" # doit(lut_spectral_file, file_res, kdis_name, kdis_dir, excluded_dataset=["normed_ext_coeff"]) # ptcles = "cams_aer" # lut_spectral_file = "/rfs/proj/geosds/input/artdeco_static/"+ptcles+".h5" # file_res = "/scratch/geosds/opt_"+ptcles+"_msg3_seviri_visswir.h5" # doit(lut_spectral_file, file_res, kdis_name, kdis_dir, excluded_dataset=["normed_ext_coeff"]) # kdis_name = "ckdmip_FCI" # # kdis_dir = "/rfs/proj/pykdis/output/"+kdis_name+"/" # kdis_dir = "/rfs/proj/geosds/input/artdeco_static/"+kdis_name+"/" # ptcles = "baum_ice_asc" # lut_spectral_file = "/rfs/proj/geosds/input/artdeco_static/opt_"+ptcles+".h5" # file_res = "/scratch/geosds/opt_"+ptcles+"_fci.h5" # doit(lut_spectral_file, file_res, kdis_name, kdis_dir, excluded_dataset=["normed_ext_coeff"]) # ptcles = "baum_ice_ghm" # lut_spectral_file = "/rfs/proj/geosds/input/artdeco_static/opt_"+ptcles+".h5" # file_res = "/scratch/geosds/opt_"+ptcles+"_fci.h5" # doit(lut_spectral_file, file_res, kdis_name, kdis_dir, excluded_dataset=["normed_ext_coeff"]) # ptcles = "libradtran_liquid" # lut_spectral_file = "/rfs/proj/geosds/input/artdeco_static/opt_"+ptcles+".h5" # file_res = "/scratch/geosds/opt_"+ptcles+"_fci.h5" # doit(lut_spectral_file, file_res, kdis_name, kdis_dir, excluded_dataset=["normed_ext_coeff"]) # ptcles = "cams_aer" # lut_spectral_file = "/rfs/proj/geosds/input/artdeco_static/"+ptcles+".h5" # file_res = "/scratch/geosds/opt_"+ptcles+"_fci.h5" # doit(lut_spectral_file, file_res, kdis_name, kdis_dir, excluded_dataset=["normed_ext_coeff"]) # kdis_name = "ckdmip_METimage" # kdis_dir = "/rfs/proj/geosds/input/artdeco_static/"+kdis_name+"/" # ptcles = "baum_ice_asc" # lut_spectral_file = "/rfs/proj/geosds/input/artdeco_static/opt_"+ptcles+".h5" # file_res = "/scratch/geosds/opt_"+ptcles+"_metimage.h5" # doit(lut_spectral_file, file_res, kdis_name, kdis_dir, excluded_dataset=["normed_ext_coeff"]) # ptcles = "baum_ice_ghm" # lut_spectral_file = "/rfs/proj/geosds/input/artdeco_static/opt_"+ptcles+".h5" # file_res = "/scratch/geosds/opt_"+ptcles+"_metimage.h5" # doit(lut_spectral_file, file_res, kdis_name, kdis_dir, excluded_dataset=["normed_ext_coeff"]) # ptcles = "libradtran_liquid" # lut_spectral_file = "/rfs/proj/geosds/input/artdeco_static/opt_"+ptcles+".h5" # file_res = "/scratch/geosds/opt_"+ptcles+"_metimage.h5" # doit(lut_spectral_file, file_res, kdis_name, kdis_dir, excluded_dataset=["normed_ext_coeff"]) # ptcles = "cams_aer" # lut_spectral_file = "/rfs/proj/geosds/input/artdeco_static/"+ptcles+".h5" # file_res = "/scratch/geosds/opt_"+ptcles+"_metimage.h5" # doit(lut_spectral_file, file_res, kdis_name, kdis_dir, excluded_dataset=["normed_ext_coeff"]) # kdis_name = "ckdmip_msg3_seviri" # kdis_dir = "/scratch/pykdis/output/"+kdis_name+"/" # ptcles = "opac" # lut_spectral_file = "/rfs/proj/artdeco_lib/opt/opt_"+ptcles+".h5" # file_res = "/scratch/geosds/opt_"+ptcles+"_msg3_seviri.h5" # doit(lut_spectral_file, file_res, kdis_name, kdis_dir, excluded_dataset=["normed_ext_coeff"]) # kdis_name = "ckdmip_msg3_seviri" # kdis_dir = "/scratch/pykdis/output/"+kdis_name+"/" # ptcles = "baum_ice_ghm" # lut_spectral_file = "/rfs/proj/artdeco_lib/opt/opt_"+ptcles+".h5" # file_res = "/scratch/geosds/opt_"+ptcles+"_msg3_seviri.h5" # doit(lut_spectral_file, file_res, kdis_name, kdis_dir, excluded_dataset=["normed_ext_coeff"]) # kdis_name = "ckdmip_FCI" # kdis_dir = "/scratch/pykdis/output/"+kdis_name+"/" # ptcles = "opac" # lut_spectral_file = "/rfs/proj/artdeco_lib/opt/opt_"+ptcles+".h5" # file_res = "/scratch/geosds/opt_"+ptcles+"_fci.h5" # doit(lut_spectral_file, file_res, kdis_name, kdis_dir, excluded_dataset=["normed_ext_coeff"])