snapshot...

modules/aeolus/aeolus_amv.py

@@ -671,7 +671,7 @@ 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):
+def analyze(prof_da, amvs_da, prof_locs_da, dist_threshold=5.0):
     # sort amvs by distance to profile
     dst = amvs_da.sel(num_params='dist_to_prof')
     s_i = np.argsort(dst, axis=1)
@@ -686,14 +686,11 @@ def analyze(prof_da, amvs_da, prof_locs_da):
     for k in range(prof_da.shape[0]):
         prof_da[k] = prof_da[k, s_i[k, ::-1]]
 
+    # analyze single cloud layer profiles
     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]
+    print('number of one layer profs: ', one_lyr_profs.shape[0])
 
     hgt_vld = one_lyr_amvs.sel(num_params='H_3D') > 0
     hgt_vld = hgt_vld.values
@@ -701,19 +698,64 @@ def analyze(prof_da, amvs_da, prof_locs_da):
 
     one_lyr_amvs = one_lyr_amvs[has_vld]
     one_lyr_profs = one_lyr_profs[has_vld]
+    print('number of one layer profs with at least one valid AMV height', one_lyr_profs.shape[0])
 
     # 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 = (np.logical_and(dst > 0.0, dst < dist_threshold)).values
     close = close.sum(1) > 0
 
     one_lyr_amvs = one_lyr_amvs[close]
     one_lyr_profs = one_lyr_profs[close]
+    print('number of one layer profs with at least one AMV within 5km: ', one_lyr_profs.shape[0])
 
+    cnt = 0
     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)
+        dst = one_lyr_amvs[k, :, ].sel(num_params='dist_to_prof')
+        b = np.logical_and(dst > 0.0, dst < dist_threshold)
+        h_3d = one_lyr_amvs[k, :].sel(num_params='H_3D')
+        h_3d = h_3d[b]
+        vld = h_3d > 0
+        h_3d = h_3d[vld]
+        if len(h_3d) > 0:
+            in_lyr = np.logical_and(h_3d > one_lyr_profs[k, 0, 0], h_3d < one_lyr_profs[k, 0, 1])
+            cnt += np.sum(in_lyr)
+
+    print(cnt)
+
+    # Do calculations for single and multi-layer combined
+    hgt_vld = amvs_da.sel(num_params='H_3D') > 0
+    hgt_vld = hgt_vld.values
+    has_vld = hgt_vld.sum(1) > 0
+
+    amvs_da = amvs_da[has_vld]
+    prof_da = prof_da[has_vld]
+    prof_locs_da = prof_locs_da[has_vld]
+    print('number of profs with at least one valid AMV height', prof_da.shape[0])
 
+    # compare profile highest cloud layer top to closest within 5km
+    dst = amvs_da.sel(num_params='dist_to_prof')
+    close = (np.logical_and(dst > 0.0, dst < dist_threshold)).values
+    close = close.sum(1) > 0
+
+    amvs_da = amvs_da[close]
+    prof_da = prof_da[close]
+    prof_locs_da = prof_locs_da[close]
+    print('number of profs with at least one AMV within 5km: ', prof_da.shape[0])
+
+    cnt = 0
+    for k in range(prof_da.shape[0]):
+        dst = amvs_da[k, :, ].sel(num_params='dist_to_prof')
+        b = np.logical_and(dst > 0.0, dst < dist_threshold)
+        h_3d = amvs_da[k, :].sel(num_params='H_3D')
+        h_3d = h_3d[b]
+        vld = h_3d > 0
+        h_3d = h_3d[vld]
+        if len(h_3d) > 0:
+            nlevs = prof_locs_da[k].values.astype(int)[3]
+            for j in range(nlevs):
+                in_lyr = np.logical_and(h_3d > prof_da[k, j, 0], h_3d < prof_da[k, j, 1])
+                cnt += np.sum(in_lyr)
+
+    print(cnt)
     return one_lyr_profs, one_lyr_amvs