diff --git a/modules/util/romio.py b/modules/util/romio.py deleted file mode 100644 index ca6b58ca62c5b5de87c588b740422e32df4b1191..0000000000000000000000000000000000000000 --- a/modules/util/romio.py +++ /dev/null @@ -1,64 +0,0 @@ -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() diff --git a/modules/util/sss1day_FMB_py3.py b/modules/util/sss1day_FMB_py3.py deleted file mode 100644 index f05ebb2e9a390ba0267c1db387c16503d99ddc03..0000000000000000000000000000000000000000 --- a/modules/util/sss1day_FMB_py3.py +++ /dev/null @@ -1,308 +0,0 @@ -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()