import os
import numpy
# 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, logical_and
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
wind_err_thresh = 50
check_wind_err = True
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')
# new field for Mie scattering ratio
# sr0 = coda.fetch(product, 'mie_wind_prod_conf_data', -1, 'mie_wind_qc', 'fitting_mie_sr')
#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
if check_wind_err:
totlevs = numpy.sum(logical_and(mie_valid[i, :] > 0, logical_and(rid[i, :] > 0, mie_err[i, :] < wind_err_thresh)))
else:
totlevs = numpy.sum(logical_and(mie_valid[i, :] > 0, rid[i, :] > 0))
if totlevs > 0:
if typ_id[i] == 1: # cloudy type Mie winds
print(yyyy, mm, dd, hh, min, sec, '%7.2f %7.2f %2i' % (float(lonrid[i]), float(latrid[i]), int(totlevs)), file=f60)
if check_wind_err:
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('%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)
else:
for m in range(nlevm):
if rid[i, m] > 0:
if mie_valid[i, m] > 0.0:
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
if check_wind_err:
totlevs = numpy.sum(logical_and(ray_valid[i, :] > 0, logical_and(rid[i, :] > 0, ray_err[i, :] < wind_err_thresh)))
else:
totlevs = numpy.sum(logical_and(ray_valid[i, :] > 0, rid[i, :] > 0))
# for clear sky Rayleigh winds ----------
if totlevs > 0:
if typ_id[i] == 2:
print(yyyy, mm, dd, hh, min, sec, '%7.2f %7.2f %2i' % (float(lonrid[i]), float(latrid[i]), int(totlevs)), file=f66)
if check_wind_err:
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('%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)
else:
for m in range(nlev):
if rid[i, m] > 0:
if ray_valid[i, m] > 0.0:
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()