From 14b4ec685b1a942506679f7af870f9de05c7f439 Mon Sep 17 00:00:00 2001 From: gula Date: Sun, 18 Oct 2020 22:05:41 +0200 Subject: [PATCH 01/20] add modification to use Okubo-Weiss (and vorticity) instead of sla. --- src/py_eddy_tracker/dataset/grid.py | 60 ++++++++++++++++++++++++----- src/py_eddy_tracker/eddy_feature.py | 8 ++-- 2 files changed, 55 insertions(+), 13 deletions(-) diff --git a/src/py_eddy_tracker/dataset/grid.py b/src/py_eddy_tracker/dataset/grid.py index e87ba747..674a4ee6 100644 --- a/src/py_eddy_tracker/dataset/grid.py +++ b/src/py_eddy_tracker/dataset/grid.py @@ -578,6 +578,7 @@ def eddy_identification( precision=None, force_height_unit=None, force_speed_unit=None, + vorticity_name='vrt', mle=1, ): """ @@ -616,7 +617,11 @@ def eddy_identification( self.units(grid_height) if force_height_unit is None else force_height_unit ) units = UnitRegistry() - in_h_unit = units.parse_expression(h_units) + if grid_height in ['ow']: + in_h_unit = None + else: + in_h_unit = units.parse_expression(h_units) + if in_h_unit is not None: factor, _ = in_h_unit.to("m").to_tuple() logger.info( @@ -629,14 +634,31 @@ def eddy_identification( # Get h grid data = self.grid(grid_height).astype("f8") + + if grid_height in ['ow']: + # Get vorticity as an aditional field (to identify cyc/acyc) + vrt = self.grid(vorticity_name, indexs=dict(xi_u=slice(None),eta_v=slice(None))).astype('f8') + + if vorticity_name=='vrt': vrt = self.psi2rho(vrt) + + print(vrt.shape,data.shape) + # In case of a reduce mask if len(data.mask.shape) == 0 and not data.mask: data.mask = zeros(data.shape, dtype="bool") + data.mask[isnan(data)] = 1 # we remove noisy information if precision is not None: data = (data / precision).round() * precision # Compute levels for ssh - z_min, z_max = data.min(), data.max() + + if grid_height in ['ow']: + #z_min, z_max = -2e-10, -0.5e-10#-0.2 * data.std() + z_min, z_max = -0.3 * data.std(),-0.2 * data.std() + else: + z_min, z_max = data.min(), data.max() + print('init z_min, z_max, step',z_min, z_max, step) + print(data.mask.sum()) d_z = z_max - z_min data_tmp = data[~data.mask] epsilon = 0.001 # in % @@ -654,8 +676,9 @@ def eddy_identification( z_max_p, ) z_min, z_max = z_min_p, z_max_p - + print('step1 z_min, z_max, step',z_min, z_max, step) #debug levels = arange(z_min - z_min % step, z_max - z_max % step + 2 * step, step) + print('levels',levels) #debug # Get x and y values x, y = self.x_c, self.y_c @@ -682,7 +705,10 @@ def eddy_identification( a_and_c = list() for anticyclonic_search in [True, False]: eddies = list() - iterator = 1 if anticyclonic_search else -1 + if grid_height in ['ow']: + iterator = -1 + else: + iterator = 1 if anticyclonic_search else -1 # Loop over each collection for coll_ind, coll in enumerate(self.contours.iter(step=iterator)): @@ -724,10 +750,17 @@ def eddy_identification( continue # Test to know cyclone or anticyclone - if has_value( + if grid_height in ['ow']: + #acyc_not_cyc = vrt[i_x, i_y] <= 0 + if has_value( + vrt, i_x_in, i_y_in, 0, below=not anticyclonic_search + ): + continue + else: + if has_value( data, i_x_in, i_y_in, cvalues, below=anticyclonic_search ): - continue + continue # FIXME : Maybe limit max must be replace with a maximum of surface if ( @@ -744,7 +777,7 @@ def eddy_identification( data, anticyclonic_search=anticyclonic_search, level=self.contours.levels[corrected_coll_index], - step=step, + step=step, grid_height=grid_height, mle=mle, ) # If we have a valid amplitude @@ -980,6 +1013,7 @@ def get_amplitude( anticyclonic_search=True, level=None, step=None, + grid_height='sla', mle=1, ): # Instantiate Amplitude object @@ -996,10 +1030,16 @@ def get_amplitude( mle=mle, ) - if anticyclonic_search: - reset_centroid = amp.all_pixels_above_h0(level) + if grid_height in ['ow']: + if anticyclonic_search: + reset_centroid = amp.all_pixels_below_h0(level,grid_height='ow') + else: + reset_centroid = amp.all_pixels_below_h0(level,grid_height='ow') else: - reset_centroid = amp.all_pixels_below_h0(level) + if anticyclonic_search: + reset_centroid = amp.all_pixels_above_h0(level) + else: + reset_centroid = amp.all_pixels_below_h0(level) return reset_centroid, amp diff --git a/src/py_eddy_tracker/eddy_feature.py b/src/py_eddy_tracker/eddy_feature.py index f039b097..66176dc9 100644 --- a/src/py_eddy_tracker/eddy_feature.py +++ b/src/py_eddy_tracker/eddy_feature.py @@ -87,7 +87,7 @@ def within_amplitude_limits(self): """ return self.interval_min <= self.amplitude - def all_pixels_below_h0(self, level): + def all_pixels_below_h0(self, level, grid_height='sla'): """ Check CSS11 criterion 1: The SSH values of all of the pixels are below a given SSH threshold for cyclonic eddies. @@ -120,8 +120,10 @@ def all_pixels_below_h0(self, level): nb_real_extrema = ( (level - self.grid_extract.data[lmi_i, lmi_j]) >= self.interval_min ).sum() - if nb_real_extrema > self.mle: - return False + + if grid_height!='ow': + if nb_real_extrema > self.mle: + return False index = self.grid_extract.data[lmi_i, lmi_j].argmin() i, j = lmi_i[index], lmi_j[index] self.amplitude = abs(self.grid_extract.data[i, j] - self.h_0) From c5449bee483dd6963756cfe5c57a7f30cc2c230c Mon Sep 17 00:00:00 2001 From: gula Date: Mon, 18 Jan 2021 22:36:37 +0100 Subject: [PATCH 02/20] use nanstd to compute OW levels --- src/py_eddy_tracker/dataset/grid.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/py_eddy_tracker/dataset/grid.py b/src/py_eddy_tracker/dataset/grid.py index 674a4ee6..08482820 100644 --- a/src/py_eddy_tracker/dataset/grid.py +++ b/src/py_eddy_tracker/dataset/grid.py @@ -32,6 +32,7 @@ round_, nanmean, exp, + nanstd, mean as np_mean, ) from datetime import datetime @@ -654,7 +655,7 @@ def eddy_identification( if grid_height in ['ow']: #z_min, z_max = -2e-10, -0.5e-10#-0.2 * data.std() - z_min, z_max = -0.3 * data.std(),-0.2 * data.std() + z_min, z_max = -0.3 * nanstd(data),-0.2 * nanstd(data) else: z_min, z_max = data.min(), data.max() print('init z_min, z_max, step',z_min, z_max, step) From 26f2eba166b46e93a4d43c3875e8987145ad9723 Mon Sep 17 00:00:00 2001 From: jeremy collin Date: Mon, 8 Feb 2021 15:58:14 +0100 Subject: [PATCH 03/20] added z_min, z_max parameters for eddy_detection --- src/py_eddy_tracker/dataset/grid.py | 57 +++++++++++++++-------------- 1 file changed, 29 insertions(+), 28 deletions(-) diff --git a/src/py_eddy_tracker/dataset/grid.py b/src/py_eddy_tracker/dataset/grid.py index 674a4ee6..a37b9cd1 100644 --- a/src/py_eddy_tracker/dataset/grid.py +++ b/src/py_eddy_tracker/dataset/grid.py @@ -32,6 +32,7 @@ round_, nanmean, exp, + nanstd, mean as np_mean, ) from datetime import datetime @@ -571,6 +572,8 @@ def eddy_identification( uname, vname, date, + z_min=None, + z_max=None, step=0.005, shape_error=55, sampling=50, @@ -589,6 +592,8 @@ def eddy_identification( :param str vname: Grid name of v speed component :param datetime.datetime date: Date which will be store in object to date data :param float,int step: Height between two layers in m + :param float z_min : minimum value foe eddy detection + :param float z_max : maximum value for eddy detection :param float,int shape_error: Maximal error allow for outter contour in % :param int sampling: Sampling of contour and speed profile :param (int,int),None pixel_limit: @@ -640,8 +645,6 @@ def eddy_identification( vrt = self.grid(vorticity_name, indexs=dict(xi_u=slice(None),eta_v=slice(None))).astype('f8') if vorticity_name=='vrt': vrt = self.psi2rho(vrt) - - print(vrt.shape,data.shape) # In case of a reduce mask if len(data.mask.shape) == 0 and not data.mask: @@ -650,35 +653,33 @@ def eddy_identification( # we remove noisy information if precision is not None: data = (data / precision).round() * precision - # Compute levels for ssh - - if grid_height in ['ow']: - #z_min, z_max = -2e-10, -0.5e-10#-0.2 * data.std() - z_min, z_max = -0.3 * data.std(),-0.2 * data.std() - else: - z_min, z_max = data.min(), data.max() - print('init z_min, z_max, step',z_min, z_max, step) - print(data.mask.sum()) - d_z = z_max - z_min - data_tmp = data[~data.mask] - epsilon = 0.001 # in % - z_min_p, z_max_p = ( - percentile(data_tmp, epsilon), - percentile(data_tmp, 100 - epsilon), - ) - d_zp = z_max_p - z_min_p - if d_z / d_zp > 2: - logger.warning( - "Maybe some extrema are present zmin %f (m) and zmax %f (m) will be replace by %f and %f", - z_min, - z_max, - z_min_p, - z_max_p, + # Compute levels for ssh/okubo + if z_min is None or z_max is None: + if grid_height in ['ow']: + z_min, z_max = -0.3 * 1e-9, -0.2 * 1e-9 + #z_min, z_max = -0.3 * nanstd(data),-0.2 * nanstd(data) + else: + z_min, z_max = data.min(), data.max() + d_z = z_max - z_min + data_tmp = data[~data.mask] + epsilon = 0.001 # in % + z_min_p, z_max_p = ( + percentile(data_tmp, epsilon), + percentile(data_tmp, 100 - epsilon), ) - z_min, z_max = z_min_p, z_max_p + d_zp = z_max_p - z_min_p + if d_z / d_zp > 2: + logger.warning( + "Maybe some extrema are present zmin %f (m) and zmax %f (m) will be replace by %f and %f", + z_min, + z_max, + z_min_p, + z_max_p, + ) + z_min, z_max = z_min_p, z_max_p + print('step1 z_min, z_max, step',z_min, z_max, step) #debug levels = arange(z_min - z_min % step, z_max - z_max % step + 2 * step, step) - print('levels',levels) #debug # Get x and y values x, y = self.x_c, self.y_c From e65566b3acccebe1f2135be31034dbaa028d5ceb Mon Sep 17 00:00:00 2001 From: jeremy collin Date: Tue, 9 Feb 2021 17:55:05 +0100 Subject: [PATCH 04/20] detection: changed time dtype in netcdf to float32 --- src/py_eddy_tracker/__init__.py | 2 +- src/py_eddy_tracker/dataset/grid.py | 1 - 2 files changed, 1 insertion(+), 2 deletions(-) diff --git a/src/py_eddy_tracker/__init__.py b/src/py_eddy_tracker/__init__.py index 04324267..9621c7d8 100644 --- a/src/py_eddy_tracker/__init__.py +++ b/src/py_eddy_tracker/__init__.py @@ -100,7 +100,7 @@ def parse_args(self, *args, **kwargs): attr_name="time", nc_name="time", old_nc_name=["j1"], - nc_type="int32", + nc_type="float32", nc_dims=("obs",), nc_attr=dict( standard_name="time", diff --git a/src/py_eddy_tracker/dataset/grid.py b/src/py_eddy_tracker/dataset/grid.py index 09d56683..2c237625 100644 --- a/src/py_eddy_tracker/dataset/grid.py +++ b/src/py_eddy_tracker/dataset/grid.py @@ -617,7 +617,6 @@ def eddy_identification( self.init_speed_coef(uname, vname) # Get unit of h grid - h_units = ( self.units(grid_height) if force_height_unit is None else force_height_unit ) From dffc8500205a27a7542fab68a8c806a767073609 Mon Sep 17 00:00:00 2001 From: jeremy collin Date: Mon, 1 Mar 2021 11:40:46 +0100 Subject: [PATCH 05/20] tracking: date to datetime [hourly] --- src/scripts/EddyTracking | 25 ++++++++++++++----------- 1 file changed, 14 insertions(+), 11 deletions(-) diff --git a/src/scripts/EddyTracking b/src/scripts/EddyTracking index a1c268bc..806a327e 100644 --- a/src/scripts/EddyTracking +++ b/src/scripts/EddyTracking @@ -11,6 +11,7 @@ from os import mkdir from re import compile as re_compile from os.path import join as join_path from numpy import bytes_, empty, unique +from numpy import dtype as npdtype from netCDF4 import Dataset from datetime import datetime from glob import glob @@ -39,17 +40,17 @@ def browse_dataset_in( filenames = bytes_(glob(full_path)) dataset_list = empty( - len(filenames), dtype=[("filename", "S500"), ("date", "datetime64[D]"),] + len(filenames), dtype=[("filename", "S500"), ("datetime", npdtype('M8[us]')),] ) dataset_list["filename"] = filenames logger.info("%s grids available", dataset_list.shape[0]) mode_attrs = False if "(" not in date_regexp: - logger.debug("Attrs date : %s", date_regexp) + logger.debug("Attrs datetime : %s", date_regexp) mode_attrs = date_regexp.strip().split(":") else: - logger.debug("Pattern date : %s", date_regexp) + logger.debug("Pattern datetime : %s", date_regexp) for item in dataset_list: str_date = None @@ -65,11 +66,13 @@ def browse_dataset_in( str_date = result.groups()[0] if str_date is not None: - item["date"] = datetime.strptime(str_date, date_model).date() + #try: + # item["date"] = datetime.strptime(str_date, date_model).date() + item["datetime"] = datetime.strptime(str_date, date_model) - dataset_list.sort(order=["date", "filename"]) + dataset_list.sort(order=["datetime", "filename"]) - steps = unique(dataset_list["date"][1:] - dataset_list["date"][:-1]) + steps = unique(dataset_list["datetime"][1:] - dataset_list["datetime"][:-1]) if len(steps) > 1: raise Exception("Several days steps in grid dataset %s" % steps) @@ -80,11 +83,11 @@ def browse_dataset_in( if start_date is not None or end_date is not None: logger.info( "Available grid from %s to %s", - dataset_list[0]["date"], - dataset_list[-1]["date"], + dataset_list[0]["datetime"], + dataset_list[-1]["datetime"], ) logger.info("Filtering grid by time %s, %s", start_date, end_date) - mask = (dataset_list["date"] >= start_date) * (dataset_list["date"] <= end_date) + mask = (dataset_list["datetime"] >= start_date) * (dataset_list["datetime"] <= end_date) dataset_list = dataset_list[mask] return dataset_list @@ -171,14 +174,14 @@ if __name__ == "__main__": files_model=None, files=CONFIG["PATHS"]["FILES_PATTERN"], date_regexp=".*_([0-9]*?).[nz].*", - date_model="%Y%m%d", + date_model="%Y%m%d%H", ) else: DATASET_LIST = browse_dataset_in( data_dir=dirname(CONFIG["PATHS"]["FILES_PATTERN"]), files_model=basename(CONFIG["PATHS"]["FILES_PATTERN"]), date_regexp=".*_([0-9]*?).[nz].*", - date_model="%Y%m%d", + date_model="%Y%m%d%H", ) if BLANK_PERIOD > 0: From b69799c8e09d875a4babdfd2999779464df14ee7 Mon Sep 17 00:00:00 2001 From: jeremy collin Date: Mon, 1 Mar 2021 15:38:58 +0100 Subject: [PATCH 06/20] tracking: skiped len(steps)>1 exception --- src/scripts/EddyTracking | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/scripts/EddyTracking b/src/scripts/EddyTracking index 806a327e..ddaeae71 100644 --- a/src/scripts/EddyTracking +++ b/src/scripts/EddyTracking @@ -74,7 +74,8 @@ def browse_dataset_in( steps = unique(dataset_list["datetime"][1:] - dataset_list["datetime"][:-1]) if len(steps) > 1: - raise Exception("Several days steps in grid dataset %s" % steps) + print("Several days steps in grid dataset %s" % steps) + #raise Exception("Several days steps in grid dataset %s" % steps) if sub_sampling_step != 1: logger.info("Grid subsampling %d", sub_sampling_step) From d6b7a86641b2d8ddbc73a808a33e7f20bd86ad2d Mon Sep 17 00:00:00 2001 From: Louise BELEY Date: Tue, 2 Mar 2021 11:00:35 +0000 Subject: [PATCH 07/20] indentation vrt fix --- src/py_eddy_tracker/dataset/grid.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/py_eddy_tracker/dataset/grid.py b/src/py_eddy_tracker/dataset/grid.py index 2c237625..bb6fabcd 100644 --- a/src/py_eddy_tracker/dataset/grid.py +++ b/src/py_eddy_tracker/dataset/grid.py @@ -642,8 +642,7 @@ def eddy_identification( if grid_height in ['ow']: # Get vorticity as an aditional field (to identify cyc/acyc) vrt = self.grid(vorticity_name, indexs=dict(xi_u=slice(None),eta_v=slice(None))).astype('f8') - - if vorticity_name=='vrt': vrt = self.psi2rho(vrt) + if vorticity_name=='vrt': vrt = self.psi2rho(vrt) # In case of a reduce mask if len(data.mask.shape) == 0 and not data.mask: From 3147088aa239abeb041849218a64da085374d25e Mon Sep 17 00:00:00 2001 From: gula Date: Wed, 14 Apr 2021 11:48:24 +0200 Subject: [PATCH 08/20] solve merge issue (interval instead of step) --- src/py_eddy_tracker/dataset/grid.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/py_eddy_tracker/dataset/grid.py b/src/py_eddy_tracker/dataset/grid.py index b44ff38b..7ca743a0 100644 --- a/src/py_eddy_tracker/dataset/grid.py +++ b/src/py_eddy_tracker/dataset/grid.py @@ -811,7 +811,7 @@ def eddy_identification( data, anticyclonic_search=anticyclonic_search, level=self.contours.levels[corrected_coll_index], - step=step, grid_height=grid_height, + interval=step, grid_height=grid_height, mle=mle, **kwargs, ) @@ -1039,7 +1039,7 @@ def get_amplitude( data, anticyclonic_search=True, level=None, - step=None, + interval=None, grid_height='sla', mle=1, **kwargs From 38f7bf287bca370bf24f8c66e0c6a42196e00764 Mon Sep 17 00:00:00 2001 From: gula Date: Wed, 14 Apr 2021 14:30:26 +0200 Subject: [PATCH 09/20] correct issues with datetime after merge --- src/py_eddy_tracker/appli/eddies.py | 25 +-- src/py_eddy_tracker/dataset/grid.py | 1 - src/scripts/EddyTracking | 274 ++++++++++++++++++++++++++++ 3 files changed, 288 insertions(+), 12 deletions(-) create mode 100644 src/scripts/EddyTracking diff --git a/src/py_eddy_tracker/appli/eddies.py b/src/py_eddy_tracker/appli/eddies.py index d5a727f9..05c62b5c 100644 --- a/src/py_eddy_tracker/appli/eddies.py +++ b/src/py_eddy_tracker/appli/eddies.py @@ -12,7 +12,9 @@ from re import compile as re_compile from netCDF4 import Dataset +from datetime import datetime from numpy import bincount, bytes_, empty, in1d, unique +from numpy import dtype as npdtype from yaml import safe_load from .. import EddyParser @@ -241,7 +243,7 @@ def browse_dataset_in( len(filenames), dtype=[ ("filename", "S500"), - ("date", "datetime64[D]"), + ("datetime", npdtype('M8[us]')), ], ) dataset_list["filename"] = filenames @@ -249,10 +251,10 @@ def browse_dataset_in( logger.info("%s grids available", dataset_list.shape[0]) mode_attrs = False if "(" not in date_regexp: - logger.debug("Attrs date : %s", date_regexp) + logger.debug("Attrs datetime : %s", date_regexp) mode_attrs = date_regexp.strip().split(":") else: - logger.debug("Pattern date : %s", date_regexp) + logger.debug("Pattern datetime : %s", date_regexp) for item in dataset_list: str_date = None @@ -268,13 +270,14 @@ def browse_dataset_in( str_date = result.groups()[0] if str_date is not None: - item["date"] = datetime.strptime(str_date, date_model).date() + item["datetime"] = datetime.strptime(str_date, date_model) - dataset_list.sort(order=["date", "filename"]) + dataset_list.sort(order=["datetime", "filename"]) - steps = unique(dataset_list["date"][1:] - dataset_list["date"][:-1]) + steps = unique(dataset_list["datetime"][1:] - dataset_list["datetime"][:-1]) if len(steps) > 1: - raise Exception("Several days steps in grid dataset %s" % steps) + print("Several days steps in grid dataset %s" % steps) + #raise Exception("Several days steps in grid dataset %s" % steps) if sub_sampling_step != 1: logger.info("Grid subsampling %d", sub_sampling_step) @@ -283,11 +286,11 @@ def browse_dataset_in( if start_date is not None or end_date is not None: logger.info( "Available grid from %s to %s", - dataset_list[0]["date"], - dataset_list[-1]["date"], + dataset_list[0]["datetime"], + dataset_list[-1]["datetime"], ) logger.info("Filtering grid by time %s, %s", start_date, end_date) - mask = (dataset_list["date"] >= start_date) * (dataset_list["date"] <= end_date) + mask = (dataset_list["datetime"] >= start_date) * (dataset_list["datetime"] <= end_date) dataset_list = dataset_list[mask] return dataset_list @@ -304,7 +307,7 @@ def track( correspondances_only=False, **kw_c, ): - kw = dict(date_regexp=".*_([0-9]*?).[nz].*", date_model="%Y%m%d") + kw = dict(date_regexp=".*_([0-9]*?).[nz].*", date_model="%Y%m%d%H") if isinstance(pattern, list): kw.update(dict(data_dir=None, files_model=None, files=pattern)) else: diff --git a/src/py_eddy_tracker/dataset/grid.py b/src/py_eddy_tracker/dataset/grid.py index 7ca743a0..b5a513fd 100644 --- a/src/py_eddy_tracker/dataset/grid.py +++ b/src/py_eddy_tracker/dataset/grid.py @@ -1039,7 +1039,6 @@ def get_amplitude( data, anticyclonic_search=True, level=None, - interval=None, grid_height='sla', mle=1, **kwargs diff --git a/src/scripts/EddyTracking b/src/scripts/EddyTracking new file mode 100644 index 00000000..ddaeae71 --- /dev/null +++ b/src/scripts/EddyTracking @@ -0,0 +1,274 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- +""" +Track eddy with Identification file produce with EddyIdentification +""" +from py_eddy_tracker import EddyParser +from yaml import load as yaml_load +from py_eddy_tracker.tracking import Correspondances +from os.path import exists, dirname, basename +from os import mkdir +from re import compile as re_compile +from os.path import join as join_path +from numpy import bytes_, empty, unique +from numpy import dtype as npdtype +from netCDF4 import Dataset +from datetime import datetime +from glob import glob +import logging + + +logger = logging.getLogger("pet") + + +def browse_dataset_in( + data_dir, + files_model, + date_regexp, + date_model, + start_date=None, + end_date=None, + sub_sampling_step=1, + files=None, +): + pattern_regexp = re_compile(".*/" + date_regexp) + if files is not None: + filenames = bytes_(files) + else: + full_path = join_path(data_dir, files_model) + logger.info("Search files : %s", full_path) + filenames = bytes_(glob(full_path)) + + dataset_list = empty( + len(filenames), dtype=[("filename", "S500"), ("datetime", npdtype('M8[us]')),] + ) + dataset_list["filename"] = filenames + + logger.info("%s grids available", dataset_list.shape[0]) + mode_attrs = False + if "(" not in date_regexp: + logger.debug("Attrs datetime : %s", date_regexp) + mode_attrs = date_regexp.strip().split(":") + else: + logger.debug("Pattern datetime : %s", date_regexp) + + for item in dataset_list: + str_date = None + if mode_attrs: + with Dataset(item["filename"].decode("utf-8")) as h: + if len(mode_attrs) == 1: + str_date = getattr(h, mode_attrs[0]) + else: + str_date = getattr(h.variables[mode_attrs[0]], mode_attrs[1]) + else: + result = pattern_regexp.match(str(item["filename"])) + if result: + str_date = result.groups()[0] + + if str_date is not None: + #try: + # item["date"] = datetime.strptime(str_date, date_model).date() + item["datetime"] = datetime.strptime(str_date, date_model) + + dataset_list.sort(order=["datetime", "filename"]) + + steps = unique(dataset_list["datetime"][1:] - dataset_list["datetime"][:-1]) + if len(steps) > 1: + print("Several days steps in grid dataset %s" % steps) + #raise Exception("Several days steps in grid dataset %s" % steps) + + if sub_sampling_step != 1: + logger.info("Grid subsampling %d", sub_sampling_step) + dataset_list = dataset_list[::sub_sampling_step] + + if start_date is not None or end_date is not None: + logger.info( + "Available grid from %s to %s", + dataset_list[0]["datetime"], + dataset_list[-1]["datetime"], + ) + logger.info("Filtering grid by time %s, %s", start_date, end_date) + mask = (dataset_list["datetime"] >= start_date) * (dataset_list["datetime"] <= end_date) + + dataset_list = dataset_list[mask] + return dataset_list + + +def usage(): + """Usage + """ + # Run using: + parser = EddyParser("Tool to use identification step to compute tracking") + parser.add_argument("yaml_file", help="Yaml file to configure py-eddy-tracker") + parser.add_argument("--correspondance_in", help="Filename of saved correspondance") + parser.add_argument("--correspondance_out", help="Filename to save correspondance") + parser.add_argument( + "--save_correspondance_and_stop", + action="store_true", + help="Stop tracking after correspondance computation," + " merging can be done with EddyFinalTracking", + ) + parser.add_argument( + "--zarr", action="store_true", help="Output will be wrote in zarr" + ) + parser.add_argument("--unraw", action="store_true", help="Load unraw data") + parser.add_argument( + "--blank_period", + type=int, + default=0, + help="Nb of detection which will not use at the end of the period", + ) + args = parser.parse_args() + + # Read yaml configuration file + with open(args.yaml_file, "r") as stream: + config = yaml_load(stream) + if args.correspondance_in is not None and not exists(args.correspondance_in): + args.correspondance_in = None + return ( + config, + args.save_correspondance_and_stop, + args.correspondance_in, + args.correspondance_out, + args.blank_period, + args.zarr, + not args.unraw, + ) + + +if __name__ == "__main__": + ( + CONFIG, + SAVE_STOP, + CORRESPONDANCES_IN, + CORRESPONDANCES_OUT, + BLANK_PERIOD, + ZARR, + RAW, + ) = usage() + # Create output directory + SAVE_DIR = CONFIG["PATHS"].get("SAVE_DIR", None) + if SAVE_DIR is not None and not exists(SAVE_DIR): + mkdir(SAVE_DIR) + + YAML_CORRESPONDANCES_OUT = CONFIG["PATHS"].get("CORRESPONDANCES_OUT", None) + if CORRESPONDANCES_IN is None: + CORRESPONDANCES_IN = CONFIG["PATHS"].get("CORRESPONDANCES_IN", None) + if CORRESPONDANCES_OUT is None: + CORRESPONDANCES_OUT = YAML_CORRESPONDANCES_OUT + if YAML_CORRESPONDANCES_OUT is None and CORRESPONDANCES_OUT is None: + CORRESPONDANCES_OUT = "{path}/{sign_type}_correspondances.nc" + + CLASS = None + if "CLASS" in CONFIG: + CLASSNAME = CONFIG["CLASS"]["CLASS"] + CLASS = getattr( + __import__(CONFIG["CLASS"]["MODULE"], globals(), locals(), CLASSNAME), + CLASSNAME, + ) + + NB_VIRTUAL_OBS_MAX_BY_SEGMENT = int(CONFIG.get("VIRTUAL_LENGTH_MAX", 0)) + + if isinstance(CONFIG["PATHS"]["FILES_PATTERN"], list): + DATASET_LIST = browse_dataset_in( + data_dir=None, + files_model=None, + files=CONFIG["PATHS"]["FILES_PATTERN"], + date_regexp=".*_([0-9]*?).[nz].*", + date_model="%Y%m%d%H", + ) + else: + DATASET_LIST = browse_dataset_in( + data_dir=dirname(CONFIG["PATHS"]["FILES_PATTERN"]), + files_model=basename(CONFIG["PATHS"]["FILES_PATTERN"]), + date_regexp=".*_([0-9]*?).[nz].*", + date_model="%Y%m%d%H", + ) + + if BLANK_PERIOD > 0: + DATASET_LIST = DATASET_LIST[:-BLANK_PERIOD] + logger.info("Last %d files will be pop", BLANK_PERIOD) + + START_TIME = datetime.now() + logger.info("Start tracking on %d files", len(DATASET_LIST)) + + NB_OBS_MIN = int(CONFIG.get("TRACK_DURATION_MIN", 14)) + if NB_OBS_MIN > len(DATASET_LIST): + raise Exception( + "Input file number (%s) is shorter than TRACK_DURATION_MIN (%s)." + % (len(DATASET_LIST), NB_OBS_MIN) + ) + + CORRESPONDANCES = Correspondances( + datasets=DATASET_LIST["filename"], + virtual=NB_VIRTUAL_OBS_MAX_BY_SEGMENT, + class_method=CLASS, + previous_correspondance=CORRESPONDANCES_IN, + ) + CORRESPONDANCES.track() + logger.info("Track finish") + + logger.info("Start merging") + DATE_START, DATE_STOP = CORRESPONDANCES.period + DICT_COMPLETION = dict( + date_start=DATE_START, + date_stop=DATE_STOP, + date_prod=START_TIME, + path=SAVE_DIR, + sign_type=CORRESPONDANCES.current_obs.sign_legend, + ) + + CORRESPONDANCES.save(CORRESPONDANCES_OUT, DICT_COMPLETION) + if SAVE_STOP: + exit() + + # Merge correspondance, only do if we stop and store just after compute of correspondance + CORRESPONDANCES.prepare_merging() + + logger.info( + "Longer track saved have %d obs", CORRESPONDANCES.nb_obs_by_tracks.max() + ) + logger.info( + "The mean length is %d observations before filtering", + CORRESPONDANCES.nb_obs_by_tracks.mean(), + ) + + CORRESPONDANCES.get_unused_data(raw_data=RAW).write_file( + path=SAVE_DIR, filename="%(path)s/%(sign_type)s_untracked.nc", zarr_flag=ZARR + ) + + SHORT_CORRESPONDANCES = CORRESPONDANCES._copy() + SHORT_CORRESPONDANCES.shorter_than(size_max=NB_OBS_MIN) + + CORRESPONDANCES.longer_than(size_min=NB_OBS_MIN) + + FINAL_EDDIES = CORRESPONDANCES.merge(raw_data=RAW) + SHORT_TRACK = SHORT_CORRESPONDANCES.merge(raw_data=RAW) + + # We flag obs + if CORRESPONDANCES.virtual: + FINAL_EDDIES["virtual"][:] = FINAL_EDDIES["time"] == 0 + FINAL_EDDIES.filled_by_interpolation(FINAL_EDDIES["virtual"] == 1) + SHORT_TRACK["virtual"][:] = SHORT_TRACK["time"] == 0 + SHORT_TRACK.filled_by_interpolation(SHORT_TRACK["virtual"] == 1) + + # Total running time + FULL_TIME = datetime.now() - START_TIME + logger.info("Mean duration by loop : %s", FULL_TIME / (len(DATASET_LIST) - 1)) + logger.info("Duration : %s", FULL_TIME) + + logger.info( + "Longer track saved have %d obs", CORRESPONDANCES.nb_obs_by_tracks.max() + ) + logger.info( + "The mean length is %d observations after filtering", + CORRESPONDANCES.nb_obs_by_tracks.mean(), + ) + + FINAL_EDDIES.write_file(path=SAVE_DIR, zarr_flag=ZARR) + SHORT_TRACK.write_file( + filename="%(path)s/%(sign_type)s_track_too_short.nc", + path=SAVE_DIR, + zarr_flag=ZARR, + ) + From e479827a9950ba18ac4f89faf94c3c2d732af10f Mon Sep 17 00:00:00 2001 From: gula Date: Thu, 6 May 2021 12:08:04 +0200 Subject: [PATCH 10/20] update grid for zooming with roms outputs and ow --- src/py_eddy_tracker/dataset/grid.py | 14 +++++++++----- 1 file changed, 9 insertions(+), 5 deletions(-) diff --git a/src/py_eddy_tracker/dataset/grid.py b/src/py_eddy_tracker/dataset/grid.py index b5a513fd..1d23fac0 100644 --- a/src/py_eddy_tracker/dataset/grid.py +++ b/src/py_eddy_tracker/dataset/grid.py @@ -385,6 +385,7 @@ def load(self): """ x_name, y_name = self.coordinates with Dataset(self.filename) as h: + self.x_dim = h.variables[x_name].dimensions self.y_dim = h.variables[y_name].dimensions @@ -665,11 +666,11 @@ def eddy_identification( precision /= factor # Get ssh grid - data = self.grid(grid_height).astype("f8") + data = 1e10*self.grid(grid_height).astype("f8") if grid_height in ['ow']: # Get vorticity as an aditional field (to identify cyc/acyc) - vrt = self.grid(vorticity_name, indexs=dict(xi_u=slice(None),eta_v=slice(None))).astype('f8') + vrt = self.grid(vorticity_name, indexs=self.indexs).astype('f8') if vorticity_name=='vrt': vrt = self.psi2rho(vrt) # In case of a reduced mask @@ -706,7 +707,7 @@ def eddy_identification( z_min, z_max = z_min_p, z_max_p print('step1 z_min, z_max, step',z_min, z_max, step) #debug - levels = arange(z_min - z_min % step, z_max - z_max % step + 2 * step, step) + levels = arange(z_min - z_min % step, z_max - z_max % step + step, step) # Get x and y values x, y = self.x_c, self.y_c @@ -803,7 +804,6 @@ def eddy_identification( # Here the considered contour passed shape_error test, masked_pixels test, # values strictly above (AEs) or below (CEs) the contour, number_pixels test) - # Compute amplitude reset_centroid, amp = self.get_amplitude( contour, @@ -815,6 +815,7 @@ def eddy_identification( mle=mle, **kwargs, ) + # If we have a valid amplitude if (not amp.within_amplitude_limits()) or (amp.amplitude == 0): contour.reject = 4 @@ -985,6 +986,7 @@ def get_uavg( all_contours.iter(start=level_start + step, step=step) ): level_contour = coll.get_nearest_path_bbox_contain_pt(centlon_e, centlat_e) + # Leave loop if no contours at level if level_contour is None: break @@ -1078,6 +1080,7 @@ class UnRegularGridDataset(GridDataset): def load(self): """Load variable (data)""" x_name, y_name = self.coordinates + with Dataset(self.filename) as h: self.x_dim = h.variables[x_name].dimensions self.y_dim = h.variables[y_name].dimensions @@ -1091,7 +1094,8 @@ def load(self): self.y_c = self.vars[y_name] self.init_pos_interpolator() - + + @property def bounds(self): """Give bounds""" From f0833e0e39c0cd2fc803ada17623812b4c6c8d4d Mon Sep 17 00:00:00 2001 From: gula Date: Thu, 6 May 2021 12:09:46 +0200 Subject: [PATCH 11/20] Create detection_gigatl6.py --- examples/detection_gigatl6.py | 124 ++++++++++++++++++++++++++++++++++ 1 file changed, 124 insertions(+) create mode 100644 examples/detection_gigatl6.py diff --git a/examples/detection_gigatl6.py b/examples/detection_gigatl6.py new file mode 100644 index 00000000..7d151083 --- /dev/null +++ b/examples/detection_gigatl6.py @@ -0,0 +1,124 @@ +from datetime import datetime, timedelta +from py_eddy_tracker import start_logger +from py_eddy_tracker.dataset.grid import UnRegularGridDataset +from numpy import zeros, arange, float +from netCDF4 import Dataset + +# for plotting +from matplotlib import pyplot as plt + +class RomsDataset(UnRegularGridDataset): + + __slots__ = list() + + def init_speed_coef(self, uname, vname): + # xi_u and eta_v must be specified because this dimension are not use in lon/lat + u = self.grid(uname, indexs=self.indexs) + v = self.grid(vname, indexs=self.indexs) + u = self.rho_2d(u.T).T + v = self.rho_2d(v) + self._speed_norm = (v ** 2 + u ** 2) ** .5 + + + @staticmethod + def rho_2d(x): + """Transformation to have u or v on same grid than h + """ + M, Lp = x.shape + new_x = zeros((M + 1, Lp)) + new_x[1:-1] = .5 * (x[:-1] + x[1:]) + new_x[0] = new_x[1] + new_x[-1] = new_x[-2] + return new_x + + @staticmethod + def psi2rho(var_psi): + """Transformation to have vrt on same grid than h + """ + M, L = var_psi.shape + Mp=M+1; Lp=L+1 + Mm=M-1; Lm=L-1 + var_rho = zeros((Mp, Lp)) + var_rho[1:M,1:L]=0.25*(var_psi[0:Mm,0:Lm]+var_psi[0:Mm,1:L]+var_psi[1:M,0:Lm]+var_psi[1:M,1:L]) + var_rho[0,:]=var_rho[1,:] + var_rho[Mp-1,:]=var_rho[M-1,:] + var_rho[:,0]=var_rho[:,1] + var_rho[:,Lp-1]=var_rho[:,L-1] + return var_rho + + + + +if __name__ == '__main__': + start_logger().setLevel('DEBUG') + + # Pick a depth + s_rho = 5 # 1000 m + depths = arange(-3500, 0, 500) # Make sure it is the same than used to generate the 'horizontal_section' files. + + depth = depths[s_rho] + + # Time loop + for time in range(7000, 7020): + + + infiletime = time%10 + filetime = time - time%10 + + # Using times: + grid_name = './GIGATL6/gigatl6_1h_horizontal_section.' + '{0:05}'.format(filetime) + '.nc' + + # or using dates + realyear_origin = datetime(2004,1,15) + date1 = realyear_origin + timedelta(days=float(time)/2.) + date2 = date1 + timedelta(days=float(4.5)) + + filedate = '{0:04}'.format(date1.year)+'-'+\ + '{0:02}'.format(date1.month) + '-'+\ + '{0:02}'.format(date1.day)+ '-'+\ + '{0:04}'.format(date2.year)+'-'+\ + '{0:02}'.format(date2.month) + '-' +\ + '{0:02}'.format(date2.day) + + print(filedate) + + lon_name, lat_name = 'lon', 'lat' + + + h = RomsDataset(grid_name, lon_name, lat_name, + indexs=dict(time=infiletime, + eta_rho=slice(500, 800), + xi_rho=slice(500, 800), + eta_v=slice(500, 800-1), + xi_u=slice(500, 800-1), + s_rho=s_rho) + ) + + # Must be set with time of grid + date = date1 + + # Identification every 2 mm + a, c = h.eddy_identification('ow', 'u', 'v', date, z_min = -10, z_max = -0.1, step = 0.05, pixel_limit=(10, 2000), shape_error=40, force_height_unit='m',force_speed_unit='m/s',vorticity_name='vrt') + + filename = 'gigatl6_1h_horizontal_section_' + '{0:04}'.format(-depth) + '_' + + with Dataset(date.strftime('Anticyclonic_' + filename + '%Y%m%d%H.nc'), 'w') as h: + a.to_netcdf(h) + with Dataset(date.strftime('Cyclonic_' + filename + '%Y%m%d%H.nc'), 'w') as h: + c.to_netcdf(h) + + # PLOT + + + fig = plt.figure(figsize=(15,7)) + ax = fig.add_axes([.03,.03,.94,.94]) + ax.set_title('Eddies detected -- Cyclonic(red) and Anticyclonic(blue)') + #ax.set_ylim(-75,75) + ax.set_xlim(250,360) + ax.set_aspect('equal') + a.display(ax, color='b', linewidth=.5) + c.display(ax, color='r', linewidth=.5) + ax.grid() + fig.savefig('eddies_' + date.strftime( filename + '%Y%m%d%H') +'.png') + + From b02273c0eb90a67a60ebf768574904bcbf9f09dd Mon Sep 17 00:00:00 2001 From: gula Date: Fri, 7 May 2021 23:25:39 +0200 Subject: [PATCH 12/20] update code --- examples/detection_gigatl6.py | 98 ++++++++++++++++++++--------- src/py_eddy_tracker/dataset/grid.py | 13 ++-- 2 files changed, 77 insertions(+), 34 deletions(-) diff --git a/examples/detection_gigatl6.py b/examples/detection_gigatl6.py index 7d151083..06880290 100644 --- a/examples/detection_gigatl6.py +++ b/examples/detection_gigatl6.py @@ -13,8 +13,11 @@ class RomsDataset(UnRegularGridDataset): def init_speed_coef(self, uname, vname): # xi_u and eta_v must be specified because this dimension are not use in lon/lat + print(self.indexs) u = self.grid(uname, indexs=self.indexs) v = self.grid(vname, indexs=self.indexs) + print('u.shape',u.shape) + u = self.rho_2d(u.T).T v = self.rho_2d(v) self._speed_norm = (v ** 2 + u ** 2) ** .5 @@ -46,59 +49,94 @@ def psi2rho(var_psi): var_rho[:,Lp-1]=var_rho[:,L-1] return var_rho +########################## +varname = 'zeta' +########################## if __name__ == '__main__': start_logger().setLevel('DEBUG') # Pick a depth - s_rho = 5 # 1000 m - depths = arange(-3500, 0, 500) # Make sure it is the same than used to generate the 'horizontal_section' files. - - depth = depths[s_rho] + + if varname == 'zeta': + s_rho = -1 + depth = 0 + else: + s_rho = 5 # 1000 m + depths = arange(-3500, 0, 500) # Make sure it is the same than used to generate the 'horizontal_section' files. + depth = depths[s_rho] # Time loop - for time in range(7000, 7020): + for time in range(34800, 34801): - - infiletime = time%10 - filetime = time - time%10 + dtfile = 120 + infiletime = time%dtfile + filetime = time - time%dtfile # Using times: - grid_name = './GIGATL6/gigatl6_1h_horizontal_section.' + '{0:05}'.format(filetime) + '.nc' + #grid_name = './GIGATL6/gigatl6_1h_horizontal_section.' + '{0:05}'.format(filetime) + '.nc' - # or using dates - realyear_origin = datetime(2004,1,15) - date1 = realyear_origin + timedelta(days=float(time)/2.) - date2 = date1 + timedelta(days=float(4.5)) - - filedate = '{0:04}'.format(date1.year)+'-'+\ - '{0:02}'.format(date1.month) + '-'+\ - '{0:02}'.format(date1.day)+ '-'+\ - '{0:04}'.format(date2.year)+'-'+\ - '{0:02}'.format(date2.month) + '-' +\ - '{0:02}'.format(date2.day) + if varname =='ow': + grid_name = './GIGATL6/gigatl6_1h_horizontal_section.' + '{0:05}'.format(filetime) + '.nc' + + + elif varname == 'zeta': + + # or using dates + realyear_origin = datetime(2004,1,15) + date1 = realyear_origin + timedelta(days=float(filetime)/24.) + date2 = date1 + timedelta(days=float(4.5)) + + filedate = '{0:04}'.format(date1.year)+'-'+\ + '{0:02}'.format(date1.month) + '-'+\ + '{0:02}'.format(date1.day)+ '-'+\ + '{0:04}'.format(date2.year)+'-'+\ + '{0:02}'.format(date2.month) + '-' +\ + '{0:02}'.format(date2.day) - print(filedate) - - lon_name, lat_name = 'lon', 'lat' - + print(filedate) + grid_name = './GIGATL6/GIGATL6_1h_inst_surf_' + filedate + '.nc' + + + # Identification + if varname=='zeta': + lon_name, lat_name = 'nav_lon_rho', 'nav_lat_rho' + elif varname=='ow': + lon_name, lat_name = 'lon', 'lat' + + # domain grid points: x1, x2, y1, y2 + x1, x2 = 0, 1500 + y1, y2 = 0, 2000 + h = RomsDataset(grid_name, lon_name, lat_name, indexs=dict(time=infiletime, - eta_rho=slice(500, 800), - xi_rho=slice(500, 800), - eta_v=slice(500, 800-1), - xi_u=slice(500, 800-1), + eta_rho=slice(y1, y2), + xi_rho=slice(x1, x2), + eta_v=slice(y1, y2-1), + xi_u=slice(x1, x2-1), + y_rho=slice(y1, y2), + x_rho=slice(x1, x2), + y_v=slice(y1, y2-1), + y_u=slice(y1, y2), + x_u=slice(x1, x2-1), + x_v=slice(x1, x2), s_rho=s_rho) ) # Must be set with time of grid date = date1 - # Identification every 2 mm - a, c = h.eddy_identification('ow', 'u', 'v', date, z_min = -10, z_max = -0.1, step = 0.05, pixel_limit=(10, 2000), shape_error=40, force_height_unit='m',force_speed_unit='m/s',vorticity_name='vrt') + # Identification + if varname=='zeta': + z_min = -2 ; z_max = 1; step = 0.02 + elif varname=='ow': + z_min = -10; z_max = -0.1; step = 0.05 + + + a, c = h.eddy_identification(varname, 'u', 'v', date, z_min = z_min, z_max = z_max, step = step, pixel_limit=(10, 2000), shape_error=40, force_height_unit='m',force_speed_unit='m/s',vorticity_name='vrt') filename = 'gigatl6_1h_horizontal_section_' + '{0:04}'.format(-depth) + '_' diff --git a/src/py_eddy_tracker/dataset/grid.py b/src/py_eddy_tracker/dataset/grid.py index 1d23fac0..3331667d 100644 --- a/src/py_eddy_tracker/dataset/grid.py +++ b/src/py_eddy_tracker/dataset/grid.py @@ -665,8 +665,11 @@ def eddy_identification( if precision is not None: precision /= factor - # Get ssh grid - data = 1e10*self.grid(grid_height).astype("f8") + # Get ssh or ow + if grid_height in ['ow']: + data = 1e10 * self.grid(grid_height).astype("f8") + else: + data = self.grid(grid_height).astype("f8") if grid_height in ['ow']: # Get vorticity as an aditional field (to identify cyc/acyc) @@ -706,9 +709,10 @@ def eddy_identification( ) z_min, z_max = z_min_p, z_max_p - print('step1 z_min, z_max, step',z_min, z_max, step) #debug + logger.warning('z_min, z_max, step are: %f, %f, %f.',z_min, z_max, step) #debug levels = arange(z_min - z_min % step, z_max - z_max % step + step, step) - + print('levels used are:',levels) + # Get x and y values x, y = self.x_c, self.y_c @@ -740,6 +744,7 @@ def eddy_identification( else: iterator = 1 if anticyclonic_search else -1 + # Loop over each collection for coll_ind, coll in enumerate(self.contours.iter(step=iterator)): corrected_coll_index = coll_ind From a2815c59f0d1664773d36d731e85d326635107e9 Mon Sep 17 00:00:00 2001 From: gula Date: Tue, 18 May 2021 21:54:40 +0200 Subject: [PATCH 13/20] add datarmor example --- examples/{ => datarmor}/detection_gigatl6.py | 97 +++++++++++++++----- examples/datarmor/run_detection_gigatl6.pbs | 47 ++++++++++ 2 files changed, 120 insertions(+), 24 deletions(-) rename examples/{ => datarmor}/detection_gigatl6.py (66%) create mode 100755 examples/datarmor/run_detection_gigatl6.pbs diff --git a/examples/detection_gigatl6.py b/examples/datarmor/detection_gigatl6.py similarity index 66% rename from examples/detection_gigatl6.py rename to examples/datarmor/detection_gigatl6.py index 06880290..54b108e5 100644 --- a/examples/detection_gigatl6.py +++ b/examples/datarmor/detection_gigatl6.py @@ -7,6 +7,9 @@ # for plotting from matplotlib import pyplot as plt +import sys,os,shutil + + class RomsDataset(UnRegularGridDataset): __slots__ = list() @@ -51,42 +54,82 @@ def psi2rho(var_psi): ########################## -varname = 'zeta' +varname = sys.argv[1] + +print("varname id %s", varname) + +#varname = 'zeta' ########################## if __name__ == '__main__': start_logger().setLevel('DEBUG') + + + + #### + # create case-specific folder + folder = varname + '/' + sys.argv[2] + '/' + + try: + os.mkdir(folder) + except OSError: + print ("Directory %s already exists" % folder) + + + # copy script in folder + shutil.copy('detection_gigatl6.py',folder) + - # Pick a depth + # Pick a depth/isopycnal if varname == 'zeta': s_rho = -1 depth = 0 - else: - s_rho = 5 # 1000 m - depths = arange(-3500, 0, 500) # Make sure it is the same than used to generate the 'horizontal_section' files. - depth = depths[s_rho] + elif varname == 'ow': + #isopycnal + grid_name = './iso/gigatl6_1h_isopycnal_section.01440.nc' + nc = Dataset(grid_name,'r') + isopycnals = nc.variables['isopycnal'][:] + nc.close() + + s_rho = 1 # + depth = isopycnals[s_rho] + + # Time loop - for time in range(34800, 34801): + #for time in range(34800, 34801): + for time in range(1440, 6000): + + + # hourly data + #dtfile = 1 + # 12-hourly + dtfile = 12 + + tfile = 5*24//dtfile + + ########### + realyear_origin = datetime(2004,1,15) + date = realyear_origin + timedelta(days=float(time)*dtfile/24.) + + infiletime = time%tfile + filetime = time - time%tfile - dtfile = 120 - infiletime = time%dtfile - filetime = time - time%dtfile # Using times: #grid_name = './GIGATL6/gigatl6_1h_horizontal_section.' + '{0:05}'.format(filetime) + '.nc' if varname =='ow': - grid_name = './GIGATL6/gigatl6_1h_horizontal_section.' + '{0:05}'.format(filetime) + '.nc' + #grid_name = './GIGATL6/gigatl6_1h_horizontal_section.' + '{0:05}'.format(filetime) + '.nc' + grid_name = './iso/gigatl6_1h_isopycnal_section.' + '{0:05}'.format(filetime) + '.nc' elif varname == 'zeta': # or using dates - realyear_origin = datetime(2004,1,15) - date1 = realyear_origin + timedelta(days=float(filetime)/24.) + date1 = realyear_origin + timedelta(days=float(filetime)*dtfile/24.) date2 = date1 + timedelta(days=float(4.5)) filedate = '{0:04}'.format(date1.year)+'-'+\ @@ -98,7 +141,7 @@ def psi2rho(var_psi): print(filedate) - grid_name = './GIGATL6/GIGATL6_1h_inst_surf_' + filedate + '.nc' + grid_name = './HIS/GIGATL6_1h_inst_surf_' + filedate + '.nc' # Identification @@ -108,11 +151,14 @@ def psi2rho(var_psi): lon_name, lat_name = 'lon', 'lat' # domain grid points: x1, x2, y1, y2 - x1, x2 = 0, 1500 + x1, x2 = 400, 1200 + y1, y2 = 400, 1200 + + x1, x2 = 0, 1600 y1, y2 = 0, 2000 h = RomsDataset(grid_name, lon_name, lat_name, - indexs=dict(time=infiletime, + indexs=dict(time=infiletime,time_counter=infiletime, eta_rho=slice(y1, y2), xi_rho=slice(x1, x2), eta_v=slice(y1, y2-1), @@ -126,23 +172,26 @@ def psi2rho(var_psi): s_rho=s_rho) ) - # Must be set with time of grid - date = date1 + # Identification if varname=='zeta': - z_min = -2 ; z_max = 1; step = 0.02 + z_min = -2 ; z_max = 1.5; step = 0.02 elif varname=='ow': - z_min = -10; z_max = -0.1; step = 0.05 + # ow is multiplied by 1e10 + z_min = -1; z_max = -0.01; step = 0.01 a, c = h.eddy_identification(varname, 'u', 'v', date, z_min = z_min, z_max = z_max, step = step, pixel_limit=(10, 2000), shape_error=40, force_height_unit='m',force_speed_unit='m/s',vorticity_name='vrt') - filename = 'gigatl6_1h_horizontal_section_' + '{0:04}'.format(-depth) + '_' + + #### + + filename = 'gigatl6_1h_' + varname + '_'+ '{0:04}'.format(depth) + '_' - with Dataset(date.strftime('Anticyclonic_' + filename + '%Y%m%d%H.nc'), 'w') as h: + with Dataset(date.strftime(folder + 'Anticyclonic_' + filename + '%Y%m%d%H.nc'), 'w') as h: a.to_netcdf(h) - with Dataset(date.strftime('Cyclonic_' + filename + '%Y%m%d%H.nc'), 'w') as h: + with Dataset(date.strftime(folder + 'Cyclonic_' + filename + '%Y%m%d%H.nc'), 'w') as h: c.to_netcdf(h) # PLOT @@ -157,6 +206,6 @@ def psi2rho(var_psi): a.display(ax, color='b', linewidth=.5) c.display(ax, color='r', linewidth=.5) ax.grid() - fig.savefig('eddies_' + date.strftime( filename + '%Y%m%d%H') +'.png') + fig.savefig(folder + 'eddies_' + date.strftime( filename + '%Y%m%d%H') +'.png') diff --git a/examples/datarmor/run_detection_gigatl6.pbs b/examples/datarmor/run_detection_gigatl6.pbs new file mode 100755 index 00000000..b969254d --- /dev/null +++ b/examples/datarmor/run_detection_gigatl6.pbs @@ -0,0 +1,47 @@ +#!/bin/bash +#PBS -q omp +#PBS -l ncpus=2 +#PBS -l mem=65g +#PBS -l walltime=24:00:00 +# COMMENT +# Tina ODAKA 06.03.2017 +# example for using datarmor with 'queue omp' in a personalised way +# By default, if you use 'omp' it use 1 node with 2 openmp threads +# here we show you example of it using 10 openmp threads (specify ncpus=10) +# and it use 10g of memory in this example (specifly mem=10g) +# and you can run max for 30 minutes with walltime setting +# To optimised multi thread, you need to place a command 'omplace' before your 'exe' +# Here, in this example, we suppose that your code 'exe' is compiled with +# intel compiler. +# for the users who use gnu compiler, plz use omp.10g.gnu.pbs +# +# cd to the directory you submitted your job + +cd $PBS_O_WORKDIR +qstat -f $PBS_JOBID + +##################### + + + +source /usr/share/Modules/3.2.10/init/bash +module load anaconda-py3.6/4.7.12 + +##################### + +#export PYTHONPATH=$HOME/Python_Modules_p3 +#setenv PYTHONPATH $HOME/Python_Modules_p3 +source activate eddy-track + + + +##################### +run='run00' +variable='ow' + +mkdir $variable/$run + +python detection_gigatl6.py $variable $run >& $variable/$run/output_detection_$run + + + From 88fcdc8a13018ab6d3108090dbb7ca8e99b216fa Mon Sep 17 00:00:00 2001 From: gula Date: Wed, 19 May 2021 18:11:22 +0200 Subject: [PATCH 14/20] update for gigatl1 --- examples/datarmor/detection_gigatl1.py | 234 +++++++++++++++++++++++++ src/py_eddy_tracker/dataset/grid.py | 10 +- 2 files changed, 242 insertions(+), 2 deletions(-) create mode 100644 examples/datarmor/detection_gigatl1.py diff --git a/examples/datarmor/detection_gigatl1.py b/examples/datarmor/detection_gigatl1.py new file mode 100644 index 00000000..2fec2b85 --- /dev/null +++ b/examples/datarmor/detection_gigatl1.py @@ -0,0 +1,234 @@ +from datetime import datetime, timedelta +from py_eddy_tracker import start_logger +from py_eddy_tracker.dataset.grid import UnRegularGridDataset +from numpy import zeros, arange, float +from netCDF4 import Dataset + +# for plotting +from matplotlib import pyplot as plt + +import sys,os,shutil + + +class RomsDataset(UnRegularGridDataset): + + __slots__ = list() + + def init_speed_coef(self, uname, vname): + # xi_u and eta_v must be specified because this dimension are not use in lon/lat + print(self.indexs) + u = self.grid(uname, indexs=self.indexs) + v = self.grid(vname, indexs=self.indexs) + print('u.shape',u.shape) + + u = self.rho_2d(u.T).T + v = self.rho_2d(v) + self._speed_norm = (v ** 2 + u ** 2) ** .5 + + + @staticmethod + def rho_2d(x): + """Transformation to have u or v on same grid than h + """ + M, Lp = x.shape + new_x = zeros((M + 1, Lp)) + new_x[1:-1] = .5 * (x[:-1] + x[1:]) + new_x[0] = new_x[1] + new_x[-1] = new_x[-2] + return new_x + + @staticmethod + def psi2rho(var_psi): + """Transformation to have vrt on same grid than h + """ + M, L = var_psi.shape + Mp=M+1; Lp=L+1 + Mm=M-1; Lm=L-1 + var_rho = zeros((Mp, Lp)) + var_rho[1:M,1:L]=0.25*(var_psi[0:Mm,0:Lm]+var_psi[0:Mm,1:L]+var_psi[1:M,0:Lm]+var_psi[1:M,1:L]) + var_rho[0,:]=var_rho[1,:] + var_rho[Mp-1,:]=var_rho[M-1,:] + var_rho[:,0]=var_rho[:,1] + var_rho[:,Lp-1]=var_rho[:,L-1] + return var_rho + +########################## + +varname = sys.argv[1] + +print("varname id %s", varname) + +#varname = 'zeta' + +simul = 'gigatl6' + +########################## + +if __name__ == '__main__': + start_logger().setLevel('DEBUG') + + + + #### + # create case-specific folder + folder = varname + '/' + sys.argv[2] + '/' + + try: + os.mkdir(folder) + except OSError: + print ("Directory %s already exists" % folder) + + + # copy script in folder + shutil.copy('detection_gigatl6.py',folder) + + + + #### + # define simulation + + if 'gigatl1' in simul: + folder_simulation = '~/megatl/GIGATL1/GIGATL1_1h/SURF/final/' + elif 'gigatl6' in simul: + folder_simulation = './HIS/' + + # Pick a depth/isopycnal + + if varname == 'zeta': + s_rho = -1 + depth = 0 + elif varname == 'ow': + #isopycnal + filename = './iso/gigatl6_1h_isopycnal_section.01440.nc' + nc = Dataset(filename,'r') + isopycnals = nc.variables['isopycnal'][:] + nc.close() + + s_rho = 1 # + depth = isopycnals[s_rho] + + + + # Time loop + #for time in range(34800, 34801): + for time in range(1440, 6000): + + + # hourly data + #dtfile = 1 + # 12-hourly + dtfile = 12 + + tfile = 5*24//dtfile + + ########### + realyear_origin = datetime(2004,1,15) + date = realyear_origin + timedelta(days=float(time)*dtfile/24.) + + infiletime = time%tfile + filetime = time - time%tfile + + + # Using times: + #filename = './GIGATL6/gigatl6_1h_horizontal_section.' + '{0:05}'.format(filetime) + '.nc' + + if varname =='ow': + #filename = './GIGATL6/gigatl6_1h_horizontal_section.' + '{0:05}'.format(filetime) + '.nc' + filename = './iso/gigatl6_1h_isopycnal_section.' + '{0:05}'.format(filetime) + '.nc' + + + elif varname == 'zeta': + + + # or using dates + date1 = realyear_origin + timedelta(days=float(filetime)*dtfile/24.) + date2 = date1 + timedelta(days=float(4.5)) + + if 'gigatl1' in simul: + filedate = '{0:04}'.format(date1.year)+'-'+\ + '{0:02}'.format(date1.month) + '-'+\ + '{0:02}'.format(date1.day) + + filename = folder_simulation + 'gigatl1_surf.' + filedate + '.nc' + gridname = '~/megatl/GIGATL1/gigatl1_grd.nc' + + elif 'gigatl6' in simul: + filedate = '{0:04}'.format(date1.year)+'-'+\ + '{0:02}'.format(date1.month) + '-'+\ + '{0:02}'.format(date1.day)+ '-'+\ + '{0:04}'.format(date2.year)+'-'+\ + '{0:02}'.format(date2.month) + '-' +\ + '{0:02}'.format(date2.day) + + print(filedate) + + filename = './HIS/GIGATL6_1h_inst_surf_' + filedate + '.nc' + gridname = filename + + # Identification + if varname=='zeta' and 'gigatl6' in simul: + lon_name, lat_name = 'nav_lon_rho', 'nav_lat_rho' + elif varname=='ow': + lon_name, lat_name = 'lon', 'lat' + else: + lon_name, lat_name = 'lon_rho', 'lat_rho' + + # domain grid points: x1, x2, y1, y2 + x1, x2 = 400, 1200 + y1, y2 = 400, 1200 + + x1, x2 = 0, 1600 + y1, y2 = 0, 2000 + + h = RomsDataset(filename, lon_name, lat_name, gridname = gridname, + indexs=dict(time=infiletime,time_counter=infiletime, + eta_rho=slice(y1, y2), + xi_rho=slice(x1, x2), + eta_v=slice(y1, y2-1), + xi_u=slice(x1, x2-1), + y_rho=slice(y1, y2), + x_rho=slice(x1, x2), + y_v=slice(y1, y2-1), + y_u=slice(y1, y2), + x_u=slice(x1, x2-1), + x_v=slice(x1, x2), + s_rho=s_rho) + ) + + + + # Identification + if varname=='zeta': + z_min = -2 ; z_max = 1.5; step = 0.02 + elif varname=='ow': + # ow is multiplied by 1e10 + z_min = -1; z_max = -0.01; step = 0.01 + + + a, c = h.eddy_identification(varname, 'u', 'v', date, z_min = z_min, z_max = z_max, step = step, pixel_limit=(10, 2000), shape_error=40, force_height_unit='m',force_speed_unit='m/s',vorticity_name='vrt') + + + #### + + filename = 'gigatl6_1h_' + varname + '_'+ '{0:04}'.format(depth) + '_' + + with Dataset(date.strftime(folder + 'Anticyclonic_' + filename + '%Y%m%d%H.nc'), 'w') as h: + a.to_netcdf(h) + with Dataset(date.strftime(folder + 'Cyclonic_' + filename + '%Y%m%d%H.nc'), 'w') as h: + c.to_netcdf(h) + + # PLOT + + + fig = plt.figure(figsize=(15,7)) + ax = fig.add_axes([.03,.03,.94,.94]) + ax.set_title('Eddies detected -- Cyclonic(red) and Anticyclonic(blue)') + #ax.set_ylim(-75,75) + ax.set_xlim(250,360) + ax.set_aspect('equal') + a.display(ax, color='b', linewidth=.5) + c.display(ax, color='r', linewidth=.5) + ax.grid() + fig.savefig(folder + 'eddies_' + date.strftime( filename + '%Y%m%d%H') +'.png') + + diff --git a/src/py_eddy_tracker/dataset/grid.py b/src/py_eddy_tracker/dataset/grid.py index 3331667d..49ccdd1a 100644 --- a/src/py_eddy_tracker/dataset/grid.py +++ b/src/py_eddy_tracker/dataset/grid.py @@ -255,6 +255,7 @@ class GridDataset(object): "y_dim", "coordinates", "filename", + "gridname", "dimensions", "indexs", "variables_description", @@ -270,12 +271,13 @@ class GridDataset(object): N = 1 def __init__( - self, filename, x_name, y_name, centered=None, indexs=None, unset=False + self, filename, x_name, y_name, gridname=None, centered=None, indexs=None, unset=False ): """ :param str filename: Filename to load :param str x_name: Name of longitude coordinates :param str y_name: Name of latitude coordinates + :param str gridname: Name of file containing lon,lat (filename is used if not defined) :param bool,None centered: Allow to know how coordinates could be used with pixel :param dict indexs: A dictionary that sets indexes to use for non-coordinate dimensions :param bool unset: Set to True to create an empty grid object without file @@ -292,6 +294,10 @@ def __init__( self.centered = centered self.contours = None self.filename = filename + if gridname is not None: + self.gridname = gridname + else: + self.gridname = filename self.coordinates = x_name, y_name self.vars = dict() self.indexs = dict() if indexs is None else indexs @@ -384,7 +390,7 @@ def load(self): Get coordinates and setup coordinates function """ x_name, y_name = self.coordinates - with Dataset(self.filename) as h: + with Dataset(self.gridname) as h: self.x_dim = h.variables[x_name].dimensions self.y_dim = h.variables[y_name].dimensions From 52422b899bd4222f4bbd42cd71b0bb898cc45d5c Mon Sep 17 00:00:00 2001 From: gula Date: Mon, 5 Jul 2021 18:36:37 +0200 Subject: [PATCH 15/20] add gaussian filtering in eddy_Detection --- src/py_eddy_tracker/dataset/grid.py | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/src/py_eddy_tracker/dataset/grid.py b/src/py_eddy_tracker/dataset/grid.py index 49ccdd1a..f6f40f87 100644 --- a/src/py_eddy_tracker/dataset/grid.py +++ b/src/py_eddy_tracker/dataset/grid.py @@ -613,6 +613,8 @@ def eddy_identification( force_speed_unit=None, vorticity_name='vrt', mle=1, + filtering=False, + filtering_scale=30, **kwargs, ): """ @@ -676,7 +678,7 @@ def eddy_identification( data = 1e10 * self.grid(grid_height).astype("f8") else: data = self.grid(grid_height).astype("f8") - + if grid_height in ['ow']: # Get vorticity as an aditional field (to identify cyc/acyc) vrt = self.grid(vorticity_name, indexs=self.indexs).astype('f8') @@ -687,9 +689,13 @@ def eddy_identification( data.mask = zeros(data.shape, dtype="bool") data.mask[isnan(data)] = 1 + if filtering and grid_height in ['zeta']: + data = data - gaussian_filter(data, filtering_scale) + # we remove noisy information if precision is not None: data = (data / precision).round() * precision + # Compute levels for ssh/okubo if z_min is None or z_max is None: if grid_height in ['ow']: From 32579cbb1203d1b7f100ae693ed8fa2a7476d949 Mon Sep 17 00:00:00 2001 From: gula Date: Mon, 5 Jul 2021 21:27:45 +0200 Subject: [PATCH 16/20] update from AntSimi version --- src/py_eddy_tracker/appli/eddies.py | 17 +++++++---- src/py_eddy_tracker/dataset/grid.py | 44 ++++++++++++++++++----------- src/py_eddy_tracker/eddy_feature.py | 2 +- 3 files changed, 41 insertions(+), 22 deletions(-) diff --git a/src/py_eddy_tracker/appli/eddies.py b/src/py_eddy_tracker/appli/eddies.py index 05c62b5c..e1b195ff 100644 --- a/src/py_eddy_tracker/appli/eddies.py +++ b/src/py_eddy_tracker/appli/eddies.py @@ -17,7 +17,7 @@ from numpy import dtype as npdtype from yaml import safe_load -from .. import EddyParser +from .. import EddyParser, identify_time from ..observations.observation import EddiesObservations, reverse_index from ..observations.tracking import TrackEddiesObservations from ..tracking import Correspondances @@ -165,7 +165,12 @@ def eddies_tracking(): parser.add_argument( "--zarr", action="store_true", help="Output will be wrote in zarr" ) - parser.add_argument("--unraw", action="store_true", help="Load unraw data") + parser.add_argument( + "--unraw", + action="store_true", + help="Load unraw data, use only for netcdf." + "If unraw is active, netcdf is loaded without apply scalefactor and add_offset.", + ) parser.add_argument( "--blank_period", type=int, @@ -354,11 +359,13 @@ def track( short_c = c._copy() short_c.shorter_than(size_max=nb_obs_min) - c.longer_than(size_min=nb_obs_min) - - long_track = c.merge(raw_data=raw) short_track = short_c.merge(raw_data=raw) + if c.longer_than(size_min=nb_obs_min) is False: + long_track = short_track.empty_dataset() + else: + long_track = c.merge(raw_data=raw) + # We flag obs if c.virtual: long_track["virtual"][:] = long_track["time"] == 0 diff --git a/src/py_eddy_tracker/dataset/grid.py b/src/py_eddy_tracker/dataset/grid.py index f6f40f87..2dad3bcd 100644 --- a/src/py_eddy_tracker/dataset/grid.py +++ b/src/py_eddy_tracker/dataset/grid.py @@ -821,6 +821,7 @@ def eddy_identification( # Here the considered contour passed shape_error test, masked_pixels test, # values strictly above (AEs) or below (CEs) the contour, number_pixels test) + # Compute amplitude reset_centroid, amp = self.get_amplitude( contour, @@ -832,7 +833,6 @@ def eddy_identification( mle=mle, **kwargs, ) - # If we have a valid amplitude if (not amp.within_amplitude_limits()) or (amp.amplitude == 0): contour.reject = 4 @@ -1003,7 +1003,6 @@ def get_uavg( all_contours.iter(start=level_start + step, step=step) ): level_contour = coll.get_nearest_path_bbox_contain_pt(centlon_e, centlat_e) - # Leave loop if no contours at level if level_contour is None: break @@ -1097,7 +1096,6 @@ class UnRegularGridDataset(GridDataset): def load(self): """Load variable (data)""" x_name, y_name = self.coordinates - with Dataset(self.filename) as h: self.x_dim = h.variables[x_name].dimensions self.y_dim = h.variables[y_name].dimensions @@ -1111,8 +1109,7 @@ def load(self): self.y_c = self.vars[y_name] self.init_pos_interpolator() - - + @property def bounds(self): """Give bounds""" @@ -2332,13 +2329,26 @@ def from_netcdf_cube(cls, filename, x_name, y_name, t_name, heigth=None): @classmethod def from_netcdf_list(cls, filenames, t, x_name, y_name, indexs=None, heigth=None): new = cls() - for i, t in enumerate(t): - d = RegularGridDataset(filenames[i], x_name, y_name, indexs=indexs) + for i, _t in enumerate(t): + filename = filenames[i] + logger.debug(f"load file {i:02d}/{len(t)} t={_t} : {filename}") + d = RegularGridDataset(filename, x_name, y_name, indexs=indexs) if heigth is not None: d.add_uv(heigth) - new.datasets.append((t, d)) + new.datasets.append((_t, d)) return new + def shift_files(self, t, filename, heigth=None, **rgd_kwargs): + """Add next file to the list and remove the oldest""" + + self.datasets = self.datasets[1:] + + d = RegularGridDataset(filename, **rgd_kwargs) + if heigth is not None: + d.add_uv(heigth) + self.datasets.append((t, d)) + logger.debug(f"shift and adding i={len(self.datasets)} t={t} : {filename}") + def interp(self, grid_name, t, lons, lats, method="bilinear"): """ Compute z over lons, lats @@ -2498,6 +2508,7 @@ def advect( else: mask_particule += isnan(x) + isnan(y) while True: + logger.debug(f"advect : t={t}") if (backward and t <= t1) or (not backward and t >= t1): t0, u0, v0, m0 = t1, u1, v1, m1 t1, d1 = generator.__next__() @@ -2524,25 +2535,21 @@ def advect( yield t, x, y def get_next_time_step(self, t_init): - first = True for i, (t, dataset) in enumerate(self.datasets): if t < t_init: continue - if first: - first = False - yield self.datasets[i - 1] + + logger.debug(f"i={i}, t={t}, dataset={dataset}") yield t, dataset def get_previous_time_step(self, t_init): - first = True i = len(self.datasets) for t, dataset in reversed(self.datasets): i -= 1 if t > t_init: continue - if first: - first = False - yield self.datasets[i + 1] + + logger.debug(f"i={i}, t={t}, dataset={dataset}") yield t, dataset @@ -2616,6 +2623,11 @@ def get_uv_quad(i0, j0, u, v, m, nb_x=0): i1, j1 = i0 + 1, j0 + 1 if nb_x != 0: i1 %= nb_x + i_max, j_max = m.shape + + if i1 >= i_max or j1 >= j_max: + return True, nan, nan, nan, nan, nan, nan, nan, nan + if m[i0, j0] or m[i0, j1] or m[i1, j0] or m[i1, j1]: return True, nan, nan, nan, nan, nan, nan, nan, nan # Extract value for u and v diff --git a/src/py_eddy_tracker/eddy_feature.py b/src/py_eddy_tracker/eddy_feature.py index 64855206..5c4760b3 100644 --- a/src/py_eddy_tracker/eddy_feature.py +++ b/src/py_eddy_tracker/eddy_feature.py @@ -65,7 +65,7 @@ def __init__( :param float contour_height: :param array data: :param float interval: - :param int mle: maximum number of local maxima in contour + :param int mle: maximum number of local extrema in contour :param int nb_step_min: number of intervals to consider an eddy :param int nb_step_to_be_mle: number of intervals to be considered as an another maxima """ From eef3eccbe1e3115120bd84326cce68b7c487cff5 Mon Sep 17 00:00:00 2001 From: gula Date: Thu, 15 Jul 2021 19:30:08 +0200 Subject: [PATCH 17/20] Update grid.py --- src/py_eddy_tracker/dataset/grid.py | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/src/py_eddy_tracker/dataset/grid.py b/src/py_eddy_tracker/dataset/grid.py index 2dad3bcd..ec691a7e 100644 --- a/src/py_eddy_tracker/dataset/grid.py +++ b/src/py_eddy_tracker/dataset/grid.py @@ -613,7 +613,7 @@ def eddy_identification( force_speed_unit=None, vorticity_name='vrt', mle=1, - filtering=False, + filtering=None, filtering_scale=30, **kwargs, ): @@ -689,8 +689,12 @@ def eddy_identification( data.mask = zeros(data.shape, dtype="bool") data.mask[isnan(data)] = 1 - if filtering and grid_height in ['zeta']: - data = data - gaussian_filter(data, filtering_scale) + if filtering is not None and grid_height in ['zeta']: + if filtering is 'highpass': + data = data - gaussian_filter(data, filtering_scale) + if filtering is 'lowpass': + data = data - gaussian_filter(data, filtering_scale) + # we remove noisy information if precision is not None: From d65270218c458772555cbc4e00fe9052367ae72e Mon Sep 17 00:00:00 2001 From: gula Date: Thu, 15 Jul 2021 19:32:41 +0200 Subject: [PATCH 18/20] Update grid.py --- src/py_eddy_tracker/dataset/grid.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/py_eddy_tracker/dataset/grid.py b/src/py_eddy_tracker/dataset/grid.py index ec691a7e..70a34d2b 100644 --- a/src/py_eddy_tracker/dataset/grid.py +++ b/src/py_eddy_tracker/dataset/grid.py @@ -693,7 +693,7 @@ def eddy_identification( if filtering is 'highpass': data = data - gaussian_filter(data, filtering_scale) if filtering is 'lowpass': - data = data - gaussian_filter(data, filtering_scale) + data = gaussian_filter(data, filtering_scale) # we remove noisy information From 842748a7888b708651d2a01ae76ff5985ac90ebe Mon Sep 17 00:00:00 2001 From: gula Date: Thu, 15 Jul 2021 19:46:03 +0200 Subject: [PATCH 19/20] Update grid.py --- src/py_eddy_tracker/dataset/grid.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/py_eddy_tracker/dataset/grid.py b/src/py_eddy_tracker/dataset/grid.py index 70a34d2b..c65cec30 100644 --- a/src/py_eddy_tracker/dataset/grid.py +++ b/src/py_eddy_tracker/dataset/grid.py @@ -692,8 +692,8 @@ def eddy_identification( if filtering is not None and grid_height in ['zeta']: if filtering is 'highpass': data = data - gaussian_filter(data, filtering_scale) - if filtering is 'lowpass': - data = gaussian_filter(data, filtering_scale) + elif filtering is 'lowpass': + data[:] = gaussian_filter(data, filtering_scale) # we remove noisy information From c690c552d884da59ad0865f1fd5afbc6ba56491d Mon Sep 17 00:00:00 2001 From: gula Date: Thu, 15 Jul 2021 20:53:52 +0200 Subject: [PATCH 20/20] Update grid.py correction of grid.py to allow for different mle values --- src/py_eddy_tracker/dataset/grid.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/py_eddy_tracker/dataset/grid.py b/src/py_eddy_tracker/dataset/grid.py index c65cec30..41a638f1 100644 --- a/src/py_eddy_tracker/dataset/grid.py +++ b/src/py_eddy_tracker/dataset/grid.py @@ -1073,6 +1073,8 @@ def get_amplitude( contour_height=contour_height, # All grid data=data, + # nb of allowed extrema + mle=mle, **kwargs, )