Skip to content
Snippets Groups Projects
Commit 0c8882bf authored by tomrink's avatar tomrink
Browse files

snapshot...

parent 8ae8c5a8
No related branches found
No related tags found
No related merge requests found
import xarray as xr
import numpy as np
from util.geos_nav import GEOSNavigation
from netCDF4 import Dataset
import pickle
from util.util import haversine_np
def romio_to_fgf(filename, goes_e_w='EAST'):
ds = xr.open_dataset(filename, engine='cfgrib')
lons = ds['longitude']
lats = ds['latitude']
lons = lons.values
lats = lats.values
lat_array, lon_array = np.meshgrid(lats, lons, indexing='ij')
nav = GEOSNavigation()
if goes_e_w == 'EAST':
nav = GEOSNavigation()
elif goes_e_w == 'WEST':
nav = GEOSNavigation(sub_lon=-137.0)
cc, ll = nav.earth_to_lc_s(lon_array.flatten(), lat_array.flatten())
tlons, tlats = nav.lc_to_earth(cc, ll)
dist = haversine_np(lon_array.flatten(), lat_array.flatten(), tlons, tlats)
ok = np.invert(np.isnan(dist))
print(np.average(dist[ok]))
print(np.histogram(dist[ok]))
cc = cc.reshape(lon_array.shape)
ll = ll.reshape(lon_array.shape)
f = open('/home/rink/elems.pkl', 'wb')
pickle.dump(cc, f)
f.close()
f = open('/home/rink/lines.pkl', 'wb')
pickle.dump(ll, f)
f.close()
ds.close()
return cc, ll
def create_file(filename, fgf_elem, fgf_line):
dim_x_len = fgf_elem.shape[1]
dim_y_len = fgf_elem.shape[0]
rootgrp = Dataset(filename, 'w', format='NETCDF4')
dim_x = rootgrp.createDimension('romio_x', size=dim_x_len)
dim_y = rootgrp.createDimension('romio_y', size=dim_y_len)
romio_fgf_elems = rootgrp.createVariable('romio_fgf_elems', 'f4', ['romio_y', 'romio_x'])
romio_fgf_lines = rootgrp.createVariable('romio_fgf_lines', 'f4', ['romio_y', 'romio_x'])
romio_fgf_elems[:, :] = fgf_elem[:, :]
romio_fgf_lines[:, :] = fgf_line[:, :]
rootgrp.close()
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()
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment