Skip to content
Snippets Groups Projects
Commit 1cce267f authored by rink's avatar rink
Browse files

initial commit

parent 3b2327c8
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
modules/util/sss1day_FMB_py3.py 0 → 100644
+ 314
− 0
View file @ 1cce267f
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()
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment