docs/source/tracking_template_object_compare.rst

+Template Object Compare
+=======================
+
+Template for tracking, where the tracking strategy is solely based on pairwise comparison of object descriptions from consecutive timesteps.
+
+
+.. autoclass:: enstools.feature.tracking.template_object_compare.TrackingCompareTemplate
+

enstools/feature/identification/__init__.py

 # init file for identification
-from .identification import IdentificationTechnique
\ No newline at end of file
+from .identification import IdentificationStrategy
\ No newline at end of file

enstools/feature/identification/_proto_gen/storm_pb2.py

+# -*- coding: utf-8 -*-
+# Generated by the protocol buffer compiler.  DO NOT EDIT!
+# source: tmpz751s9fm
+"""Generated protocol buffer code."""
+from google.protobuf.internal import builder as _builder
+from google.protobuf import descriptor as _descriptor
+from google.protobuf import descriptor_pool as _descriptor_pool
+from google.protobuf import symbol_database as _symbol_database
+# @@protoc_insertion_point(imports)
+
+_sym_db = _symbol_database.Default()
+
+
+
+
+DESCRIPTOR = _descriptor_pool.Default().AddSerializedFile(b'\n\x0btmpz751s9fm\"\x1f\n\x03Pos\x12\x0b\n\x03lat\x18\x01 \x02(\x02\x12\x0b\n\x03lon\x18\x02 \x02(\x02\"B\n\nProperties\x12\x14\n\x0cmin_pressure\x18\x01 \x02(\x02\x12\x1e\n\x10min_pressure_pos\x18\x02 \x02(\x0b\x32\x04.Pos\")\n\tVoxelData\x12\r\n\x05index\x18\x01 \x03(\r\x12\r\n\x05value\x18\x02 \x01(\x02\"-\n\nVertexData\x12\x10\n\x08position\x18\x01 \x03(\x02\x12\r\n\x05value\x18\x02 \x01(\x02\"+\n\x08\x46\x61\x63\x65\x44\x61ta\x12\x10\n\x08vertices\x18\x01 \x03(\r\x12\r\n\x05value\x18\x02 \x01(\x02\"C\n\x13VoxelRepresentation\x12\x0c\n\x04\x64\x65sc\x18\x01 \x01(\t\x12\x1e\n\nvoxel_data\x18\x02 \x03(\x0b\x32\n.VoxelData\"D\n\x12LineRepresentation\x12\x0c\n\x04\x64\x65sc\x18\x01 \x01(\t\x12 \n\x0bvertex_data\x18\x02 \x03(\x0b\x32\x0b.VertexData\"f\n\x16\x42oundaryRepresentation\x12\x0c\n\x04\x64\x65sc\x18\x01 \x01(\t\x12 \n\x0bvertex_data\x18\x02 \x03(\x0b\x32\x0b.VertexData\x12\x1c\n\tface_data\x18\x03 \x03(\x0b\x32\t.FaceData\"2\n\tGraphNode\x12\x0c\n\x04time\x18\x01 \x02(\t\x12\x17\n\x06object\x18\x02 \x02(\x0b\x32\x07.Object\"K\n\x0fGraphConnection\x12\x1a\n\x06parent\x18\x01 \x02(\x0b\x32\n.GraphNode\x12\x1c\n\x08\x63hildren\x18\x02 \x03(\x0b\x32\n.GraphNode\".\n\x0bObjectGraph\x12\x1f\n\x05\x65\x64ges\x18\x01 \x03(\x0b\x32\x10.GraphConnection\"\xb4\x01\n\x06Object\x12\n\n\x02id\x18\x01 \x02(\x05\x12%\n\x08line_rep\x18\x02 \x03(\x0b\x32\x13.LineRepresentation\x12\'\n\tvoxel_rep\x18\x03 \x01(\x0b\x32\x14.VoxelRepresentation\x12-\n\x0c\x62oundary_rep\x18\x04 \x03(\x0b\x32\x17.BoundaryRepresentation\x12\x1f\n\nproperties\x18\x05 \x01(\x0b\x32\x0b.Properties\"8\n\x08Timestep\x12\x12\n\nvalid_time\x18\x01 \x01(\t\x12\x18\n\x07objects\x18\x02 \x03(\x0b\x32\x07.Object\"\x99\x01\n\x0cTrackableSet\x12\x11\n\tinit_time\x18\x01 \x01(\t\x12\x0e\n\x06member\x18\x02 \x01(\r\x12\r\n\x05level\x18\x03 \x01(\x02\x12\x1c\n\ttimesteps\x18\x04 \x03(\x0b\x32\t.Timestep\x12\x1b\n\x05graph\x18\x05 \x01(\x0b\x32\x0c.ObjectGraph\x12\x1c\n\x06tracks\x18\x06 \x03(\x0b\x32\x0c.ObjectGraph\"_\n\x12\x44\x61tasetDescription\x12\x0c\n\x04name\x18\x01 \x01(\t\x12\x0c\n\x04\x66ile\x18\x02 \x01(\t\x12\x10\n\x08run_time\x18\x03 \x02(\t\x12\x1b\n\x04sets\x18\x04 \x03(\x0b\x32\r.TrackableSet')
+
+_builder.BuildMessageAndEnumDescriptors(DESCRIPTOR, globals())
+_builder.BuildTopDescriptorsAndMessages(DESCRIPTOR, 'tmpz751s9fm_pb2', globals())
+if _descriptor._USE_C_DESCRIPTORS == False:
+
+  DESCRIPTOR._options = None
+  _POS._serialized_start=15
+  _POS._serialized_end=46
+  _PROPERTIES._serialized_start=48
+  _PROPERTIES._serialized_end=114
+  _VOXELDATA._serialized_start=116
+  _VOXELDATA._serialized_end=157
+  _VERTEXDATA._serialized_start=159
+  _VERTEXDATA._serialized_end=204
+  _FACEDATA._serialized_start=206
+  _FACEDATA._serialized_end=249
+  _VOXELREPRESENTATION._serialized_start=251
+  _VOXELREPRESENTATION._serialized_end=318
+  _LINEREPRESENTATION._serialized_start=320
+  _LINEREPRESENTATION._serialized_end=388
+  _BOUNDARYREPRESENTATION._serialized_start=390
+  _BOUNDARYREPRESENTATION._serialized_end=492
+  _GRAPHNODE._serialized_start=494
+  _GRAPHNODE._serialized_end=544
+  _GRAPHCONNECTION._serialized_start=546
+  _GRAPHCONNECTION._serialized_end=621
+  _OBJECTGRAPH._serialized_start=623
+  _OBJECTGRAPH._serialized_end=669
+  _OBJECT._serialized_start=672
+  _OBJECT._serialized_end=852
+  _TIMESTEP._serialized_start=854
+  _TIMESTEP._serialized_end=910
+  _TRACKABLESET._serialized_start=913
+  _TRACKABLESET._serialized_end=1066
+  _DATASETDESCRIPTION._serialized_start=1068
+  _DATASETDESCRIPTION._serialized_end=1163
+# @@protoc_insertion_point(module_scope)

enstools/feature/identification/_proto_gen/threshold_pb2.py

 # -*- coding: utf-8 -*-
 # Generated by the protocol buffer compiler.  DO NOT EDIT!
+<<<<<<< HEAD
+# source: tmpwdhwhb1j
+=======
 # source: tmpncjfrx7c
+>>>>>>> origin/aew_kitweather
 """Generated protocol buffer code."""
+from google.protobuf.internal import builder as _builder
 from google.protobuf import descriptor as _descriptor
-from google.protobuf import message as _message
-from google.protobuf import reflection as _reflection
+from google.protobuf import descriptor_pool as _descriptor_pool
 from google.protobuf import symbol_database as _symbol_database
 # @@protoc_insertion_point(imports)
 
@@ -13,6 +17,9 @@ _sym_db = _symbol_database.Default()
 
 
 
+<<<<<<< HEAD
+DESCRIPTOR = _descriptor_pool.Default().AddSerializedFile(b'\n\x0btmpwdhwhb1j\"6\n\tMaskArray\x12\r\n\x05shape\x18\x01 \x03(\x05\x12\x0c\n\x04\x64\x61ta\x18\x02 \x02(\x0c\x12\x0c\n\x04zlib\x18\x03 \x02(\x08\"u\n\nProperties\x12\x18\n\x04mask\x18\x01 \x02(\x0b\x32\n.MaskArray\x12\x17\n\x0fouter_threshold\x18\x02 \x02(\x02\x12\x17\n\x0finner_threshold\x18\x03 \x02(\x02\x12\x1b\n\x13\x63omparison_operator\x18\x04 \x02(\t\")\n\tVoxelData\x12\r\n\x05index\x18\x01 \x03(\r\x12\r\n\x05value\x18\x02 \x01(\x02\"-\n\nVertexData\x12\x10\n\x08position\x18\x01 \x03(\x02\x12\r\n\x05value\x18\x02 \x01(\x02\"+\n\x08\x46\x61\x63\x65\x44\x61ta\x12\x10\n\x08vertices\x18\x01 \x03(\r\x12\r\n\x05value\x18\x02 \x01(\x02\"C\n\x13VoxelRepresentation\x12\x0c\n\x04\x64\x65sc\x18\x01 \x01(\t\x12\x1e\n\nvoxel_data\x18\x02 \x03(\x0b\x32\n.VoxelData\"D\n\x12LineRepresentation\x12\x0c\n\x04\x64\x65sc\x18\x01 \x01(\t\x12 \n\x0bvertex_data\x18\x02 \x03(\x0b\x32\x0b.VertexData\"f\n\x16\x42oundaryRepresentation\x12\x0c\n\x04\x64\x65sc\x18\x01 \x01(\t\x12 \n\x0bvertex_data\x18\x02 \x03(\x0b\x32\x0b.VertexData\x12\x1c\n\tface_data\x18\x03 \x03(\x0b\x32\t.FaceData\"2\n\tGraphNode\x12\x0c\n\x04time\x18\x01 \x02(\t\x12\x17\n\x06object\x18\x02 \x02(\x0b\x32\x07.Object\"K\n\x0fGraphConnection\x12\x1a\n\x06parent\x18\x01 \x02(\x0b\x32\n.GraphNode\x12\x1c\n\x08\x63hildren\x18\x02 \x03(\x0b\x32\n.GraphNode\".\n\x0bObjectGraph\x12\x1f\n\x05\x65\x64ges\x18\x01 \x03(\x0b\x32\x10.GraphConnection\"\xb4\x01\n\x06Object\x12\n\n\x02id\x18\x01 \x02(\x05\x12%\n\x08line_rep\x18\x02 \x03(\x0b\x32\x13.LineRepresentation\x12\'\n\tvoxel_rep\x18\x03 \x01(\x0b\x32\x14.VoxelRepresentation\x12-\n\x0c\x62oundary_rep\x18\x04 \x03(\x0b\x32\x17.BoundaryRepresentation\x12\x1f\n\nproperties\x18\x05 \x01(\x0b\x32\x0b.Properties\"8\n\x08Timestep\x12\x12\n\nvalid_time\x18\x01 \x01(\t\x12\x18\n\x07objects\x18\x02 \x03(\x0b\x32\x07.Object\"\x99\x01\n\x0cTrackableSet\x12\x11\n\tinit_time\x18\x01 \x01(\t\x12\x0e\n\x06member\x18\x02 \x01(\r\x12\r\n\x05level\x18\x03 \x01(\x02\x12\x1c\n\ttimesteps\x18\x04 \x03(\x0b\x32\t.Timestep\x12\x1b\n\x05graph\x18\x05 \x01(\x0b\x32\x0c.ObjectGraph\x12\x1c\n\x06tracks\x18\x06 \x03(\x0b\x32\x0c.ObjectGraph\"_\n\x12\x44\x61tasetDescription\x12\x0c\n\x04name\x18\x01 \x01(\t\x12\x0c\n\x04\x66ile\x18\x02 \x01(\t\x12\x10\n\x08run_time\x18\x03 \x02(\t\x12\x1b\n\x04sets\x18\x04 \x03(\x0b\x32\r.TrackableSet')
+=======
 DESCRIPTOR = _descriptor.FileDescriptor(
   name='tmpncjfrx7c',
   package='',
@@ -832,6 +839,41 @@ DatasetDescription = _reflection.GeneratedProtocolMessageType('DatasetDescriptio
   # @@protoc_insertion_point(class_scope:DatasetDescription)
   })
 _sym_db.RegisterMessage(DatasetDescription)
-
-
+>>>>>>> origin/aew_kitweather
+
+_builder.BuildMessageAndEnumDescriptors(DESCRIPTOR, globals())
+_builder.BuildTopDescriptorsAndMessages(DESCRIPTOR, 'tmpwdhwhb1j_pb2', globals())
+if _descriptor._USE_C_DESCRIPTORS == False:
+
+  DESCRIPTOR._options = None
+  _MASKARRAY._serialized_start=15
+  _MASKARRAY._serialized_end=69
+  _PROPERTIES._serialized_start=71
+  _PROPERTIES._serialized_end=188
+  _VOXELDATA._serialized_start=190
+  _VOXELDATA._serialized_end=231
+  _VERTEXDATA._serialized_start=233
+  _VERTEXDATA._serialized_end=278
+  _FACEDATA._serialized_start=280
+  _FACEDATA._serialized_end=323
+  _VOXELREPRESENTATION._serialized_start=325
+  _VOXELREPRESENTATION._serialized_end=392
+  _LINEREPRESENTATION._serialized_start=394
+  _LINEREPRESENTATION._serialized_end=462
+  _BOUNDARYREPRESENTATION._serialized_start=464
+  _BOUNDARYREPRESENTATION._serialized_end=566
+  _GRAPHNODE._serialized_start=568
+  _GRAPHNODE._serialized_end=618
+  _GRAPHCONNECTION._serialized_start=620
+  _GRAPHCONNECTION._serialized_end=695
+  _OBJECTGRAPH._serialized_start=697
+  _OBJECTGRAPH._serialized_end=743
+  _OBJECT._serialized_start=746
+  _OBJECT._serialized_end=926
+  _TIMESTEP._serialized_start=928
+  _TIMESTEP._serialized_end=984
+  _TRACKABLESET._serialized_start=987
+  _TRACKABLESET._serialized_end=1140
+  _DATASETDESCRIPTION._serialized_start=1142
+  _DATASETDESCRIPTION._serialized_end=1237
 # @@protoc_insertion_point(module_scope)

enstools/feature/identification/african_easterly_waves/configuration.py

@@ -8,30 +8,42 @@ import numpy as np
 # latS = -35
 # lonW = -100
 # lonE = 45
-data_lat = (-35, 35)
+data_lat = (-15, 35)
 data_lon = (-100, 45)
 
-aew_clim_dir = '/lsdf/MOD/Gruppe_Transregio/Gruppe_Knippertz/kitweather/data/era5/cv_clim_era5.nc' # TODO # '/home/christoph/phd/data/framework_example_ds/aew/' # '/project/meteo/w2w/C3/fischer/belanger/aew_clim/' # 
+aew_clim_dir = '/lsdf/MOD/Gruppe_Transregio/Gruppe_Knippertz/kitweather/data/era5/cv_clim_era5.nc' # 'C:\\Users\\Christoph\\phd\\data\\enstools-feature\\cv_clim_era5.nc' # '/home/christoph/phd/data/aew/clim/cv_clim_era5.nc' # '/home/christoph/phd/data/framework_example_ds/aew/' # '/project/meteo/w2w/C3/fischer/belanger/aew_clim/' #
+in_files = '/home/ws/he7273/phd_all/data/coll_oper/jja2021/jja2021.nc' # 'C:\\Users\\Christoph\\phd\\data\\enstools-feature\\2008_sum_uv.nc' # '/home/christoph/phd/data/framework_example_ds/aew/cv_aug_08.nc'
+
+out_dir = join('/home/ws/he7273/phd_all/data/aew/out/') # '/project/meteo/w2w/C3/fischer/belanger/out/'
+plot_dir = join('/home/ws/he7273/phd_all/data/aew/plots/') # '/project/meteo/w2w/C3/fischer/belanger/plots/'
+
 timedelta_ana = np.timedelta64(7, 'D')
 
 # aew_kitweather_ecmwf_dir = '/lsdf/MOD/Gruppe_Transregio/Gruppe_Knippertz/kitweather/data/ecmwf-hres/forecasts/'
 # num_files_ecmwf = 41 # files in finished forecast data
 # file_prefix_ecmwf = 'ecmwf-hres_latlon_1.0deg_'
 # file_prefix_ecmwf_rain = 'ecmwf-hres_latlon_0.4deg_'
-plot_dir_ecmwf = '/home/iconeps/plots/aew/ecmwf-hres/' # '/project/meteo/w2w/C3/fischer/belanger/plots/'
 
 # aew_kitweather_icon_dir = '/lsdf/MOD/Gruppe_Transregio/Gruppe_Knippertz/kitweather/data/icon-global-det/forecasts/'
 # num_files_icon = 31 # files in finished forecast data
 # file_prefix_icon = 'icon-global-det_latlon_1.0deg_'
 # file_prefix_icon_rain = 'icon-global-det_latlon_0.25deg_'
+plot_dir_ecmwf = '/home/iconeps/plots/aew/ecmwf-hres/' # '/project/meteo/w2w/C3/fischer/belanger/plots/'
+plot_dir_ecmwf_prefix = 'aew_ecmwf-hres_'
+plot_dir_ecmwf_ens = '/home/iconeps/plots/aew/ecmwf-ens/' # '/project/meteo/w2w/C3/fischer/belanger/plots/'
+plot_dir_ecmwf_ens_prefix = 'aew_ecmwf-ens_'
 plot_dir_icon = '/home/iconeps/plots/aew/icon-global-det/' # '/project/meteo/w2w/C3/fischer/belanger/plots/'
+plot_dir_icon_prefix = 'aew_icon-global-det_'
 plot_dir_gfs = '/home/iconeps/plots/aew/gfs/'
+plot_dir_gfs_prefix = 'aew_gfs_'
+plot_dir_gfs_ens = '/home/iconeps/plots/aew/gefs/' # '/project/meteo/w2w/C3/fischer/belanger/plots/'
+plot_dir_gfs_ens_prefix = 'aew_gefs_'
+
+ens_rain_threshold = 1.0 # display prob of exceeding X mm hr-1
 
 latest_run_dir = '/home/iconeps/plots/aew/'
 latest_run_info_file_name = 'latest_aew_run.txt'
 
-in_files = '/home/christoph/phd/data/aew/cv/2008cv.nc' # '/home/christoph/phd/data/framework_example_ds/aew/cv_aug_08.nc'
-
 # out_dir = join(expanduser("~") + '/phd/data/aew/out/') # '/project/meteo/w2w/C3/fischer/belanger/out/'
 
 cv_data_dir = '/project/meteo/w2w/C2/athul/data_athul/AEW_cv_data/' # reference where ALL cv data is (only for clim calc.)
@@ -49,16 +61,18 @@ def get_clim_file():
 
 
 # wave area to be extracted: at least one point of trough needs to be in this range
-wave_filter_lat = (0, 30)
+wave_filter_lat = (-5, 30) # TODO see bottom method -> more adaptive on lons.
 wave_filter_lon = (-110, 55)
 
-levels = [70000]  # 700 hPa
+levels = [70000.0]  # 700 hPa
+# u_dim = 'u700_rea'
+# v_dim = 'v700_rea'
 
 # time of interest, if None all
 # june-oct is AEW season
-start_date = None # '2008-09-01T00:00' #  # '2008-08-01T00:00'
-end_date = None # '2008-09-30T00:00' # None # '2008-08-03T00:00'
 
+start_date = None # '2008-08-01T00:00' #  # '2008-08-01T00:00'
+end_date = None # '2008-08-15T00:00' # None # '2008-08-03T00:00'
 
 # Algorithm parameters
 # max u wind (m/s) (0 = only keep west-propagating). Belanger: 2.5; Berry: 0.0
@@ -91,3 +105,17 @@ avg_speed_min_deg_per_h = avg_speed_min_m_per_s * 3.6 * 0.00914
 # at 10°N we have in longitude direction 0.00914 degrees/km (360/(40,000*cos(10deg)))
 speed_range_m_per_s = [3.0, 12.0] # [5,10], but be more gentle with polygons.
 speed_deg_per_h = [-m_per_s * 3.6 * 0.00914 for m_per_s in speed_range_m_per_s]  # negative -> westward [-5, -10]
+
+def is_point_in_area(lon, lat):
+
+    if lon < -65:
+        if 5 < lat < 28:
+            return True
+    elif lon < -50:
+        if 3 < lat < 27:
+            return True
+    else:
+        if 0 < lat < 26:
+            return True
+    
+    return False
\ No newline at end of file

enstools/feature/identification/african_easterly_waves/filtering.py

 
-# filter: keep current wavetrough if:
-#   troughs are within a certain spatial window
-#   length of the trough < threshold
-def keep_wavetrough(properties, config):
+"""
 
-    in_area = False
-    for line_pt in properties.line_pts:
 
-        # check if any point is outside filtering area
-        if (line_pt.lon > config.wave_filter_lon[0] and line_pt.lon < config.wave_filter_lon[1]
-                and line_pt.lat > config.wave_filter_lat[0] and line_pt.lat < config.wave_filter_lat[1]):
-            in_area = True
-
-    if not in_area: # no point of line segment in our area
-        return False
-
-    if properties.length_deg <= config.degree_len_thr: # too small
-        return False
-
-    return True
+"""
\ No newline at end of file

enstools/feature/identification/african_easterly_waves/identification.py

-from enstools.feature.identification import IdentificationTechnique
+from enstools.feature.identification import IdentificationStrategy
 import xarray as xr
 import numpy as np
 import os, sys
@@ -7,14 +7,15 @@ import metpy.calc as mpcalc
 from .util import calc_adv
 from matplotlib import pyplot as plt
 import cartopy.crs as ccrs
-from .filtering import keep_wavetrough
+
 from .processing import populate_object, compute_cv
 from skimage.draw import line_aa
-from enstools.feature.util.enstools_utils import get_u_var, get_v_var, get_vertical_dim, get_longitude_dim, get_latitude_dim
+from enstools.feature.util.enstools_utils import get_u_var, get_v_var, get_vertical_dim, get_longitude_dim, get_latitude_dim, get_init_time_dim, get_valid_time_dim
 import threading
 from skimage.draw import line
 
-class AEWIdentification(IdentificationTechnique):
+
+class AEWIdentification(IdentificationStrategy):
 
     def __init__(self, wt_out_file=False, wt_traj_dir=None, cv='cv', year_summer=None, month=None, **kwargs):
         """
@@ -32,7 +33,7 @@ class AEWIdentification(IdentificationTechnique):
         self.config = cfg  # config
         self.config.out_traj_dir = wt_traj_dir
         self.config.cv_name = cv
-        
+
         if year_summer is not None:
             if month is not None:
                 m_str = str(month).zfill(2)
@@ -41,7 +42,7 @@ class AEWIdentification(IdentificationTechnique):
             else:
                 self.config.start_date = str(year_summer) + '-06-01T00:00'
                 self.config.end_date = str(year_summer) + '-10-31T00:00'
-        
+
         self.config.out_wt = wt_out_file
         if wt_out_file:
             self.config.sum_over_all = True
@@ -49,7 +50,7 @@ class AEWIdentification(IdentificationTechnique):
         pass
 
     def precompute(self, dataset: xr.Dataset, **kwargs):
-        print("Precompute for PV identification...")
+        print("Precompute for AEW identification...")
 
         plt.switch_backend('agg')  # this is thread safe matplotlib but cant display.
 
@@ -58,13 +59,28 @@ class AEWIdentification(IdentificationTechnique):
         lat_range = self.config.data_lat
         clim_file = self.config.get_clim_file()
 
+        u_name = self.config.u_dim if hasattr(self.config, 'u_dim') and (self.config.u_dim is not None) else get_u_var(dataset)
+        v_name = self.config.v_dim if hasattr(self.config, 'v_dim') and (self.config.v_dim is not None) else get_v_var(dataset)
+
+        if u_name is None or v_name is None:
+            print("Could not locate u and v fields in dataset. Needed to compute advection terms.")
+            exit()
+
+        # restrict dimensions only to the ones present in u/v
+        all_dims = [d for d in dataset.dims]
+        for d in all_dims:
+            if d not in dataset[u_name].dims:
+                dataset = dataset.drop_dims(d)
+
         level_str = get_vertical_dim(dataset)
         lat_str = get_latitude_dim(dataset)
         lon_str = get_longitude_dim(dataset)
+        init_time_str = get_init_time_dim(dataset)
+        valid_time_str = get_valid_time_dim(dataset)
 
         if os.path.isfile(clim_file):
             cv_clim = xr.open_dataset(clim_file)
-            
+
         else:
             # generate: need all 40y of CV data.
             print("Climatology file not found. Computing climatology...")
@@ -72,31 +88,48 @@ class AEWIdentification(IdentificationTechnique):
             from .climatology import compute_climatology
             cv_clim = compute_climatology(self.config)
             cv_clim.to_netcdf(clim_file)
-        
+
+        lat_str_clim = get_latitude_dim(cv_clim)
+        lon_str_clim = get_longitude_dim(cv_clim)
         cv_clim = cv_clim.sel(
-            **{lat_str: slice(lat_range[0], lat_range[1])},
-            **{lon_str: slice(lon_range[0], lon_range[1])})
-        
+            **{lat_str_clim: slice(lat_range[0], lat_range[1])},
+            **{lon_str_clim: slice(lon_range[0], lon_range[1])})
+
         # --------------- SUBSET DATA ACCORDING TO CFG
         start_date_dt = np.datetime64(self.config.start_date) if self.config.start_date is not None else None
         end_date_dt = np.datetime64(self.config.end_date) if self.config.end_date is not None else None
-        
+
         # if data is lon=0..360, change it to -180..180
         dataset.coords[lon_str] = (dataset.coords[lon_str] + 180) % 360 - 180
         dataset = dataset.sortby(dataset[lon_str])
 
+        filter_time_str = init_time_str if init_time_str is not None else valid_time_str
+
+
         # get the data we want to investigate
         dataset = dataset.sortby(lat_str) # in case of descending
         dataset = dataset.sel(
             **{lat_str: slice(lat_range[0], lat_range[1])},
             **{lon_str: slice(lon_range[0], lon_range[1])},
-            time=slice(start_date_dt, end_date_dt))
-        
+            **{filter_time_str: slice(start_date_dt, end_date_dt)})
+
+        if len(dataset.time.values) == 0:
+            print("Given start and end time leads to no data to process.")
+            exit(1)
+
         # dataset = dataset.expand_dims('level')
         # level_str = 'level'
-        if level_str is not None:
-            dataset = dataset.sel(**{level_str: self.config.levels})
-
+        if level_str in dataset[u_name].dims:
+            dataset = dataset.sel(**{level_str: self.config.levels}) # 3-D wind field, select levels
+        elif level_str is not None:
+            dataset = dataset.sel(**{level_str: self.config.levels}) # 2-D, reduce 3-D field too
+            dataset[u_name] = dataset[u_name].expand_dims('level')
+            dataset[v_name] = dataset[v_name].expand_dims('level')
+        else:
+            dataset[u_name] = dataset[u_name].expand_dims(level=self.config.levels)
+            dataset[v_name] = dataset[v_name].expand_dims(level=self.config.levels)
+            level_str = 'level'
+        
         # rename cv_clim dimensions to be same as in data.
         cv_clim = cv_clim.rename({'lat': lat_str, 'lon': lon_str})
         if 'plev' in cv_clim.dims and 'plev' != level_str:
@@ -104,17 +137,16 @@ class AEWIdentification(IdentificationTechnique):
             cv_clim = cv_clim.rename({'plev': level_str})
             cv_clim = cv_clim.assign_coords({level_str: cv_clim[level_str] / 100})
 
+        # also only use levels also present in data
+        cv_clim = cv_clim.sel({level_str: dataset[level_str].values})
+
         if self.config.cv_name not in dataset.data_vars:
             print("Curvature Vorticity not found, trying to compute it out of U and V...")
-            u_name = get_u_var(dataset)
-            v_name = get_v_var(dataset)
-
             dataset = compute_cv(dataset, u_name, v_name, self.config.cv_name)
 
         # make dataset to 2.5 (or same as cv_clim)
         dataset = dataset.interp({lat_str: cv_clim.coords[lat_str], lon_str: cv_clim.coords[lon_str]})
 
-
         # make sure that lat and lon are last two dimensions
         if lat_str not in dataset[self.config.cv_name].coords.dims[-2:] or lon_str not in dataset[
                                                                                               self.config.cv_name].coords.dims[
@@ -124,8 +156,8 @@ class AEWIdentification(IdentificationTechnique):
 
         # --------------- DO NUMPY PARALLELIZED STUFF: CREATE TROUGH MASKS
 
-        u = dataset.u if 'u' in dataset.data_vars else dataset.U
-        v = dataset.v if 'v' in dataset.data_vars else dataset.V
+        u = dataset[u_name]
+        v = dataset[v_name]
         cv = dataset[self.config.cv_name]
         # smooth CV with kernel
 
@@ -190,6 +222,8 @@ class AEWIdentification(IdentificationTechnique):
     def identify(self, data_chunk: xr.Dataset, **kwargs):
         objs = []
         trough_mask_cur = data_chunk.trough_mask
+        if np.isnan(trough_mask_cur).all():
+            return data_chunk, objs
 
         def clip(tup, mint, maxt):
             return np.clip(tup, mint, maxt)
@@ -212,36 +246,41 @@ class AEWIdentification(IdentificationTechnique):
 
         id_ = 1
         for path in paths:
+
             # get new object, set id
             o = self.get_new_object()
             o.id = id_
             # populate it
-            populate_object(o.properties, path)
+            populate_object(o.properties, path, self.config)
 
             # add to objects if keep
-            if keep_wavetrough(o.properties, self.config):
-                objs.append(o)
-                id_ += 1
-                if not self.config.out_wt:
-                    continue
+            if not self.keep_wavetrough(o.properties):
+                continue
+
+            objs.append(o)
+            id_ += 1
 
-                for v_idx in range(len(path.vertices) - 1):
-                    start_lonlat = path.vertices[v_idx][0], path.vertices[v_idx][1]
-                    end_lonlat = path.vertices[v_idx + 1][0], path.vertices[v_idx + 1][1]
+            # if wavetrough out dataset, gen lines
+            if not self.config.out_wt:
+                continue
 
-                    start_idx = ((start_lonlat[0] - min_lon) / (max_lon - min_lon) * lons,
-                                 (start_lonlat[1] - min_lat) / (max_lat - min_lat) * lats)
-                    # start_idx = clip(start_idx, (0, 0), (lons, lats))
+            for v_idx in range(len(path.vertices) - 1):
+                start_lonlat = path.vertices[v_idx][0], path.vertices[v_idx][1]
+                end_lonlat = path.vertices[v_idx + 1][0], path.vertices[v_idx + 1][1]
 
-                    end_idx = ((end_lonlat[0] - min_lon) / (max_lon - min_lon) * lons,
-                               (end_lonlat[1] - min_lat) / (max_lat - min_lat) * lats)
-                    # end_idx = clip(end_idx, (0, 0), (lons, lats))
+                start_idx = ((start_lonlat[0] - min_lon) / (max_lon - min_lon) * lons,
+                             (start_lonlat[1] - min_lat) / (max_lat - min_lat) * lats)
+                # start_idx = clip(start_idx, (0, 0), (lons, lats))
 
-                    rr, cc, val = line_aa(int(start_idx[0]), int(start_idx[1]), int(end_idx[0]), int(end_idx[1]))
-                    rr = clip(rr, 0, lons - 1)
-                    cc = clip(cc, 0, lats - 1)
+                end_idx = ((end_lonlat[0] - min_lon) / (max_lon - min_lon) * lons,
+                           (end_lonlat[1] - min_lat) / (max_lat - min_lat) * lats)
+                # end_idx = clip(end_idx, (0, 0), (lons, lats))
 
-                    wt.data[cc, rr] = np.where(np.greater(val, wt.data[cc, rr]), val, wt.data[cc, rr])
+                rr, cc, val = line_aa(int(start_idx[0]), int(start_idx[1]), int(end_idx[0]), int(end_idx[1]))
+                rr = clip(rr, 0, lons - 1)
+                cc = clip(cc, 0, lats - 1)
+
+                wt.data[cc, rr] = np.where(np.greater(val, wt.data[cc, rr]), val, wt.data[cc, rr])
 
         return data_chunk, objs
 
@@ -267,7 +306,7 @@ class AEWIdentification(IdentificationTechnique):
             level_str = get_vertical_dim(dataset)
             if level_str is not None:
                 dataset = dataset.squeeze(drop=True)
-                
+
             if self.config.sum_over_all:
                 dataset['wavetroughs'] = dataset.wavetroughs.sum(dim='time')
 
@@ -276,7 +315,7 @@ class AEWIdentification(IdentificationTechnique):
             if not os.path.exists(self.config.out_traj_dir):
                 os.makedirs(self.config.out_traj_dir)
 
-            assert(len(data_desc.sets) == 1) # TODO assert one set. maybe expand at some point
+            assert (len(data_desc.sets) == 1)  # TODO assert one set. maybe expand at some point
             desc_set = data_desc.sets[0]
             desc_times = desc_set.timesteps
 
@@ -288,27 +327,29 @@ class AEWIdentification(IdentificationTechnique):
                 lon_list = []
                 lat_list = []
                 pres_list = []
-                max_pts_in_wt = -1 # TODO what if no wts
-                for o in ts.objects: # get lons and lats
+                max_pts_in_wt = -1  # TODO what if no wts
+                for o in ts.objects:  # get lons and lats
                     pt_list = o.properties.line_pts
                     lon_list.append(np.array([pt.lon for pt in pt_list]))
                     lat_list.append(np.array([pt.lat for pt in pt_list]))
                     pres_list.append(np.array([850.0 for pt in pt_list]))
                     max_pts_in_wt = max(max_pts_in_wt, len(lon_list[-1]))
                 # go again and fill with NaNs at end
-                for i in range(len(lon_list)): # get lons and lats
+                for i in range(len(lon_list)):  # get lons and lats
                     lon_list[i] = np.pad(lon_list[i], (0, max_pts_in_wt - len(lon_list[i])), mode='constant',
                                          constant_values=np.nan)
                     lat_list[i] = np.pad(lat_list[i], (0, max_pts_in_wt - len(lat_list[i])), mode='constant',
                                          constant_values=np.nan)
                     pres_list[i] = np.pad(pres_list[i], (0, max_pts_in_wt - len(pres_list[i])), mode='constant',
-                                        constant_values=np.nan)
+                                          constant_values=np.nan)
 
-                dataset_wt = dataset_wt.expand_dims(time=np.arange(0, max_pts_in_wt).astype(dtype=float)) # fake traj time
+                dataset_wt = dataset_wt.expand_dims(
+                    time=np.arange(0, max_pts_in_wt).astype(dtype=float))  # fake traj time
                 dataset_wt = dataset_wt.expand_dims(ensemble=[0])
                 dataset_wt = dataset_wt.expand_dims(trajectory=np.arange(1, len(ts.objects) + 1))
 
-                lons = xr.DataArray(np.zeros((1, len(ts.objects), max_pts_in_wt)), dims=("ensemble", "trajectory", "time"))
+                lons = xr.DataArray(np.zeros((1, len(ts.objects), max_pts_in_wt)),
+                                    dims=("ensemble", "trajectory", "time"))
                 lons.attrs['standard_name'] = "longitude"
                 lons.attrs['long_name'] = "longitude"
                 lons.attrs['units'] = "degrees_east"
@@ -341,13 +382,44 @@ class AEWIdentification(IdentificationTechnique):
                 dataset_wt['time'].attrs['long_name'] = "time"
                 dataset_wt['time'].attrs['units'] = "hours since " + ts.valid_time.replace('T', ' ')
                 dataset_wt['time'].attrs['trajectory_starttime'] = ts.valid_time.replace('T', ' ')
-                dataset_wt['time'].attrs['forecast_inittime'] = ts.valid_time.replace('T', ' ') # '2006-09-01 12:00:00' # TODO ts.valid_time.replace('T', ' ')
-                
-                out_path = self.config.out_traj_dir + ts.valid_time.replace(':','_') + '.nc'
+                dataset_wt['time'].attrs['forecast_inittime'] = ts.valid_time.replace('T',
+                                                                                      ' ')  # '2006-09-01 12:00:00' # TODO ts.valid_time.replace('T', ' ')
+
+                out_path = self.config.out_traj_dir + ts.valid_time.replace(':', '_') + '.nc'
                 dataset_wt.to_netcdf(out_path)
 
         return dataset, data_desc
 
+    # filter: keep current wavetrough if:
+    #   troughs are within a certain spatial window
+    #   length of the trough < threshold
+    def keep_wavetrough(self, properties):
+        """
+        Called for each wavetrough, check if kept based on filtering heuristics:
+        - WT requires any point of wavetrough in config.wave_filter range
+        - WT requires minimum length threshold
+
+        Parameters
+        ----------
+        properties
+
+        Returns
+        -------
+        True if kept
+        """
+
+        in_area = False
+        for line_pt in properties.line_pts:
+
+            # check if any point is outside filtering area
+            if (self.config.wave_filter_lon[0] < line_pt.lon < self.config.wave_filter_lon[1]
+                    and self.config.wave_filter_lat[0] < line_pt.lat < self.config.wave_filter_lat[1]):
+                in_area = True
 
+        if not in_area:  # no point of line segment in our area
+            return False
 
+        if properties.length_deg <= self.config.degree_len_thr:  # too small
+            return False
 
+        return True

enstools/feature/identification/african_easterly_waves/processing.py

@@ -106,14 +106,129 @@ def compute_cv(dataset, u_str, v_str, cv_str):
     return dataset
 
 
-def populate_object(obj_props, path):
+import math
+from enstools.feature.util.enstools_utils import get_u_var, get_v_var, get_vertical_dim, get_longitude_dim, \
+    get_latitude_dim
+
+
+
+# calculates the dx's and dy's for each grid cell
+# takes list of latitudes and longitudes as input and returns field with dimensions len(lats) x len(lons)
+def calculate_dx_dy(lats, lons):
+    lat_res = lats[1] - lats[0]
+    lon_res = lons[1] - lons[0]
+    nlat = len(lats)
+    nlon = len(lons)
+    lon_percent = (lon_res * nlon) / 360.0  # percentage of global longitudes. e.g. if 110W to 70E -> 0.5
+    lat_percent = (lat_res * nlat) / 180.0  # percentage of global latitudes
+    Re = 6378100
+
+    # field filled with the latitude value at each grid point
+    lat_f = np.tile(lats, (nlon, 1)).transpose()
+    # field filled with the longitude value at each grid point
+    lon_f = np.tile(lons, (nlat, 1))
+    # dx: circumference at current latitude  *  percentage of longitudes in dataset  /  amount of lon grid points
+    dx = np.cos(lat_f * math.pi / 180.0) * 2.0 * math.pi * Re * lon_percent / nlon
+    # dy: constant everywhere: half circumference  *  latitude percentage /  amount of lat points
+    dy = (lat_f * 0) + lat_percent * math.pi * Re / nlat
+
+    return dx, dy
+
+
+def compute_cv(dataset, u_str, v_str, cv_str):
+    xr.set_options(keep_attrs=True)  # assign coords keeps attributes
+
+    lon_str = get_longitude_dim(dataset)
+    lat_str = get_latitude_dim(dataset)
+
+    lats = dataset.coords[lat_str].data
+    lons = dataset.coords[lon_str].data
+
+    # be careful if your data is non-continuous by the chosen window (e.g. +120E..-120W)
+    # reorder axis for efficient numpy broadcasting
+    dataset = dataset.sortby(lat_str)
+    dataset = dataset.transpose(..., lat_str, lon_str)
+
+    dataset[cv_str] = xr.zeros_like(dataset[u_str], dtype=float)
+    dataset[cv_str].attrs = {'long_name': "Curvature Vorticity", 'units': "s**-1"}
+
+    # calculate dx and dy for each cell in grid in meters
+    dx, dy = calculate_dx_dy(lats, lons)
+
+    # Relative Vorticity = dv/dx - du/dy, use central differences
+    u_arr = dataset[u_str].data
+    v_arr = dataset[v_str].data
+    ndims = u_arr.ndim
+
+    # roll axes so we have a reference to each cell's neighbours
+    # ax[-1] = lon, ax[-2] = lat
+    # works as expected if longitude band is full 360 degrees. Rolling does exactly that.
+    v_xp = np.roll(v_arr, -1, axis=ndims - 1)  # roll -1 to get +1. v_xp = v(x+1)
+    v_xm = np.roll(v_arr, 1, axis=ndims - 1)
+    v_yp = np.roll(v_arr, -1, axis=ndims - 2)
+    v_ym = np.roll(v_arr, 1, axis=ndims - 2)
+
+    u_xp = np.roll(u_arr, -1, axis=ndims - 1)  # roll -1 to get +1
+    u_xm = np.roll(u_arr, 1, axis=ndims - 1)
+    u_yp = np.roll(u_arr, -1, axis=ndims - 2)
+    u_ym = np.roll(u_arr, 1, axis=ndims - 2)
+
+    # central differences:  (V(x+1)-V(x-1)) / 2x - (U(y+1)-U(y-1)) / 2y
+    RV = ((v_xp - v_xm) / (2 * dx)) - ((u_yp - u_ym) / (2 * dy))
+    # results verified to ERA5 vo-data
+
+    # split into shear + curvature: RV = -dV/dn + V/R
+    # shear vorticity:
+    # rate of change of wind speed in the direction of flow
+    # -dV/dn
+
+    # wind direction of this cell
+    ref_angle = np.arctan2(v_arr, u_arr)
+
+    ### METHOD 1 ### from analytical standpoint
+
+    sin_angle = np.sin(ref_angle)
+    cos_angle = np.cos(ref_angle)
+    # here is where the magic happens...
+    # to cartesian: dV/dn = - dV/dx * sin(phi) + dV/dy * cos(phi) ## n is normal to V, split it into x,y respecting the
+    # n-direction, where phi is the rotation angle of the natural coordinate system, this also is the angle of the
+    # wind vector
+    #
+    # we get the magnitude of shear as the projection of dV/dn onto the vector itself (direction e_t)
+    # with e_t = cos(phi)*e_x + sin(phi)*e_y , so:
+    # dV/dn * e_t = du/dx * (-sin*cos) + du/dy cos^2 + dv/dx *(-sin^2) + dv/dy (cos*sin)
+    SV = (u_xp - u_xm) / (2 * dx) * (-sin_angle * cos_angle) + \
+         (u_yp - u_ym) / (2 * dy) * (cos_angle ** 2) + \
+         (v_xp - v_xm) / (2 * dx) * (-(sin_angle ** 2)) + \
+         (v_yp - v_ym) / (2 * dy) * (sin_angle * cos_angle)
+    SV = -SV  # -dV/dn
+
+    # remainder is CV
+    CV = RV - SV
+    dataset[cv_str].values = CV
+
+    xr.set_options(keep_attrs='default')  # assign coords keeps attributes
+    return dataset
+
+
+def populate_object(obj_props, path, cfg):
     # obj_props.num_nodes = len(path)
     # fill the properties defined in the .proto file.
 
+    # first, remove vertices out of area
+    filtered_vertices = []
+    for v_idx, v in enumerate(path.vertices):
+        lon = v[0]
+        lat = v[1]
+
+        # dont use vertex not fulfilling area restriction
+        if cfg.is_point_in_area(lon, lat):
+            filtered_vertices.append(v)
+
     # vertices of path
     min_lat, max_lat, min_lon, max_lon = 90.0, -90.0, 180.0, -180.0
     dist_deg = 0.0
-    for v_idx, v in enumerate(path.vertices):
+    for v_idx, v in enumerate(filtered_vertices):
         line_pt = obj_props.line_pts.add()
         line_pt.lon = v[0]
         line_pt.lat = v[1]
@@ -129,7 +244,7 @@ def populate_object(obj_props, path):
 
         if v_idx > 0:
             dist_deg = dist_deg + (
-                        ((path.vertices[v_idx - 1][0] - v[0]) ** 2 + (path.vertices[v_idx - 1][1] - v[1]) ** 2) ** 0.5)
+                        ((filtered_vertices[v_idx - 1][0] - v[0]) ** 2 + (filtered_vertices[v_idx - 1][1] - v[1]) ** 2) ** 0.5)
 
     # bounding box
     obj_props.bb.min.lat = min_lat

enstools/feature/identification/identification.proto

@@ -43,11 +43,11 @@ message GraphNode {
 }
 
 /* Connections of the tracking graph. A connection is defined as one node (parent, the key object),
-   and all to its connected nodes. (list of childs, typically of the consecutive timestep)
+   and all to its connected nodes. (list of children, typically of the consecutive timestep)
  */
 message GraphConnection {
   required GraphNode parent = 1;
-  repeated GraphNode childs = 2;
+  repeated GraphNode children = 2;
 }
 
 /* The tracking graph is defined by a set of above connections.

enstools/feature/identification/overlap_example/example.py

@@ -11,7 +11,7 @@ from enstools.feature.util.graph import DataGraph
 # set the pb_reference to the compiled pb2.py file (see proto_gen directory)
 pipeline = FeaturePipeline(overlap_example_pb2, processing_mode='3d')
 
-# change this to an identification technique that actually does something: existing one or implement your own
+# change this to an identification strategy that actually does something: existing one or implement your own
 i_strat = OverlapIdentificationExample()  # set the Identification strategy
 t_strat = OverlapTracking(field_name='identified_areas')  # set the tracking strategy
 
@@ -37,7 +37,7 @@ for set in od.sets:  # for each set, e.g. each ensemble member, or each level
     for track in set.tracks: # for each identified track (run generate_tracks beforehand)
         single_track_graph = DataGraph(track, overlap_example_pb2)  # track as data graph (basically subgraph of above)
         # single_track_graph.get_parents(some node...)
-        # single_track_graph.get_childs(some node...)
+        # single_track_graph.get_children(some node...)
         
 
 pipeline.save_result(description_type='json', description_path=out_graph_path,

enstools/feature/identification/overlap_example/identification.py

-from ..identification import IdentificationTechnique
+from ..identification import IdentificationStrategy
 
 import xarray as xr
 from random import randrange
 import datetime
 
 
-class OverlapIdentificationExample(IdentificationTechnique):
+class OverlapIdentificationExample(IdentificationStrategy):
 
     def __init__(self, some_parameter='', **kwargs):
         # Constructor. Called from example_template.py, parameters can be passed and set here.

enstools/feature/identification/pv_streamer/identification.py

 from enstools.feature.identification.pv_streamer.processing import *
-from enstools.feature.identification import IdentificationTechnique
+from enstools.feature.identification import IdentificationStrategy
 import xarray as xr
 import numpy as np
 import copy
@@ -7,7 +7,7 @@ from enstools.feature.identification.pv_streamer.object_desc import get_object_d
 from enstools.feature.identification.pv_streamer.projection_cdo import project_latlon_to_stereo, project_stereo_to_latlon
 
 
-class PVIdentification(IdentificationTechnique):
+class PVIdentification(IdentificationStrategy):
 
     def __init__(self, unit='pv', mode_2d_layer=None, theta_range=None, extract_containing=None, centroid_lat_thr=None,
                  out_type='stereo', **kwargs):
@@ -104,8 +104,6 @@ class PVIdentification(IdentificationTechnique):
 
     # main identify, called on spatial 2d/3d subsets from enstools-feature.
     def identify(self, spatial_stereo_ds: xr.Dataset, **kwargs):
-        print("PV identify")
-
         levels_list = spatial_stereo_ds.coords[self.config.vertical_dim].values
 
         # preprocess binary fields: remove disturbances etc. in 2d, in 3d more flexible
@@ -164,6 +162,7 @@ class PVIdentification(IdentificationTechnique):
 
         spatial_stereo_ds['streamer'].values = labeled_areas
         # generate object descriptions
+
         object_list = get_object_data(self, spatial_stereo_ds, levels_list, dist_expand, self.area_map, self.config)
 
         # filter object descriptions