snapshot...

modules/aeolus/aeolus_amv.py

@@ -634,8 +634,8 @@ def subset_by_time(match_file, dt_str_0, dt_str_1, amv_var_names):
 
     amvs = np.zeros((n_times, mx_namvs, nvars))
     profs = np.zeros((n_times, mx_nlevs, npvars))
-    amvs.fill(np.nan)
-    profs.fill(np.nan)
+    amvs.fill(-999.0)
+    profs.fill(-999.0)
 
     accum_prf = prf_idx_start
     accum_amv = amv_idx_start
@@ -671,3 +671,49 @@ def subset_by_time(match_file, dt_str_0, dt_str_1, amv_var_names):
     return prof_da, amvs_da, prof_locs_da
 
 
+def analyze(prof_da, amvs_da, prof_locs_da):
+    # sort amvs by distance to profile
+    dst = amvs_da.sel(num_params='dist_to_prof')
+    s_i = np.argsort(dst, axis=1)
+    s_i = s_i.values
+    for k in range(amvs_da.shape[0]):
+        amvs_da[k] = amvs_da[k, s_i[k]]
+
+    # sort profiles by level highest to lowest
+    top = prof_da.sel(num_params='layer_top')
+    s_i = np.argsort(top, axis=1)
+    s_i = s_i.values
+    for k in range(prof_da.shape[0]):
+        prof_da[k] = prof_da[k, s_i[k, ::-1]]
+
+    one_lyr = prof_locs_da.sel(space='num_levels') == 1
+    dst = amvs_da.sel(num_params='dist_to_prof')
+    close = np.logical_and(dst > 0.0, dst < 5.0)
+
+    # compare for one cloud layer profile
+
+    one_lyr_profs = prof_da[one_lyr]
+    one_lyr_amvs = amvs_da[one_lyr]
+
+    hgt_vld = one_lyr_amvs.sel(num_params='H_3D') > 0
+    hgt_vld = hgt_vld.values
+    has_vld = hgt_vld.sum(1) > 0
+
+    one_lyr_amvs = one_lyr_amvs[has_vld]
+    one_lyr_profs = one_lyr_profs[has_vld]
+
+    # compare profile highest cloud layer top to closest within 5km
+    dst = one_lyr_amvs.sel(num_params='dist_to_prof')
+    close = (np.logical_and(dst > 0.0, dst < 2.0)).values
+    close = close.sum(1) > 0
+
+    one_lyr_amvs = one_lyr_amvs[close]
+    one_lyr_profs = one_lyr_profs[close]
+
+    for k in range(one_lyr_profs.shape[0]):
+        dst = one_lyr_amvs[k, :, 2]
+        b = np.logical_and(dst > 0.0, dst < 2.0)
+        diffs = one_lyr_amvs[k, b, 0] - one_lyr_profs[k, 0, 1]
+        print(diffs.values, one_lyr_amvs[k, b, 2].values)
+
+    return one_lyr_profs, one_lyr_amvs