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

Initial commit

parents
No related branches found
No related tags found
No related merge requests found
Showing
with 926 additions and 0 deletions
<?xml version="1.0" encoding="UTF-8"?>
<project version="4">
<component name="ProjectRootManager" version="2" project-jdk-name="Python 3.7 (modules)" project-jdk-type="Python SDK" />
</project>
\ No newline at end of file
<?xml version="1.0" encoding="UTF-8"?>
<module type="PYTHON_MODULE" version="4">
<component name="NewModuleRootManager">
<content url="file://$MODULE_DIR$" />
<orderEntry type="jdk" jdkName="Python 3.7 (modules)" jdkType="Python SDK" />
<orderEntry type="sourceFolder" forTests="false" />
</component>
<component name="TestRunnerService">
<option name="PROJECT_TEST_RUNNER" value="Unittests" />
</component>
</module>
\ No newline at end of file
<?xml version="1.0" encoding="UTF-8"?>
<project version="4">
<component name="ProjectModuleManager">
<modules>
<module fileurl="file://$PROJECT_DIR$/.idea/modules.iml" filepath="$PROJECT_DIR$/.idea/modules.iml" />
</modules>
</component>
</project>
\ No newline at end of file
<?xml version="1.0" encoding="UTF-8"?>
<project version="4">
<component name="ChangeListManager">
<list default="true" id="6685192a-aae3-4353-baf3-c54b18376dc3" name="Default Changelist" comment="" />
<option name="EXCLUDED_CONVERTED_TO_IGNORED" value="true" />
<option name="SHOW_DIALOG" value="false" />
<option name="HIGHLIGHT_CONFLICTS" value="true" />
<option name="HIGHLIGHT_NON_ACTIVE_CHANGELIST" value="false" />
<option name="LAST_RESOLUTION" value="IGNORE" />
</component>
<component name="FileEditorManager">
<leaf SIDE_TABS_SIZE_LIMIT_KEY="300">
<file pinned="false" current-in-tab="false">
<entry file="file://$PROJECT_DIR$/aeolus/aeolus_amv.py">
<provider selected="true" editor-type-id="text-editor">
<state relative-caret-position="150">
<caret line="10" column="16" lean-forward="true" selection-start-line="10" selection-start-column="16" selection-end-line="10" selection-end-column="16" />
<folding>
<element signature="e#0#19#0" expanded="true" />
</folding>
</state>
</provider>
</entry>
</file>
<file pinned="false" current-in-tab="false">
<entry file="file://$PROJECT_DIR$/deeplearning/amv_raob.py">
<provider selected="true" editor-type-id="text-editor">
<state relative-caret-position="10442">
<caret line="881" lean-forward="true" selection-start-line="881" selection-end-line="881" />
<folding>
<element signature="e#0#15#0" expanded="true" />
</folding>
</state>
</provider>
</entry>
</file>
<file pinned="false" current-in-tab="false">
<entry file="file://$PROJECT_DIR$/aeolus/geos_nav.py">
<provider selected="true" editor-type-id="text-editor">
<state relative-caret-position="30">
<caret line="2" selection-start-line="2" selection-end-line="2" />
</state>
</provider>
</entry>
</file>
<file pinned="false" current-in-tab="true">
<entry file="file://$PROJECT_DIR$/deeplearning/cloudheight.py">
<provider selected="true" editor-type-id="text-editor">
<state relative-caret-position="240">
<caret line="16" column="17" lean-forward="true" selection-start-line="16" selection-start-column="17" selection-end-line="16" selection-end-column="17" />
<folding>
<element signature="e#0#23#0" expanded="true" />
</folding>
</state>
</provider>
</entry>
</file>
<file pinned="false" current-in-tab="false">
<entry file="file://$PROJECT_DIR$/deeplearning/__init__.py">
<provider selected="true" editor-type-id="text-editor" />
</entry>
</file>
<file pinned="false" current-in-tab="false">
<entry file="file://$PROJECT_DIR$/amv/aeolus.py">
<provider selected="true" editor-type-id="text-editor">
<state relative-caret-position="105">
<caret line="7" lean-forward="true" selection-start-line="7" selection-end-line="7" />
</state>
</provider>
</entry>
</file>
<file pinned="false" current-in-tab="false">
<entry file="file://$PROJECT_DIR$/util/geos_nav.py">
<provider selected="true" editor-type-id="text-editor">
<state relative-caret-position="150">
<caret line="10" lean-forward="true" selection-start-line="10" selection-end-line="10" />
<folding>
<element signature="e#1#19#0" expanded="true" />
</folding>
</state>
</provider>
</entry>
</file>
<file pinned="false" current-in-tab="false">
<entry file="file://$PROJECT_DIR$/amv/intercompare.py">
<provider selected="true" editor-type-id="text-editor">
<state relative-caret-position="90">
<caret line="6" column="19" lean-forward="true" selection-start-line="6" selection-start-column="19" selection-end-line="6" selection-end-column="19" />
</state>
</provider>
</entry>
</file>
<file pinned="false" current-in-tab="false">
<entry file="file://$PROJECT_DIR$/util/bar_plot.py">
<provider selected="true" editor-type-id="text-editor">
<state>
<folding>
<element signature="e#0#18#0" expanded="true" />
</folding>
</state>
</provider>
</entry>
</file>
<file pinned="false" current-in-tab="false">
<entry file="file://$PROJECT_DIR$/util/goesr_l1b.py">
<provider selected="true" editor-type-id="text-editor">
<state relative-caret-position="315">
<caret line="21" lean-forward="true" selection-start-line="21" selection-end-line="21" />
<folding>
<element signature="e#430#444#0" expanded="true" />
</folding>
</state>
</provider>
</entry>
</file>
</leaf>
</component>
<component name="IdeDocumentHistory">
<option name="CHANGED_PATHS">
<list>
<option value="$PROJECT_DIR$/util/plot.py" />
<option value="$PROJECT_DIR$/util/lon_lat_grid.py" />
<option value="$PROJECT_DIR$/util/geos_nav.py" />
<option value="$PROJECT_DIR$/amv/intercompare.py" />
<option value="$PROJECT_DIR$/amv/aeolus.py" />
<option value="$PROJECT_DIR$/aeolus/aeolus_amv.py" />
<option value="$PROJECT_DIR$/aeolus/geos_nav.py" />
<option value="$PROJECT_DIR$/deeplearning/amv_raob.py" />
<option value="$PROJECT_DIR$/deeplearning/cloudheight.py" />
</list>
</option>
</component>
<component name="ProjectFrameBounds">
<option name="x" value="462" />
<option name="y" value="38" />
<option name="width" value="1743" />
<option name="height" value="1402" />
</component>
<component name="ProjectView">
<navigator proportions="" version="1">
<foldersAlwaysOnTop value="true" />
</navigator>
<panes>
<pane id="ProjectPane">
<subPane>
<expand>
<path>
<item name="modules" type="b2602c69:ProjectViewProjectNode" />
<item name="modules" type="462c0819:PsiDirectoryNode" />
</path>
<path>
<item name="modules" type="b2602c69:ProjectViewProjectNode" />
<item name="External Libraries" type="cb654da1:ExternalLibrariesNode" />
</path>
</expand>
<select />
</subPane>
</pane>
<pane id="Scope" />
</panes>
</component>
<component name="PropertiesComponent">
<property name="settings.editor.selected.configurable" value="com.jetbrains.python.configuration.PyActiveSdkModuleConfigurable" />
</component>
<component name="RunDashboard">
<option name="ruleStates">
<list>
<RuleState>
<option name="name" value="ConfigurationTypeDashboardGroupingRule" />
</RuleState>
<RuleState>
<option name="name" value="StatusDashboardGroupingRule" />
</RuleState>
</list>
</option>
</component>
<component name="SvnConfiguration">
<configuration />
</component>
<component name="TaskManager">
<task active="true" id="Default" summary="Default task">
<changelist id="6685192a-aae3-4353-baf3-c54b18376dc3" name="Default Changelist" comment="" />
<created>1599498568302</created>
<option name="number" value="Default" />
<option name="presentableId" value="Default" />
<updated>1599498568302</updated>
</task>
<servers />
</component>
<component name="ToolWindowManager">
<frame x="462" y="38" width="1743" height="1402" extended-state="0" />
<editor active="true" />
<layout>
<window_info content_ui="combo" id="Project" order="0" visible="true" weight="0.24970964" />
<window_info id="Structure" order="1" side_tool="true" weight="0.25" />
<window_info id="Favorites" order="2" side_tool="true" />
<window_info anchor="bottom" id="Message" order="0" />
<window_info anchor="bottom" id="Find" order="1" />
<window_info anchor="bottom" id="Run" order="2" />
<window_info anchor="bottom" id="Debug" order="3" weight="0.4" />
<window_info anchor="bottom" id="Cvs" order="4" weight="0.25" />
<window_info anchor="bottom" id="Inspection" order="5" weight="0.4" />
<window_info anchor="bottom" id="TODO" order="6" />
<window_info anchor="bottom" id="Version Control" order="7" />
<window_info anchor="bottom" id="Terminal" order="8" />
<window_info anchor="bottom" id="Event Log" order="9" side_tool="true" />
<window_info active="true" anchor="bottom" id="Python Console" order="10" visible="true" weight="0.32977098" />
<window_info anchor="right" id="Commander" internal_type="SLIDING" order="0" type="SLIDING" weight="0.4" />
<window_info anchor="right" id="Ant Build" order="1" weight="0.25" />
<window_info anchor="right" content_ui="combo" id="Hierarchy" order="2" weight="0.25" />
</layout>
</component>
<component name="editorHistoryManager">
<entry file="file://$PROJECT_DIR$/util/util.py">
<provider selected="true" editor-type-id="text-editor" />
</entry>
<entry file="file://$PROJECT_DIR$/util/setup.py">
<provider selected="true" editor-type-id="text-editor">
<state relative-caret-position="15">
<caret line="1" lean-forward="true" selection-start-line="1" selection-end-line="1" />
<folding>
<element signature="e#40#105#0" expanded="true" />
</folding>
</state>
</provider>
</entry>
<entry file="file://$PROJECT_DIR$/util/plot_cm.py">
<provider selected="true" editor-type-id="text-editor">
<state>
<folding>
<element signature="e#0#25#0" expanded="true" />
</folding>
</state>
</provider>
</entry>
<entry file="file://$PROJECT_DIR$/util/plot.py">
<provider selected="true" editor-type-id="text-editor">
<state relative-caret-position="195">
<caret line="13" lean-forward="true" selection-start-line="13" selection-end-line="13" />
</state>
</provider>
</entry>
<entry file="file://$PROJECT_DIR$/util/lon_lat_grid.py">
<provider selected="true" editor-type-id="text-editor">
<state relative-caret-position="105">
<caret line="7" column="38" lean-forward="true" selection-start-line="7" selection-start-column="38" selection-end-line="7" selection-end-column="38" />
</state>
</provider>
</entry>
<entry file="file://$PROJECT_DIR$/util/goesr_l1b.py">
<provider selected="true" editor-type-id="text-editor">
<state relative-caret-position="315">
<caret line="21" lean-forward="true" selection-start-line="21" selection-end-line="21" />
<folding>
<element signature="e#430#444#0" expanded="true" />
</folding>
</state>
</provider>
</entry>
<entry file="file://$PROJECT_DIR$/util/geos_nav.py">
<provider selected="true" editor-type-id="text-editor">
<state relative-caret-position="150">
<caret line="10" lean-forward="true" selection-start-line="10" selection-end-line="10" />
<folding>
<element signature="e#1#19#0" expanded="true" />
</folding>
</state>
</provider>
</entry>
<entry file="file://$PROJECT_DIR$/aeolus/aeolus_amv.py">
<provider selected="true" editor-type-id="text-editor">
<state relative-caret-position="150">
<caret line="10" column="16" lean-forward="true" selection-start-line="10" selection-start-column="16" selection-end-line="10" selection-end-column="16" />
<folding>
<element signature="e#0#19#0" expanded="true" />
</folding>
</state>
</provider>
</entry>
<entry file="file://$PROJECT_DIR$/amv/intercompare.py">
<provider selected="true" editor-type-id="text-editor">
<state relative-caret-position="90">
<caret line="6" column="19" lean-forward="true" selection-start-line="6" selection-start-column="19" selection-end-line="6" selection-end-column="19" />
</state>
</provider>
</entry>
<entry file="file://$PROJECT_DIR$/amv/aeolus.py">
<provider selected="true" editor-type-id="text-editor">
<state relative-caret-position="105">
<caret line="7" lean-forward="true" selection-start-line="7" selection-end-line="7" />
</state>
</provider>
</entry>
<entry file="file://$PROJECT_DIR$/util/bar_plot.py">
<provider selected="true" editor-type-id="text-editor">
<state>
<folding>
<element signature="e#0#18#0" expanded="true" />
</folding>
</state>
</provider>
</entry>
<entry file="file://$PROJECT_DIR$/aeolus/geos_nav.py">
<provider selected="true" editor-type-id="text-editor">
<state relative-caret-position="30">
<caret line="2" selection-start-line="2" selection-end-line="2" />
</state>
</provider>
</entry>
<entry file="file://$PROJECT_DIR$/deeplearning/__init__.py">
<provider selected="true" editor-type-id="text-editor" />
</entry>
<entry file="file://$PROJECT_DIR$/deeplearning/amv_raob.py">
<provider selected="true" editor-type-id="text-editor">
<state relative-caret-position="10442">
<caret line="881" lean-forward="true" selection-start-line="881" selection-end-line="881" />
<folding>
<element signature="e#0#15#0" expanded="true" />
</folding>
</state>
</provider>
</entry>
<entry file="file://$PROJECT_DIR$/deeplearning/cloudheight.py">
<provider selected="true" editor-type-id="text-editor">
<state relative-caret-position="240">
<caret line="16" column="17" lean-forward="true" selection-start-line="16" selection-start-column="17" selection-end-line="16" selection-end-column="17" />
<folding>
<element signature="e#0#23#0" expanded="true" />
</folding>
</state>
</provider>
</entry>
</component>
</project>
\ No newline at end of file
import datetime, os
from datetime import timezone
import glob
import numpy as np
from netCDF4 import Dataset, Dimension, Variable
from aeolus.geos_nav import GEOSNavigation
amv_file_duration = 60 # minutes
half_width = 20 # search box centered on AEOLUS profile (FGF coordinates)
num_elems = 5424
num_lines = 5424
first_time = True
ftimes = []
flist = None
G16_lon_range = [-130.0, -30.0]
G16_lat_range = [-64.0, 64.0]
class MyGenericException(Exception):
def __init__(self, message):
self.message = message
def get_datetime(pathname):
fname = os.path.split(pathname)[1]
toks = fname.split('_')
dstr = toks[4]
tstr = toks[5]
dtstr = dstr+tstr
dto = datetime.datetime.strptime(dtstr, '%Y%j%H%M').replace(tzinfo=timezone.utc)
return dto
def get_datetime_ops(pathname):
fname = os.path.split(pathname)[1]
toks = fname.split('_')
dtstr = toks[3]
dtstr = dtstr[:-3]
dto = datetime.datetime.strptime(dtstr, 's%Y%j%H%M').replace(tzinfo=timezone.utc)
return dto
def get_file_containing_time(timestamp, files_path, file_time_span, amv_source='OPS', band='14'):
global first_time, ftimes, flist
if first_time is True:
if amv_source == 'OPS':
flist = glob.glob(files_path + 'OR_ABI-L2-DMWF*'+'C'+band+'*.nc')
dto_func = get_datetime_ops
elif amv_source == 'SSEC':
flist = glob.glob(files_path + '*WINDS_AMV_EN-'+band+'*.nc')
dto_func = get_datetime
for pname in flist: # TODO: make better with regular expressions (someday)
dto = dto_func(pname)
dto_start = dto
dto_end = dto + datetime.timedelta(minutes=file_time_span)
ftimes.append((dto_start.timestamp(), dto_end.timestamp()))
first_time = False
k = -1
for i in range(len(ftimes)):
if (timestamp >= ftimes[i][0]) and (timestamp < ftimes[i][1]):
k = i
break
if k < 0:
return None, None, None
return flist[k], ftimes[k], k
# imports the S4 NOAA output
# filename: full path as a string, '/home/user/filename'
# returns a dict: time -> list of profiles (a profile is a list of levels)
def get_aeolus_time_dict(filename, lon360=False):
time_dict = {}
with open(filename) as file:
while True:
prof_hdr = file.readline()
if not prof_hdr:
break
toks = prof_hdr.split()
yr = int(float(toks[0]))
mon = int(float(toks[1]))
dy = int(float(toks[2]))
hr = int(float(toks[3]))
mn = int(float(toks[4]))
ss = int(float(toks[5]))
lon = float(toks[6])
lat = float(toks[7])
nlevs = int(toks[8])
if lon360:
if lon < 0:
lon += 360.0
else:
if lon > 180.0:
lon -= 360.0
dto = datetime.datetime(year=yr, month=mon, day=dy, hour=hr, minute=mn, second=ss)
dto = dto.replace(tzinfo=timezone.utc)
timestamp = dto.timestamp()
prof = []
if time_dict.get(timestamp, -1) == -1:
prof_s = []
prof_s.append(prof)
time_dict[timestamp] = prof_s
else:
prof_s = time_dict.get(timestamp)
prof_s.append(prof)
for k in range(nlevs):
line = file.readline()
toks = line.split()
lvlidx = int(toks[0])
hhh = float(toks[1]) * 1000.0
hht = float(toks[2]) * 1000.0
hhb = float(toks[3]) * 1000.0
err = float(toks[4])
azm = float(toks[5])
ws = float(toks[6])
len = float(toks[7])
tup = (lat, lon, hhh, hht, hhb, azm, ws)
prof.append(tup)
file.close()
return time_dict
# make each profile at a timestamp a numpy array
def time_dict_to_nd(time_dict):
keys = list(time_dict.keys())
for key in keys:
vals = time_dict[key]
if vals is not None:
for i in range(len(vals)):
nda = np.array(vals[i])
vals[i] = nda
return time_dict
def get_search_box(nav, lon, lat):
cc, ll = nav.earth_to_lc(lon, lat)
c_rng = [cc - half_width, cc + half_width]
l_rng = [ll - half_width, ll + half_width]
if c_rng[0] < 0:
c_rng[0] = 0
if l_rng[0] < 0:
l_rng[0] = 0
if c_rng[1] >= num_elems:
c_rng[1] = num_elems - 1
if l_rng[1] >= num_lines:
l_rng[1] = num_lines - 1
return c_rng, l_rng
# aeolus_dict: time -> profiles
# amv_files_path: directory containing AMVs, '/home/user/amvdir/'
# return dict: aeolus time -> tuple (amv_lon, amv_lat, amv_pres, amv_spd, amv_dir)
def match_amvs_to_aelous(aeolus_dict, amv_files_path, amv_source='OPS', band='C14'):
nav = GEOSNavigation()
match_dict = {}
keys = list(aeolus_dict.keys())
last_f_idx = -1
for key in keys:
fname, ftime, f_idx = get_file_containing_time(key, amv_files_path, amv_file_duration, amv_source, band)
if f_idx is None:
continue
profs = aeolus_dict.get(key)
layers = profs[0]
if layers is None:
continue
lat = layers[0, 0]
lon = layers[0, 1]
# check range
if lat > G16_lat_range[1] or lat < G16_lat_range[0]:
continue
if lon > G16_lon_range[1] or lon < G16_lon_range[0]:
continue
# ------------------------------------------------
if f_idx != last_f_idx:
last_f_idx = f_idx
ds = Dataset(fname)
amv_lons = ds['Longitude'][:]
amv_lats = ds['Latitude'][:]
amv_spd = ds['Wind_Speed'][:]
amv_dir = ds['Wind_Dir'][:]
amv_pres = ds['MedianPress'][:]
cc = ds['Element'][:]
ll = ds['Line'][:]
# cc, ll = nav.earth_to_lc_s(amv_lons, amv_lats)
ds.close()
c_rng, l_rng = get_search_box(nav, lon, lat)
in_cc = np.logical_and(cc > c_rng[0], cc < c_rng[1])
in_ll = np.logical_and(ll > l_rng[0], ll < l_rng[1])
in_box = np.logical_and(in_cc, in_ll)
num_amvs = np.sum(in_box)
if num_amvs == 0:
continue
match_dict[key] = (amv_lons[in_box], amv_lats[in_box], amv_pres[in_box], amv_spd[in_box], amv_dir[in_box])
return match_dict
# full path as string filename to create, '/home/user/newfilename'
# aeolus_to_amv_dct: output from match_amvs_to_aeolus
# aeolus_dct: output from get_aeolus_time_dict
def create_file(filename, aeolus_to_amv_dct, aeolus_dct):
keys = list(aeolus_to_amv_dct.keys())
num_amvs = []
num_levs = []
times = []
namvs = 0
nlevs = 0
for key in keys:
lons, lats, pres, spd, dir = aeolus_to_amv_dct.get(key)
num_amvs.append(len(lons))
namvs += len(lons)
prof_s = aeolus_dct.get(key)
prof = prof_s[0]
num_levs.append(prof.shape[0])
nlevs += prof.shape[0]
times.append(key)
amv_per_alus = len(aeolus_to_amv_dct)
rootgrp = Dataset(filename, 'w', format='NETCDF4')
dim_amvs = rootgrp.createDimension('amvs', size=namvs)
dim_alus = rootgrp.createDimension('profs', size=nlevs)
dim_num_aeolus_prof = rootgrp.createDimension('num_aeolus_profs', size=len(aeolus_to_amv_dct))
amv_lon = rootgrp.createVariable('amv_longitude', 'f4', ['amvs'])
amv_lon.units = 'degrees east'
amv_lat = rootgrp.createVariable('amv_latitude', 'f4', ['amvs'])
amv_lat.units = 'degrees north'
amv_spd = rootgrp.createVariable('amv_spd', 'f4', ['amvs'])
amv_spd.units = 'm s-1'
amv_dir = rootgrp.createVariable('amv_dir', 'f4', ['amvs'])
amv_dir.units = 'degree'
amv_pres = rootgrp.createVariable('amv_pres', 'f4', ['amvs'])
amv_pres.units = 'hPa'
num_amvs_per_prof = rootgrp.createVariable('num_amvs_per_prof', 'i4', ['num_aeolus_profs'])
num_levs_per_prof = rootgrp.createVariable('num_levs_per_prof', 'i4', ['num_aeolus_profs'])
prof_time = rootgrp.createVariable('time', 'f4', ['num_aeolus_profs'])
prof_time.units = 'seconds since 1970-01-1 00:00:00'
prf_azm = rootgrp.createVariable('prof_azm', 'f4', ['profs'])
prf_azm.units = 'degree'
prf_spd = rootgrp.createVariable('prof_spd', 'f4', ['profs'])
prf_spd.units = 'm s-1'
prf_hht = rootgrp.createVariable('prof_hht', 'f4', ['profs'])
prf_hht.units = 'meter'
prf_hhb = rootgrp.createVariable('prof_hhb', 'f4', ['profs'])
prf_hhb.units = 'meter'
i_a = 0
i_c = 0
for idx, key in enumerate(keys):
namvs = num_amvs[idx]
nlevs = num_levs[idx]
i_b = i_a + namvs
i_d = i_c + nlevs
lons, lats, pres, spd, dir = aeolus_to_amv_dct.get(key)
amv_lon[i_a:i_b] = lons[:]
amv_lat[i_a:i_b] = lats[:]
amv_spd[i_a:i_b] = spd[:]
amv_dir[i_a:i_b] = dir[:]
amv_pres[i_a:i_b] = pres[:]
i_a += namvs
prof_s = aeolus_dct.get(key)
prof = prof_s[0]
prf_hht[i_c:i_d] = prof[:, 2]
prf_hhb[i_c:i_d] = prof[:, 3]
prf_azm[i_c:i_d] = prof[:, 5]
prf_spd[i_c:i_d] = prof[:, 6]
i_c += nlevs
num_amvs_per_prof[:] = num_amvs
num_levs_per_prof[:] = num_levs
prof_time[:] = times
rootgrp.close()
# def create_file(filename, aeolus_to_amv_dct, aeolus_dct):
# keys = list(aeolus_to_amv_dct.keys())
#
# num_amvs = []
# num_levs = []
# times = []
#
# namvs = 0
# nlevs = 0
# for key in keys:
# lons, lats, pres, spd, dir = aeolus_to_amv_dct.get(key)
# num_amvs.append(len(lons))
# namvs += len(lons)
#
# prof_s = aeolus_dct.get(key)
# prof = prof_s[0]
# num_levs.append(prof.shape[0])
# nlevs += prof.shape[0]
#
# times.append(key)
#
# amv_per_alus = len(aeolus_to_amv_dct)
# rootgrp = Dataset(filename, 'w', format='NETCDF4')
# dim_amvs = rootgrp.createDimension('amvs', size=namvs)
# dim_alus = rootgrp.createDimension('profs', size=nlevs)
# dim_num_aeolus_prof = rootgrp.createDimension('num_aeolus_profs', size=len(aeolus_to_amv_dct))
#
# amv_lon = rootgrp.createVariable('amv_longitude', 'f4', ['amvs'])
# amv_lat = rootgrp.createVariable('amv_latitude', 'f4', ['amvs'])
# amv_spd = rootgrp.createVariable('amv_spd', 'f4', ['amvs'])
# amv_dir = rootgrp.createVariable('amv_dir', 'f4', ['amvs'])
# amv_pres = rootgrp.createVariable('amv_pres', 'f4', ['amvs'])
#
# num_amvs_per_prof = rootgrp.createVariable('num_amvs_per_prof', 'i4', ['num_aeolus_profs'])
# num_levs_per_prof = rootgrp.createVariable('num_levs_per_prof', 'i4', ['num_aeolus_profs'])
# prof_time = rootgrp.createVariable('time', 'f4', ['num_aeolus_profs'])
#
# prf_azm = rootgrp.createVariable('prof_azm', 'f4', ['profs'])
# prf_spd = rootgrp.createVariable('prof_spd', 'f4', ['profs'])
# prf_hht = rootgrp.createVariable('prof_hht', 'f4', ['profs'])
# prf_hhb = rootgrp.createVariable('prof_hhb', 'f4', ['profs'])
#
# i_a = 0
# i_c = 0
# for idx, key in enumerate(keys):
# namvs = num_amvs[idx]
# nlevs = num_levs[idx]
# i_b = i_a + namvs
# i_d = i_c + nlevs
#
# lons, lats, pres, spd, dir = aeolus_to_amv_dct.get(key)
# amv_lon[i_a:i_b] = lons[:]
# amv_lat[i_a:i_b] = lats[:]
# amv_spd[i_a:i_b] = spd[:]
# amv_dir[i_a:i_b] = dir[:]
# amv_pres[i_a:i_b] = pres[:]
# i_a += namvs
#
# prof_s = aeolus_dct.get(key)
# prof = prof_s[0]
# prf_hht[i_c:i_d] = prof[:, 2]
# prf_hhb[i_c:i_d] = prof[:, 3]
# prf_azm[i_c:i_d] = prof[:, 5]
# prf_spd[i_c:i_d] = prof[:, 6]
# i_c += nlevs
#
# num_amvs_per_prof[:] = num_amvs
# num_levs_per_prof[:] = num_levs
# prof_time[:] = times
#
# rootgrp.close()
# aeolus_file: S4 NOAA txt output
# amv_files_dir: G16/17 AMV product file
# outfile: pathname for the Netcdf match file
def run(aeolus_file, amv_files_dir, outfile=None, amv_source='OPS', band='14'):
a_d = get_aeolus_time_dict(aeolus_file)
a_d = time_dict_to_nd(a_d)
m_d = match_amvs_to_aelous(a_d, amv_files_dir, amv_source, band)
if outfile is not None:
create_file(outfile, m_d, a_d)
import numpy as np
class MyGenericException(Exception):
def __init__(self, message):
self.message = message
DEG_TO_RAD = np.pi/180.0;
RAD_TO_DEG = 180.0/np.pi;
r_pol = 6356.5838 # km
r_eq = 6378.1690 # km
# AHI -------------------------------
# h = 42164.0 # barycentric height, km
# sub_lon = 140.7
# sub_lon *= DEG_TO_RAD
# scan_geom = 'GEOS'
# signs modified from static nav file to work with ABI FGF code
# CFAC = 5.588799E-05
# LFAC = -5.588799E-05
# COFF = -0.1537199
# LOFF = 0.1537199
# computed from static nav file (lon, lat) grid
# Note 0-based to the pixel center
# CFAC = 5.58924125e-05
# LFAC = -5.58810490e-05
# COFF = -1.53678977e-01
# LOFF = 1.53644345e-01
# GOES ------------------------------
h = 42164.16 # barycentric height, km
scan_geom = 'GOES'
sub_lon = -75.0
sub_lon *= DEG_TO_RAD
# official for FD
CFAC = 5.6E-05
LFAC = -5.6E-05
COFF = -0.151844
LOFF = 0.151844
# official for CONUS
#CFAC = 5.6E-05
#COFF = -0.101332
#LFAC = -5.6E-05
#LOFF = 0.128212
# computed for CLAVRx FD
#CFAC = 5.60016368e-05
#LFAC = -5.59941969e-05
#COFF = -1.51780260e-01
#LOFF = 1.51773560e-01
# 65536 = 2^16
class GEOSNavigation:
def __init__(self, sub_lon=-75.0, barycentric_height=42164.16, scan_geom='GOES', CFAC=5.6E-05, LFAC=-5.6E-05, COFF=-0.151844, LOFF=0.151844, num_elems=5424, num_lines=5424):
self.sub_lon = sub_lon * DEG_TO_RAD
self.h = barycentric_height
self.scan_geom = scan_geom
self.CFAC = CFAC
self.LFAC = LFAC
self.COFF = COFF
self.LOFF = LOFF
self.num_lines = num_lines
self.num_elems = num_elems
def earth_to_sat(self, geographic_lon, geographic_lat):
geographic_lon *= DEG_TO_RAD
geographic_lat *= DEG_TO_RAD
geocentric_lat = np.arctan(((r_pol*r_pol)/(r_eq*r_eq))*np.tan(geographic_lat))
r_earth = r_pol / np.sqrt(1.0 - ((r_eq * r_eq - r_pol * r_pol) / (r_eq * r_eq)) * np.cos(geocentric_lat) * np.cos(geocentric_lat))
r_1 = self.h - r_earth * np.cos(geocentric_lat) * np.cos(geographic_lon - self.sub_lon)
r_2 = -r_earth * np.cos(geocentric_lat) * np.sin(geographic_lon - self.sub_lon)
r_3 = r_earth * np.sin(geocentric_lat)
if r_1 > self.h:
return np.nan, np.nan
if self.scan_geom == 'GEOS':
lamda_sat = np.arctan(-r_2/r_1)
theta_sat = np.arcsin(r_3/np.sqrt(r_1*r_1 + r_2*r_2 + r_3*r_3))
elif self.scan_geom == 'GOES':
lamda_sat = np.arcsin(-r_2/np.sqrt(r_1*r_1 + r_2*r_2 + r_3*r_3))
theta_sat = np.arctan(r_3/r_1)
return lamda_sat, theta_sat
def earth_to_sat_s(self, geographic_lon, geographic_lat):
geographic_lon = geographic_lon * DEG_TO_RAD
geographic_lat = geographic_lat * DEG_TO_RAD
geocentric_lat = np.arctan(((r_pol*r_pol)/(r_eq*r_eq))*np.tan(geographic_lat))
r_earth = r_pol / np.sqrt(1.0 - ((r_eq * r_eq - r_pol * r_pol) / (r_eq * r_eq)) * (np.cos(geocentric_lat))**2)
r_1 = self.h - r_earth * np.cos(geocentric_lat) * np.cos(geographic_lon - self.sub_lon)
r_2 = -r_earth * np.cos(geocentric_lat) * np.sin(geographic_lon - self.sub_lon)
r_3 = r_earth * np.sin(geocentric_lat)
if self.scan_geom == 'GEOS':
lamda_sat = np.arctan(-r_2/r_1)
theta_sat = np.arcsin(r_3/np.sqrt(r_1*r_1 + r_2*r_2 + r_3*r_3))
elif self.scan_geom == 'GOES':
lamda_sat = np.arcsin(-r_2/np.sqrt(r_1*r_1 + r_2*r_2 + r_3*r_3))
theta_sat = np.arctan(r_3/r_1)
np.where(r_1 > self.h, lamda_sat == np.nan, lamda_sat)
np.where(r_1 > self.h, theta_sat == np.nan, theta_sat)
return lamda_sat, theta_sat
def sat_to_lc(self, lamda, theta):
# c = COFF + np.rint(lamda * (1/np.power(2, 16)) * CFAC)
# l = LOFF + np.rint(theta * (1/np.power(2, 16)) * LFAC)
# (float(2 ** 16) * (float(lc) - off)) / float(fac) * DEG_TO_RAD
cc = (lamda - self.COFF) / self.CFAC
ll = (theta - self.LOFF) / self.LFAC
cc = np.floor(cc + 0.5)
ll = np.floor(ll + 0.5)
cc = cc.astype(np.int32)
ll = ll.astype(np.int32)
return cc, ll
def earth_to_lc(self, lon, lat):
lamda, theta = self.earth_to_sat(lon, lat)
if np.isnan(lamda):
return None, None
cc, ll = self.sat_to_lc(lamda, theta)
return cc, ll
def earth_to_lc_s(self, lons, lats):
lamda, theta = self.earth_to_sat_s(lons, lats)
cc, ll = self.sat_to_lc(lamda, theta)
np.where(np.isnan(cc), cc == -1, cc)
np.where(np.isnan(ll), ll == -1, ll)
return cc, ll
def earth_to_indexs(self, lons, lats, len_x):
num = lons.shape[0]
idxs = []
for i in range(num):
cc, ll = self.earth_to_lc(lons[i], lats[i])
if cc is None:
idxs.append(-1)
else:
idx = ll * len_x + cc
idxs.append(idx)
idxs = np.array(idxs)
return idxs
This diff is collapsed.
This diff is collapsed.
File added
File added
File added
This diff is collapsed.
This diff is collapsed.
File added
File added
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment