diff --git a/modules/util/sss1day_FMB_py3.py b/modules/util/sss1day_FMB_py3.py new file mode 100644 index 0000000000000000000000000000000000000000..fa1525520f6b7c34a93a5f7e56cf02da901f21a5 --- /dev/null +++ b/modules/util/sss1day_FMB_py3.py @@ -0,0 +1,314 @@ +import os + +# Import CODA definitions + # Older dates +# os.putenv('CODA_DEFINITION', '/home/huiliu/CODA/share/coda/definitions/') +# os.putenv('CODA_DEFINITION', '/home/stevew/CODA/share/coda/definitions/') +os.putenv('CODA_DEFINITION', '/data/Personal/stevew/AEOLUS/CODA/share/coda/definitions/') + # Newer dates +#os.putenv('CODA_DEFINITION', '/home/huiliu/CODA/share/coda/definitions/AEOLUS-20190611.codadef') + +from numpy import vstack, zeros, linspace, where + +import sys +# Tell python scripts where to find CODA (next import coda) -------- +# sys.path.append('/home/huiliu/CODA/lib/python2.7/site-packages') +# sys.path.append('/home/stevew/CODA/lib/python2.7/site-packages') +sys.path.append('/data/Personal/stevew/AEOLUS/CODA/lib/python3.7/site-packages') +import coda + +import matplotlib.pyplot as plt +import glob, ast, datetime + +files = open('flist_adm.txt').read().split() +print(files) + +nfile = len(files) +allobs = 0 +ntot = 0 +ntotm = 0 +wind_err_thresh = 50 + +f66=open('./ray1day.out', 'w+') +f60=open('./mie1day.out', 'w+') + +# ------- loop over data files ------ + +for n, filename in enumerate(files): + print('Reading file: %s' %filename) + + product = coda.open(filename) + +### ------- Mie wind profile --------- +### ------- Mie wind profile --------- + +# print("Individual Mie HLOS wind points") +# latitude = coda.fetch(product, 'mie_geolocation',-1, 'windresult_geolocation/latitude_cog') +# longitude = coda.fetch(product, 'mie_geolocation',-1, 'windresult_geolocation/longitude_cog') +# mie_alt in (m) + mie_alt0 = coda.fetch(product, 'mie_geolocation',-1, 'windresult_geolocation/altitude_vcog') + mie_altt = coda.fetch(product, 'mie_geolocation',-1, 'windresult_geolocation/altitude_top') + mie_altb = coda.fetch(product, 'mie_geolocation',-1, 'windresult_geolocation/altitude_bottom') + mie_azimuth0 = coda.fetch(product, 'mie_geolocation',-1, 'windresult_geolocation/los_azimuth') + + mie_length0 = coda.fetch(product, 'mie_hloswind',-1, 'windresult/integration_length') + mie_valid0 = coda.fetch(product, 'mie_hloswind',-1, 'windresult/validity_flag') + mie_wind0 = coda.fetch(product, 'mie_hloswind',-1, 'windresult/mie_wind_velocity') + + #klukens +# mie_pppp0 = coda.fetch(product, 'mie_hloswind',-1, 'windresult/reference_pressure') +# mie_pppp0 = coda.get_field_names(product,'mie_geolocation[0]/windresult_geolocation') +# mie_pppp0 = coda.get_field_names(product,'mie_hloswind[0]/windresult') +# print "mie_pppp0 fetch all = ",mie_pppp0 + + mie_err0 = coda.fetch(product, 'mie_wind_prod_conf_data',-1, 'mie_wind_qc/hlos_error_estimate') +# mie_snr0 = coda.fetch(product, 'mie_wind_prod_conf_data',-1, 'mie_wind_qc/mie_snr') + +# ------- Rayleigh profiles --------- +# ------- Rayleigh profiles --------- + +# print("Individual Rayleight HLOS wind points") + rayleigh_azimuth0 = coda.fetch(product, 'rayleigh_geolocation',-1, 'windresult_geolocation/los_azimuth') + +# latitude = coda.fetch(product, 'rayleigh_geolocation',-1, 'windresult_geolocation/latitude_cog') +# longitude = coda.fetch(product, 'rayleigh_geolocation',-1, 'windresult_geolocation/longitude_cog') + + rayleigh_alt0 = coda.fetch(product, 'rayleigh_geolocation',-1, 'windresult_geolocation/altitude_vcog') + rayleigh_altt = coda.fetch(product, 'rayleigh_geolocation',-1, 'windresult_geolocation/altitude_top') + rayleigh_altb = coda.fetch(product, 'rayleigh_geolocation',-1, 'windresult_geolocation/altitude_bottom') + + rayleigh_wind0 = coda.fetch(product, 'rayleigh_hloswind',-1, 'windresult/rayleigh_wind_velocity') + +# ----- rayleigh data is from top (26.5km) --> bottom (24 levels) --------- +# wind_err in (m/s) + + rayleigh_err0 = coda.fetch(product, 'rayleigh_wind_prod_conf_data',-1, 'rayleigh_wind_qc/hlos_error_estimate') + rayleigh_sratio0 = coda.fetch(product, 'rayleigh_wind_prod_conf_data',-1, 'rayleigh_wind_qc/scattering_ratio') + rayleigh_wind_to_T = coda.fetch(product, 'rayleigh_hloswind',-1, 'windresult/rayleigh_wind_to_temperature') + rayleigh_wind_to_P = coda.fetch(product, 'rayleigh_hloswind',-1, 'windresult/rayleigh_wind_to_pressure') + + rayleigh_temp = coda.fetch(product, 'rayleigh_hloswind',-1, 'windresult/reference_temperature') + rayleigh_pppp = coda.fetch(product, 'rayleigh_hloswind',-1, 'windresult/reference_pressure') + + ray_length0 = coda.fetch(product, 'rayleigh_hloswind',-1, 'windresult/integration_length') + rayleigh_valid0 = coda.fetch(product, 'rayleigh_hloswind',-1, 'windresult/validity_flag') +# print(rayleigh_wind0.shape) + + +#================================================= +## ---------- Mie profile information ------------ +#================================================= + print("Mie HLOS wind profiles") + + rid = coda.fetch(product, 'mie_profile',-1, 'l2b_wind_profiles/wind_result_id_number') + rid = vstack(rid) +# print(rid.shape) + + typ_id = coda.fetch(product, 'mie_profile',-1, 'l2b_wind_profiles/Obs_Type') + + nprofm = rid.shape[0] + nlevm = rid.shape[1] + +# ----------------- the data is from top --> bottom ------------- +# NOTE: i = 0 - nprof-1 in the xrange FUNTION of PYTHON +#----------------------------------------------------------------- + + mie_azimuth = zeros(rid.shape) + mie_err = zeros(rid.shape) + mie_hhh = zeros(rid.shape) + mie_hht = zeros(rid.shape) + mie_hhb = zeros(rid.shape) + mie_wind = zeros(rid.shape) +# mie_pppp = zeros(rid.shape) #klukens + mie_valid = zeros(rid.shape) + mie_length = zeros(rid.shape) + mie_sratio = zeros(rid.shape) +# mie_snr = zeros(rid.shape) + +# wind m/s, + + mie_wind[rid !=0] = mie_wind0[rid[rid !=0]-1]*0.01 + mie_azimuth[rid !=0] = mie_azimuth0[rid[rid !=0]-1]*1.0 + + mie_valid[rid !=0] = mie_valid0[rid[rid !=0]-1] +# length in km output + mie_length[rid !=0] = mie_length0[rid[rid !=0]-1]*0.001 + +# mie_pppp[rid !=0] = mie_pppp0[rid[rid !=0]-1]*0.01 #klukens + +# wind error m/s + mie_err[rid !=0] = mie_err0[rid[rid !=0]-1]*1.0 +# mie_snr [rid !=0] = mie_snr0 [rid[rid !=0]-1]*1.0 + +# height in (km) +# add 250m shift to the height for this version of L2B data + mie_hhh[rid !=0] = mie_alt0[rid[rid !=0]-1] *0.001+0.25 + + mie_hht[rid !=0] = mie_altt[rid[rid !=0]-1] *0.001 +0.25 + mie_hhb[rid !=0] = mie_altb[rid[rid !=0]-1] *0.001 +0.25 + + latrid = coda.fetch(product, 'mie_profile',-1, 'Profile_lat_average') + lonrid = coda.fetch(product, 'mie_profile',-1, 'Profile_lon_average') + sstime = coda.fetch(product, 'mie_profile',-1, 'Profile_DateTime_Average') + + for i in range(nprofm): + timestep = datetime.datetime(2000, 1, 1)+datetime.timedelta(seconds=sstime[i]) + yyyy = timestep.strftime('%Y') + mm = timestep.strftime('%m') + dd = timestep.strftime('%d') + hh = timestep.strftime('%H') + min = timestep.strftime('%M') + sec = timestep.strftime('%S') #klukens + + totlevs = 0 + for m in range(nlevm): + if rid[i, m] > 0: + if mie_err[i,m] < wind_err_thresh: + if mie_valid[i,m] > 0.0: + totlevs = totlevs + 1 + + ntotm = ntotm + 1 + +# profile number from 0 --> nprofm +# ------------------------------------------------------ +# f60=open('./test/mie.dat.'+ str(ntotm).zfill(5), 'w+') +# ------------------------------------------------------ + if totlevs > 0: + if typ_id[i] == 1: +# cloudy type Mie winds + + # print >>f60, yyyy, mm, dd, hh, min, sec, '%7.2f %7.2f %2i' %(float(lonrid[i]), float(latrid[i]), int(totlevs)) + print(yyyy, mm, dd, hh, min, sec, '%7.2f %7.2f %2i' %(float(lonrid[i]), float(latrid[i]), int(totlevs)), file=f60) + + for m in range(nlevm): +# keep consistent with the lines of totlevs above !!! + if rid[i, m] > 0: + if mie_err[i,m] < wind_err_thresh: + if mie_valid[i,m] > 0.0: + # print >>f60, '%2i %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f' %(int(m+1), \ + # mie_hhh[i,m], mie_hht[i,m], mie_hhb[i,m], mie_err[i,m], \ + # mie_azimuth[i,m], mie_wind[i,m], mie_length[i,m]) + print('%2i %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f' %(int(m+1), + mie_hhh[i,m], mie_hht[i,m], mie_hhb[i,m], mie_err[i,m], + mie_azimuth[i,m], mie_wind[i,m], mie_length[i,m]), file=f60) + + +# ----------------------------------------------- +# Rayleigh profile information +# ----------------------------------------------- + + print("Rayleigh HLOS wind profiles") + rid = coda.fetch(product, 'rayleigh_profile',-1, 'l2b_wind_profiles/wind_result_id_number') + rid = vstack(rid) + +#========================================== +# tpye =1, cloudy, =2 clear +#========================================== + typ_id = coda.fetch(product, 'rayleigh_profile',-1, 'l2b_wind_profiles/Obs_Type') + +# profile number of Rayleigh winds profiles in this orbit + nprof = rid.shape[0] + nlev = rid.shape[1] +# print(rid.shape) + +# ----------------- the data is from top --> bottom ------------- +# NOTE: i = 0 - nprof-1 in the xrange FUNTION of PYTHON +#----------------------------------------------------------------- + + ray_azimuth = zeros(rid.shape) + ray_err = zeros(rid.shape) + wind_sens_T = zeros(rid.shape) + wind_sens_P = zeros(rid.shape) + ray_length = zeros(rid.shape) + ray_valid = zeros(rid.shape) + ref_temp = zeros(rid.shape) + ref_pppp = zeros(rid.shape) + ray_hhh = zeros(rid.shape) + ray_hht = zeros(rid.shape) + ray_hhb = zeros(rid.shape) + ray_wind = zeros(rid.shape) + ray_sratio = zeros(rid.shape) + + + ray_wind[rid !=0] = rayleigh_wind0[rid[rid !=0]-1]*0.01 + ray_azimuth[rid !=0] = rayleigh_azimuth0[rid[rid !=0]-1]*1.0 + ray_sratio [rid !=0] = rayleigh_sratio0[rid[rid !=0]-1]*1.0 + +# print(rayleigh_wind0.shape) +# print(rid.shape) +# print(rid) +# from 2D OBS fortran index rid (1-nnn), to index 0-(nnn-1) of 1D mie_wind0 +# + +# wind error m/s + ray_err[rid !=0] = rayleigh_err0[rid[rid !=0]-1]*1.0 + +# m/s/K + wind_sens_T[rid !=0] = rayleigh_wind_to_T[rid[rid !=0]-1]*0.01 +# K + ref_temp[rid !=0] = rayleigh_temp[rid[rid !=0]-1] + ref_pppp[rid !=0] = rayleigh_pppp[rid[rid !=0]-1]*0.01 + +# cm/s/hPa + wind_sens_P[rid !=0] = rayleigh_wind_to_P[rid[rid !=0]-1] *1.0e-4 + ray_length[rid !=0] = ray_length0[rid[rid !=0]-1] *0.001 + ray_valid[rid !=0] = rayleigh_valid0[rid[rid !=0]-1] + +# height in (km) +# add 250m shift to the height for this version of L2B data + ray_hhh[rid !=0] = rayleigh_alt0[rid[rid !=0]-1] *0.001 +0.25 + ray_hht[rid !=0] = rayleigh_altt[rid[rid !=0]-1] *0.001 +0.25 + ray_hhb[rid !=0] = rayleigh_altb[rid[rid !=0]-1] *0.001 +0.25 + + latrid =coda.fetch(product, 'rayleigh_profile',-1, 'Profile_lat_average') + lonrid =coda.fetch(product, 'rayleigh_profile',-1, 'Profile_lon_average') + sstime =coda.fetch(product, 'rayleigh_profile',-1, 'Profile_DateTime_Average') + + for i in range(nprof): + timestep = datetime.datetime(2000, 1, 1) + datetime.timedelta(seconds=sstime[i]) + yyyy = timestep.strftime('%Y') + mm = timestep.strftime('%m') + dd = timestep.strftime('%d') + hh = timestep.strftime('%H') + min = timestep.strftime('%M') + sec = timestep.strftime('%S') #klukens + + ntot = ntot + 1 +# print >>f4, str(rid[i, :]) + + totlevs = 0 + for m in range(nlev): + if rid[i, m] > 0: + if ray_err[i,m] < wind_err_thresh: + if ray_valid[i, m] > 0.0: + totlevs = totlevs + 1 + +# profile number from 0 --> nprof +# f66=open('./test/ray.dat.'+ str(ntot).zfill(5), 'w+') +# for clear sky Rayleigh winds ---------- + if totlevs > 0: + if typ_id[i] == 2: + + # print >>f66, yyyy, mm, dd, hh, min, sec, '%7.2f %7.2f %2i' %(float(lonrid[i]), float(latrid[i]), int(totlevs)) + print(yyyy, mm, dd, hh, min, sec, '%7.2f %7.2f %2i' %(float(lonrid[i]), float(latrid[i]), int(totlevs)), file=f66) + + for m in range(nlev): +# keep consistent with the lines of totlevs above !!! + if rid[i, m] > 0: + if ray_err[i,m] < wind_err_thresh: + if ray_valid[i, m] > 0.0: + # print >>f66, '%2i %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f' \ + # %(int(m+1), ray_hhh[i,m], ray_hht[i,m], ray_hhb[i,m],\ + # ray_err[i,m], ray_azimuth[i,m], ray_wind[i,m],\ + # ref_temp[i,m]*0.01, ref_pppp[i,m], wind_sens_T[i,m],\ + # wind_sens_P[i,m], ray_sratio [i,m], ray_length[i,m]) + print('%2i %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f' + %(int(m+1), ray_hhh[i,m], ray_hht[i,m], ray_hhb[i,m], + ray_err[i,m], ray_azimuth[i,m], ray_wind[i,m], + ref_temp[i,m]*0.01, ref_pppp[i,m], wind_sens_T[i,m], + wind_sens_P[i,m], ray_sratio [i,m], ray_length[i,m]), file=f66) + + coda.close(product) + +print('Finished data. Processed, file loop') +exit()