diff --git a/config/cfg.py b/config/cfg.py
index 0c8131dd9a6df53dc9778b05da477c81572c147a..f48bf9fd770b92cad47396b6ee6f85c8499abd52 100755
--- a/config/cfg.py
+++ b/config/cfg.py
@@ -7,37 +7,44 @@ class ExperimentConfiguration(object):
pass
exp = ExperimentConfiguration()
-exp.expname = "exp_v1.19_rr_WV_obs4-20_loc10"
+exp.expname = "exp_v1.21_P3_wbub7_WV+VIS_obs10_loc20"
exp.model_dx = 2000
exp.n_ens = 40
exp.n_nodes = 10
+#exp.nature_wrfout = '/home/fs71386/lkugler/data/sim_archive/exp_v1.19_P5+su_nat2/2008-07-30_07:00/1/wrfout_d01_%Y-%m-%d_%H:%M:%S'
+exp.nature_wrfout = '/home/fs71386/lkugler/data/sim_archive/exp_v1.19_P3_wbub7_nat/2008-07-30_12:00/1/wrfout_d01_%Y-%m-%d_%H:%M:%S'
#exp.nature_wrfout = '/home/fs71386/lkugler/data/sim_archive/exp_v1.19_Pwbub5_nat/2008-07-30_12:00/1/wrfout_d01_%Y-%m-%d_%H:%M:%S'
-exp.nature_wrfout = '/home/fs71386/lkugler/data/sim_archive/exp_v1.18_P1_nature/2008-07-30_06:00/1/wrfout_d01_%Y-%m-%d_%H:%M:%S'
+#exp.nature_wrfout = '/home/fs71386/lkugler/data/sim_archive/exp_v1.18_P1_nature/2008-07-30_06:00/1/wrfout_d01_%Y-%m-%d_%H:%M:%S'
+#exp.nature_wrfout = '/home/fs71386/lkugler/data/sim_archive/exp_v1.19_P4_nat/2008-07-30_07:00/1/wrfout_d01_%Y-%m-%d_%H:%M:%S'
#exp.input_profile = '/home/fs71386/lkugler/wrf_profiles/data/wrf/ens/2021-05-04/raso.nat.001.wrfprof'
#exp.input_profile = '/home/fs71386/lkugler/wrf_profiles/data/wrf/ens/2021-05-04/raso.nat.<iens>.wrfprof'
-exp.input_profile = '/home/fs71386/lkugler/wrf_profiles/data/wrf/ens/2021-05-04/raso.fc.<iens>.wrfprof'
-#exp.input_profile = '/home/fs71386/lkugler/wrf_profiles/data/wrf/ens/improved_pert10/raso.fc.<iens>.wrfprof'
+#exp.input_profile = '/home/fs71386/lkugler/wrf_profiles/data/wrf/ens/2021-05-04/raso.fc.<iens>.wrfprof'
+#exp.input_profile = '/home/fs71386/lkugler/data/initial_profiles/wrf/ens/large_mean_error/raso.nat.<iens>.wrfprof'
+exp.input_profile = '/gpfs/data/fs71386/lkugler/initial_profiles/wrf/ens/2022-03-31/raso.fc.<iens>.wrfprof'
+#exp.input_profile = '/gpfs/data/fs71386/lkugler/initial_profiles/wrf/ens/2022-03-31/raso.nat.<iens>.wrfprof'
+#exp.input_profile = '/gpfs/data/fs71386/lkugler/initial_profiles/wrf/ens/2022-05-18_nopert/raso.nat.<iens>.wrfprof'
+#exp.input_profile = '/gpfs/data/fs71386/lkugler/initial_profiles/wrf/ens/2022-05-18/raso.fc.<iens>.wrfprof'
# localize vertically, if it has a vertical position
# needs a horizontal scale too, to calculate the vertical normalization
# since you can not specify different vertical localizations for diff. variables
-#exp.cov_loc_vert_km_horiz_km = (2, 20)
-exp.superob_km = 20
+exp.cov_loc_vert_km_horiz_km = (2, 20)
+#exp.superob_km = 10
-n_obs = 5776 # 5776: 4km, 121: 30km, 256:16x16 (20km); 961: 10km resoltn # radar: n_obs for each observation height level
+n_obs = 961 # 22500: 2km, 5776: 4km, 121: 30km, 256:16x16 (20km); 961: 10km resoltn # radar: n_obs for each observation height level
vis = dict(plotname='VIS 0.6µm', plotunits='[1]',
kind='MSG_4_SEVIRI_BDRF', sat_channel=1, n_obs=n_obs,
error_generate=0.03, error_assimilate=0.06,
- cov_loc_radius_km=10)
+ cov_loc_radius_km=20)
wv73 = dict(plotname='Brightness temperature WV 7.3µm', plotunits='[K]',
kind='MSG_4_SEVIRI_TB', sat_channel=6, n_obs=n_obs,
- error_generate=1., error_assimilate=3.,
- cov_loc_radius_km=10)
+ error_generate=1., error_assimilate=False,
+ cov_loc_radius_km=20)
ir108 = dict(plotname='Brightness temperature IR 10.8µm', plotunits='[K]',
kind='MSG_4_SEVIRI_TB', sat_channel=9, n_obs=n_obs,
@@ -46,15 +53,15 @@ ir108 = dict(plotname='Brightness temperature IR 10.8µm', plotunits='[K]',
radar = dict(plotname='Radar reflectivity', plotunits='[dBz]',
kind='RADAR_REFLECTIVITY', n_obs=n_obs,
- error_generate=2.5, error_assimilate=5.,
- heights=np.arange(1000, 15001, 1000),
- cov_loc_radius_km=10)
+ error_generate=2.5, error_assimilate=5,
+ heights=np.arange(2000, 14001, 2000),
+ cov_loc_radius_km=20)
t = dict(plotname='Temperature', plotunits='[K]',
kind='RADIOSONDE_TEMPERATURE', n_obs=n_obs,
- error_generate=0.25, error_assimilate=0.5,
- heights=np.arange(1000, 15001, 1000),
- cov_loc_radius_km=10)
+ error_generate=0.2, error_assimilate=0.4,
+ heights=np.arange(1000, 17001, 2000),
+ cov_loc_radius_km=20)
t2m = dict(plotname='SYNOP Temperature', plotunits='[K]',
kind='SYNOP_TEMPERATURE', n_obs=n_obs,
@@ -67,8 +74,9 @@ psfc = dict(plotname='SYNOP Pressure', plotunits='[dBz]',
cov_loc_radius_km=32)
-exp.observations = [wv73]
-exp.update_vars = ['U', 'V', 'T', 'PH', 'MU', 'QVAPOR', 'QCLOUD', 'QICE', 'CLDFRA', 'PSFC']
+exp.observations = [vis, wv73]
+exp.update_vars = ['U', 'V', 'W', 'T', 'THM', 'PH', 'MU', 'QVAPOR', 'QCLOUD', 'QICE', 'PSFC']
+#exp.update_vars = ['U', 'V', 'W', 'T', 'PH', 'MU', 'QVAPOR', 'PSFC']
# directory paths depend on the name of the experiment
cluster.expname = exp.expname
diff --git a/config/clusters.py b/config/clusters.py
index 72baf5e01dbe986bc113d8306d3fbfbcccc80ecd..5b84038829bd3834617977b1f3305bd57da92709 100755
--- a/config/clusters.py
+++ b/config/clusters.py
@@ -77,7 +77,7 @@ vsc.archive_base = '/gpfs/data/fs71386/lkugler/sim_archive/'
# paths used as input
vsc.srcdir = '/gpfs/data/fs71386/lkugler/compile/WRF/WRF-4.3/run'
-vsc.dart_srcdir = '/gpfs/data/fs71386/lkugler/compile/DART/DART-9.11.9/models/wrf/work'
+vsc.dart_srcdir = '/gpfs/data/fs71386/lkugler/compile/DART/DART/models/wrf/work'
vsc.rttov_srcdir = '/gpfs/data/fs71386/lkugler/compile/RTTOV13/rtcoef_rttov13/'
vsc.scriptsdir = '/home/fs71386/lkugler/DART-WRF/dartwrf/'
diff --git a/dartwrf/obsseq.py b/dartwrf/obsseq.py
index 208510d67ec8a3cc6eecb9b937ad0e9c03cbe26d..e9340edb738347b7b45586837ef4aa57702fb251 100755
--- a/dartwrf/obsseq.py
+++ b/dartwrf/obsseq.py
@@ -1,3 +1,8 @@
+"""Read, modify and save DART obs_seq files.
+
+Not usable for creating obs_seq, since it does not know which metadata is necessary for each type
+"""
+
import os, sys, shutil, warnings
import numpy as np
import pandas as pd
diff --git a/dartwrf/utils.py b/dartwrf/utils.py
index 89392b5cb3dbe5fc2d3b57135b8cd0e86220783f..10d178a8a0fa78f65da0455b184752fef1fabb33 100755
--- a/dartwrf/utils.py
+++ b/dartwrf/utils.py
@@ -1,6 +1,10 @@
import os, sys, shutil, glob
import builtins as __builtin__
#copy = shutil.copy
+import subprocess
+
+def shell(cmd):
+ subprocess.check_output(cmd.split(' '), shell=True)
def print(*args):
__builtin__.print(*args, flush=True)
@@ -21,9 +25,6 @@ def try_remove(f):
def mkdir(path):
os.system('mkdir -p '+path)
-def mkdir_srvx8(path):
- os.system('ssh a1254888@srvx8.img.univie.ac.at mkdir -p '+path)
-
def script_to_str(path):
return open(path, 'r').read()
diff --git a/scheduler.py b/scheduler.py
index 3885d0f911c4ff87ab11aa1ceb2518270503fa0f..9367f68828b3033e27a14a4d6a52bba374095466 100755
--- a/scheduler.py
+++ b/scheduler.py
@@ -52,10 +52,10 @@ def backup_scripts():
def prepare_WRFrundir(init_time):
"""Create WRF/run directories and wrfinput files
"""
- s = my_Slurm("prep_wrfrundir", cfg_update={"time": "5", "mail-type": "BEGIN"})
- id = s.run(cluster.python+' '+cluster.scripts_rundir+'/prepare_wrfrundir.py '
- +init_time.strftime('%Y-%m-%d_%H:%M'))
- return id
+ cmd = cluster.python+' '+cluster.scripts_rundir+'/prepare_wrfrundir.py '+init_time.strftime('%Y-%m-%d_%H:%M')
+ print(cmd)
+ os.system(cmd)
+ # return id
def run_ideal(depends_on=None):
"""Run ideal for every ensemble member"""
@@ -200,7 +200,7 @@ def update_IC_from_DA(assim_time, depends_on=None):
def create_satimages(init_time, depends_on=None):
- s = my_Slurm("pRTTOV", cfg_update={"ntasks": "48", "time": "60", "nodes": "1"})
+ s = my_Slurm("pRTTOV", cfg_update={"ntasks": "48", "time": "80", "nodes": "1"})
id = s.run('/home/fs71386/lkugler/RTTOV-WRF/withmodules /home/fs71386/lkugler/RTTOV-WRF/run_init.py '+cluster.archivedir
+init_time.strftime('/%Y-%m-%d_%H:%M/'),
depends_on=[depends_on])
@@ -220,12 +220,23 @@ def gen_obsseq(depends_on=None):
return id
-def verify(depends_on=None):
- s = my_Slurm("verify-"+exp.expname, cfg_update={"time": "240", "mail-type": "FAIL,END", "ntasks": "1",
+def verify_sat(depends_on=None):
+ s = my_Slurm("verif-SAT-"+exp.expname, cfg_update={"time": "120", "mail-type": "FAIL,END", "ntasks": "1",
+ "ntasks-per-node": "1", "ntasks-per-core": "1"})
+ cmd = cluster.python_enstools+' /home/fs71386/lkugler/osse_analysis/plot_from_raw/analyze_fc.py '+exp.expname+' has_node sat verif1d FSS BS'
+ s.run(cmd, depends_on=[depends_on])
+
+def verify_wrf(depends_on=None):
+ s = my_Slurm("verif-WRF-"+exp.expname, cfg_update={"time": "180", "mail-type": "FAIL,END", "ntasks": "1",
"ntasks-per-node": "1", "ntasks-per-core": "1"})
- cmd = cluster.python_enstools+' /home/fs71386/lkugler/osse_analysis/plot_from_raw/analyze_fc.py '+exp.expname+' has_node sat wrf verif1d FSS BS'
+ cmd = cluster.python_enstools+' /home/fs71386/lkugler/osse_analysis/plot_from_raw/analyze_fc.py '+exp.expname+' has_node wrf verif1d FSS BS'
s.run(cmd, depends_on=[depends_on])
+def verify_fast(depends_on=None):
+ s = my_Slurm("verif-fast-"+exp.expname, cfg_update={"time": "30", "mail-type": "FAIL", "ntasks": "1",
+ "ntasks-per-node": "1", "ntasks-per-core": "1"})
+ cmd = cluster.python_enstools+' /home/fs71386/lkugler/osse_analysis/plot_fast/plot_single_exp.py '+exp.expname
+ s.run(cmd, depends_on=[depends_on])
################################
if __name__ == "__main__":
@@ -237,21 +248,37 @@ if __name__ == "__main__":
backup_scripts()
id = None
- init_time = dt.datetime(2008, 7, 30, 12)
- time = dt.datetime(2008, 7, 30, 12, 30)
+ if False: # warm bubble
+ prior_path_exp = '/gpfs/data/fs71386/lkugler/sim_archive/exp_v1.19_P3_wbub7_noDA'
+
+ init_time = dt.datetime(2008, 7, 30, 12)
+ time = dt.datetime(2008, 7, 30, 12,30)
+ last_assim_time = dt.datetime(2008, 7, 30, 13,30)
+ forecast_until = dt.datetime(2008, 7, 30, 18)
+
+ prepare_WRFrundir(init_time)
+ id = run_ideal(depends_on=id)
+ id = wrfinput_insert_wbubble(depends_on=id)
+
+ if True: # random
+ prior_path_exp = '/gpfs/data/fs71386/lkugler/sim_archive/exp_v1.19_P2_noDA'
+
+ init_time = dt.datetime(2008, 7, 30, 12)
+ time = dt.datetime(2008, 7, 30, 13)
+ last_assim_time = dt.datetime(2008, 7, 30, 14)
+ forecast_until = dt.datetime(2008, 7, 30, 18)
+
+ prepare_WRFrundir(init_time)
+ id = run_ideal(depends_on=id)
- id = prepare_WRFrundir(init_time)
- #id = run_ideal(depends_on=id)
#prior_path_exp = cluster.archivedir #
- #prior_path_exp = '/gpfs/data/fs71386/lkugler/sim_archive/exp_v1.19_P2_noDA'
- prior_path_exp = '/gpfs/data/fs71386/lkugler/sim_archive/exp_v1.19_P3_wbub7_noDA'
- #id = wrfinput_insert_wbubble(depends_on=id)
+ #prior_path_exp = '/gpfs/data/fs71386/lkugler/sim_archive/exp_v1.19_P5+su_noDA'
prior_init_time = init_time
prior_valid_time = time
- while time <= dt.datetime(2008, 7, 30, 13, 30):
+ while time <= last_assim_time:
# usually we take the prior from the current time
# but one could use a prior from a different time from another run
@@ -269,9 +296,9 @@ if __name__ == "__main__":
# How long shall we integrate?
timedelta_integrate = timedelta_btw_assim
output_restart_interval = timedelta_btw_assim.total_seconds()/60
- if time == dt.datetime(2008, 7, 30, 13,30): #this_forecast_init.minute in [0,]: # longer forecast every full hour
- timedelta_integrate = dt.timedelta(hours=4.5)
- output_restart_interval = 9999
+ if time == last_assim_time: #this_forecast_init.minute in [0,]: # longer forecast every full hour
+ timedelta_integrate = forecast_until - last_assim_time # dt.timedelta(hours=4)
+ output_restart_interval = 9999 #timedelta_btw_assim.total_seconds()/60 # 9999
# 3) Run WRF ensemble
id = run_ENS(begin=time, # start integration from here
@@ -291,4 +318,6 @@ if __name__ == "__main__":
prior_init_time = time - timedelta_btw_assim
id = gen_obsseq(id)
- verify(id_sat)
+ verify_sat(id_sat)
+ verify_wrf(id)
+ verify_fast(id)
diff --git a/templates/input.eval.nml b/templates/input.eval.nml
index 92212cad31ee54584af062f76d80425f4c47fe1a..e9bfa9219a51991aa00c4ae67324df789859e183 100644
--- a/templates/input.eval.nml
+++ b/templates/input.eval.nml
@@ -139,7 +139,7 @@
'../../../observations/forward_operators/obs_def_gps_mod.f90',
'../../../observations/forward_operators/obs_def_vortex_mod.f90',
'../../../observations/forward_operators/obs_def_gts_mod.f90',
- '../../../observations/forward_operators/obs_def_rttov_mod.f90',
+ '../../../observations/forward_operators/obs_def_rttov13_mod.f90',
/
&obs_kind_nml
@@ -156,10 +156,10 @@
'ACARS_TEMPERATURE',
'GPSRO_REFRACTIVITY',
'DOPPLER_RADIAL_VELOCITY',
- 'RADAR_REFLECTIVITY',
+ 'RADAR_REFLECTIVITY',
'MSG_4_SEVIRI_RADIANCE',
- 'MSG_4_SEVIRI_TB',
- 'MSG_4_SEVIRI_BDRF',
+ 'MSG_4_SEVIRI_TB',
+ 'MSG_4_SEVIRI_BDRF',
'LAND_SFC_PRESSURE',
'SYNOP_SURFACE_PRESSURE',
'SYNOP_TEMPERATURE',
@@ -214,7 +214,7 @@
'V', 'QTY_V_WIND_COMPONENT', 'TYPE_V', 'UPDATE','999',
'W', 'QTY_VERTICAL_VELOCITY', 'TYPE_W', 'UPDATE','999',
'PH', 'QTY_GEOPOTENTIAL_HEIGHT', 'TYPE_GZ', 'UPDATE','999',
- 'T', 'QTY_POTENTIAL_TEMPERATURE','TYPE_T', 'UPDATE','999',
+ 'THM', 'QTY_POTENTIAL_TEMPERATURE','TYPE_T', 'UPDATE','999',
'MU', 'QTY_PRESSURE', 'TYPE_MU', 'UPDATE','999',
'QVAPOR','QTY_VAPOR_MIXING_RATIO', 'TYPE_QV', 'UPDATE','999',
'QICE', 'QTY_ICE_MIXING_RATIO', 'TYPE_QI', 'UPDATE','999',
diff --git a/templates/input.nml b/templates/input.nml
index 31e34c1292ce18ea647b5ef2156aae13c935cd9a..9f7de60eed02a08a6bd904dbf8aa952defd9c6c0 100644
--- a/templates/input.nml
+++ b/templates/input.nml
@@ -141,7 +141,7 @@
'../../../observations/forward_operators/obs_def_gps_mod.f90',
'../../../observations/forward_operators/obs_def_vortex_mod.f90',
'../../../observations/forward_operators/obs_def_gts_mod.f90',
- '../../../observations/forward_operators/obs_def_rttov_mod.f90',
+ '../../../observations/forward_operators/obs_def_rttov13_mod.f90',
/
&obs_kind_nml
@@ -216,7 +216,7 @@
'V', 'QTY_V_WIND_COMPONENT', 'TYPE_V', 'UPDATE','999',
'W', 'QTY_VERTICAL_VELOCITY', 'TYPE_W', 'UPDATE','999',
'PH', 'QTY_GEOPOTENTIAL_HEIGHT', 'TYPE_GZ', 'UPDATE','999',
- 'T', 'QTY_POTENTIAL_TEMPERATURE','TYPE_T', 'UPDATE','999',
+ 'THM', 'QTY_POTENTIAL_TEMPERATURE','TYPE_T', 'UPDATE','999',
'MU', 'QTY_PRESSURE', 'TYPE_MU', 'UPDATE','999',
'QVAPOR','QTY_VAPOR_MIXING_RATIO', 'TYPE_QV', 'UPDATE','999',
'QICE', 'QTY_ICE_MIXING_RATIO', 'TYPE_QI', 'UPDATE','999',
diff --git a/templates/input.prioronly.nml b/templates/input.prioronly.nml
index 994e6e33cf499d236705e3e2dee8c57dd7c17d59..f5d30f4ee4ef2083ecea61e9b301c4cb6356e61a 100644
--- a/templates/input.prioronly.nml
+++ b/templates/input.prioronly.nml
@@ -139,7 +139,7 @@
'../../../observations/forward_operators/obs_def_gps_mod.f90',
'../../../observations/forward_operators/obs_def_vortex_mod.f90',
'../../../observations/forward_operators/obs_def_gts_mod.f90',
- '../../../observations/forward_operators/obs_def_rttov_mod.f90',
+ '../../../observations/forward_operators/obs_def_rttov13_mod.f90',
/
&obs_kind_nml
@@ -210,7 +210,7 @@
'V', 'QTY_V_WIND_COMPONENT', 'TYPE_V', 'UPDATE','999',
'W', 'QTY_VERTICAL_VELOCITY', 'TYPE_W', 'UPDATE','999',
'PH', 'QTY_GEOPOTENTIAL_HEIGHT', 'TYPE_GZ', 'UPDATE','999',
- 'T', 'QTY_POTENTIAL_TEMPERATURE','TYPE_T', 'UPDATE','999',
+ 'THM', 'QTY_POTENTIAL_TEMPERATURE','TYPE_T', 'UPDATE','999',
'MU', 'QTY_PRESSURE', 'TYPE_MU', 'UPDATE','999',
'QVAPOR','QTY_VAPOR_MIXING_RATIO', 'TYPE_QV', 'UPDATE','999',
'QICE', 'QTY_ICE_MIXING_RATIO', 'TYPE_QI', 'UPDATE','999',
diff --git a/templates/obs_def_rttov.IR.nml b/templates/obs_def_rttov.IR.nml
index 26c0098988fa1d76b30bced0a1af7717f785746f..09deccdd87f37096ce78fe1a24d8b7969c8767ce 100644
--- a/templates/obs_def_rttov.IR.nml
+++ b/templates/obs_def_rttov.IR.nml
@@ -48,7 +48,7 @@
aerosl_type = 1
add_clouds = .true.
ice_scheme = 1
- use_icede = .false.
+ use_icede = .true.
idg_scheme = 2
user_aer_opt_param = .false.
user_cld_opt_param = .false.
diff --git a/templates/obs_def_rttov.VIS.nml b/templates/obs_def_rttov.VIS.nml
index 2b97f672ef18403cbcbd7bce6bd8120e06eed0f7..50099921558e3ecd438fac23c97c86b24d4cb3b9 100644
--- a/templates/obs_def_rttov.VIS.nml
+++ b/templates/obs_def_rttov.VIS.nml
@@ -48,7 +48,7 @@
aerosl_type = 1
add_clouds = .true.
ice_scheme = 1
- use_icede = .false.
+ use_icede = .true.
idg_scheme = 1
user_aer_opt_param = .false.
user_cld_opt_param = .false.
diff --git a/tests/test_dart-rttov.py b/tests/test_dart-rttov.py
new file mode 100644
index 0000000000000000000000000000000000000000..11ad6f349e249e44b247f7d417b2e012b6af423e
--- /dev/null
+++ b/tests/test_dart-rttov.py
@@ -0,0 +1,154 @@
+import os, filecmp, shutil
+import numpy as np
+import datetime as dt
+import pandas as pd
+
+from config.cfg import exp, cluster
+from dartwrf import obsseq
+import dartwrf.create_obsseq as osq
+import dartwrf.assim_synth_obs as aso
+from dartwrf import wrfout_add_geo
+
+import matplotlib as mpl
+import proplot as plt
+import xarray as xr
+
+n_obs = 22500
+vis = dict(plotname='VIS 0.6µm', plotunits='[1]',
+ kind='MSG_4_SEVIRI_BDRF', sat_channel=1, n_obs=n_obs,
+ error_generate=0, error_assimilate=0.05,
+ cov_loc_radius_km=20)
+wv73 = dict(plotname='Brightness temperature WV 7.3µm', plotunits='[K]',
+ kind='MSG_4_SEVIRI_TB', sat_channel=6, n_obs=n_obs,
+ error_generate=0, error_assimilate=False,
+ cov_loc_radius_km=20)
+
+def test_rttov():
+
+ times = pd.date_range(start=dt.datetime(2008, 7, 30, 13,30),
+ end=dt.datetime(2008, 7, 30, 15,30),
+ freq=dt.timedelta(minutes=15))
+
+ case = 'exp_v1.18_P1_nature'
+ # case = 'exp_v1.19_P3_wbub7_nat'
+ for obsvar in ['VIS06']: #'WV73',
+
+ for time in times:
+
+ fname = '/home/fs71386/lkugler/data/analysis/'+case+'/dart-rttov-compare_'+obsvar+time.strftime('_%H:%M')+'-2.png'
+
+ # if os.path.isfile(fname):
+ # print(fname, 'skipping...')
+ # continue
+
+ ds = xr.open_dataset('/home/fs71386/lkugler/data/sim_archive/'+case+'/2008-07-30_06:00/1/'+time.strftime('RT2_wrfout_d01_2008-07-30_%H:%M:00.nc'))
+ #
+ da = ds[obsvar].squeeze()
+
+
+ shutil.copy('/home/fs71386/lkugler/data/sim_archive/'+case+'/2008-07-30_06:00/1/'+time.strftime('wrfout_d01_2008-07-30_%H:%M:00'),
+ cluster.dartrundir + "/wrfout_d01")
+
+ wrfout_add_geo.run(
+ cluster.dartrundir + "/../geo_em.d01.nc", cluster.dartrundir + "/wrfout_d01"
+ )
+
+
+ error_generate = np.zeros(n_obs)
+ if obsvar == 'VIS06':
+ osq.create_obsseqin_alltypes(time, [vis])
+ else:
+ osq.create_obsseqin_alltypes(time, [wv73])
+ aso.set_DART_nml()
+ aso.run_perfect_model_obs()
+
+
+ obs = obsseq.ObsSeq(cluster.dartrundir+'/obs_seq.out')
+ obs_map = obs.df.truth.values.reshape(150,150)
+
+
+ fig, ax = plt.subplots(ncols=2)
+
+ if obsvar == 'VIS06':
+ norm = mpl.colors.Normalize(vmin=0, vmax=1)
+ levels = plt.arange(0, 1, .05)
+ else:
+ norm = mpl.colors.Normalize(vmin=210, vmax=260)
+ levels = plt.arange(200, 260, 5)
+
+ ax[0].pcolormesh(obs_map, norm=norm, levels=levels)
+ ax[0].set_title('DART-RTTOV')
+ m1 = ax[1].pcolormesh(da.values[24:175,24:175], norm=norm, levels=levels)
+ ax[1].set_title('py-RTTOV')
+ fig.colorbar(m1, loc='r')
+ fig.suptitle(time.strftime('%H:%M'))
+
+ fig.savefig(fname)
+ print(fname)
+
+
+
+ fig, ax = plt.subplots()
+
+ # norm = mpl.colors.Normalize(vmin=0, vmax=1)
+ # levels = plt.arange(0, 1, .05)
+ # norm = mpl.colors.Normalize(vmin=210, vmax=260)
+ # levels = plt.arange(200, 260, 5)
+
+ diff = da.values[25:175,25:175] - obs_map
+
+ ax[0].pcolormesh(diff, norm=norm, levels=levels)
+ ax[0].set_title('pyRTTOV - DART-RTTOV')
+ fig.colorbar(m1, loc='r')
+ fig.suptitle(time.strftime('%H:%M'))
+
+ fname = '/home/fs71386/lkugler/data/analysis/'+case+'/dart-rttov-compare_'+obsvar+time.strftime('_%H:%M')+'-diff.png'
+ fig.savefig(fname)
+ print(fname)
+
+ #from IPython import embed; embed()
+
+
+ # shutil.copy(cluster.dartrundir+'/obs_seq.out',
+ # '/home/fs71386/lkugler/data/analysis/'+case+time.strftime('/obs_seq_'+obsvar+'_%H:%M.out'))
+
+
+ # fig, ax = plt.subplots()
+ # m = ax.pcolormesh(obs_map)
+ # ax.colorbar(m, loc='r')
+ # ax.set_title(time.strftime('%H:%M'))
+ # fname = '/home/fs71386/lkugler/data/analysis/exp_v1.18_P1_nature/dart-rttov-output_'+time.strftime('%H:%M')+'.png'
+ # fig.savefig(fname)
+ # print(fname)
+
+
+if __name__ == '__main__':
+ test_rttov()
+ pass
+ # case = 'test_IR+VIS'
+ # time = dt.datetime(2008,7,30,13,30)
+
+
+
+ # error_generate = np.zeros(n_obs)
+ # osq.create_obsseqin_alltypes(time, [vis])
+ # aso.set_DART_nml()
+ # aso.run_perfect_model_obs()
+
+
+ # obs = obsseq.ObsSeq(cluster.dartrundir+'/obs_seq.out')
+ # obs_map = obs.df.truth.values.reshape(150,150)
+
+
+ # ds = xr.open_dataset('/home/fs71386/lkugler/data/run_DART/test_IR+VIS/RT_wrfout_d01_cloudy.nc')
+ # da = ds['VIS06'].squeeze()
+ # pyrttov = da.values[25:175,25:175]
+
+ # diff = pyrttov - obs_map
+
+ # i, j = np.unravel_index(diff.argmax(), (150,150))
+ # print(i,j)
+
+ # print(pyrttov[i,j], obs_map[i,j])
+
+ # from IPython import embed; embed()
diff --git a/tests/test_obsseq.py b/tests/test_obsseq.py
index 516f31400e7f902408ee6a796bf7201a305c3aa9..96d7dbe013b9d5a15ca6bf6fbcf9e43ac5286a01 100644
--- a/tests/test_obsseq.py
+++ b/tests/test_obsseq.py
@@ -14,20 +14,44 @@ def test_oso():
obs = obsseq.ObsSeq(input)
# select a subset
- obs.superob(window_km=50)
+ obs.df = obs.df.superob(window_km=50)
+
+ #obs.plot(f_out="/home/fs71386/lkugler/data/analysis/map_obs_superobs_obs10-50.png")
# write to obs_seq.out in DART format
obs.to_dart(f=output)
- assert filecmp.cmp(output, output_true)
+ #assert filecmp.cmp(output, output_true)
+
+def test_plot():
+ input = cluster.scriptsdir + "/../tests/obs_seq.out_obs2"
+ obs = obsseq.ObsSeq(input)
+ # select a subset
+ obs.df = obs.df.superob(window_km=10)
+
+ obs.plot(f_out="/home/fs71386/lkugler/data/analysis/map_obs_superobs_obs2-10.png")
def test_osf():
input = cluster.scriptsdir + "/../tests/obs_seq.final"
obs = obsseq.ObsSeq(input)
- prior_Hx = obs.get_prior_Hx_matrix()
+ prior_Hx = obs.df.get_prior_Hx()
# TODO: compare with given truth
+def test_superob():
+ input = cluster.scriptsdir + "/../tests/obs_seq.out_obs2"
+ output = cluster.scriptsdir + "/../tests/obs_seq.test.out"
+ obs = obsseq.ObsSeq(input)
+
+ obs.df = obs.df.superob(window_km=10)
+
+ obs.to_dart(f=output)
+
+ obs = obsseq.ObsSeq(output)
+
+ from IPython import embed; embed()
+
if __name__ == '__main__':
- pass
\ No newline at end of file
+ test_superob()
+ pass