|  | 
|  | 1 | +# -*- coding: utf-8 -*- | 
|  | 2 | +""" | 
|  | 3 | +=========================================================================== | 
|  | 4 | +This file is part of py-eddy-tracker. | 
|  | 5 | +
 | 
|  | 6 | +    py-eddy-tracker is free software: you can redistribute it and/or modify | 
|  | 7 | +    it under the terms of the GNU General Public License as published by | 
|  | 8 | +    the Free Software Foundation, either version 3 of the License, or | 
|  | 9 | +    (at your option) any later version. | 
|  | 10 | +
 | 
|  | 11 | +    py-eddy-tracker is distributed in the hope that it will be useful, | 
|  | 12 | +    but WITHOUT ANY WARRANTY; without even the implied warranty of | 
|  | 13 | +    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the | 
|  | 14 | +    GNU General Public License for more details. | 
|  | 15 | +
 | 
|  | 16 | +    You should have received a copy of the GNU General Public License | 
|  | 17 | +    along with py-eddy-tracker.  If not, see <http://www.gnu.org/licenses/>. | 
|  | 18 | +
 | 
|  | 19 | +Copyright (c) 2014-2015 by Evan Mason | 
|  | 20 | + | 
|  | 21 | +=========================================================================== | 
|  | 22 | +
 | 
|  | 23 | +
 | 
|  | 24 | +py_eddy_tracker_amplitude.py | 
|  | 25 | +
 | 
|  | 26 | +Version 2.0.3 | 
|  | 27 | +
 | 
|  | 28 | +=========================================================================== | 
|  | 29 | +
 | 
|  | 30 | +""" | 
|  | 31 | +from py_eddy_tracker.observations import EddiesObservations, \ | 
|  | 32 | +    VirtualEddiesObservations, TrackEddiesObservations | 
|  | 33 | +from numpy import bool_, array, arange, ones, setdiff1d, zeros, uint16 | 
|  | 34 | +import logging | 
|  | 35 | + | 
|  | 36 | + | 
|  | 37 | +class Correspondances(list): | 
|  | 38 | +    """Object to store correspondances | 
|  | 39 | +    And run tracking | 
|  | 40 | +    """ | 
|  | 41 | +    UINT32_MAX = 4294967295 | 
|  | 42 | +    # Prolongation limit to 255 | 
|  | 43 | +    VIRTUAL_DTYPE = 'u1' | 
|  | 44 | +    # ID limit to 4294967295 | 
|  | 45 | +    ID_DTYPE = 'u4' | 
|  | 46 | +    # Track limit to 65535 | 
|  | 47 | +    N_DTYPE = 'u2' | 
|  | 48 | + | 
|  | 49 | + | 
|  | 50 | +    def __init__(self, datasets, virtual=0): | 
|  | 51 | +        """Initiate tracking | 
|  | 52 | +        """ | 
|  | 53 | +        # Correspondance dtype | 
|  | 54 | +        self.correspondance_dtype = [('in', 'u2'), | 
|  | 55 | +                                     ('out', 'u2'), | 
|  | 56 | +                                     ('id', self.ID_DTYPE)] | 
|  | 57 | +        # To count ID | 
|  | 58 | +        self.current_id = 0 | 
|  | 59 | +        # Dataset to iterate | 
|  | 60 | +        self.datasets = datasets | 
|  | 61 | +        self.previous2_obs = None | 
|  | 62 | +        self.previous_obs = None | 
|  | 63 | +        self.current_obs = EddiesObservations.load_from_netcdf( | 
|  | 64 | +            self.datasets[0]) | 
|  | 65 | +         | 
|  | 66 | +        # To use virtual obs | 
|  | 67 | +        # Number of obs which can prolongate real observations | 
|  | 68 | +        self.nb_virtual = virtual | 
|  | 69 | +        # Activation or not | 
|  | 70 | +        self.virtual = virtual > 0 | 
|  | 71 | +        self.virtual_obs = None | 
|  | 72 | +        self.previous_virtual_obs = None | 
|  | 73 | +        if self.virtual: | 
|  | 74 | +            # Add field to dtype to follow virtual observations | 
|  | 75 | +            self.correspondance_dtype += [ | 
|  | 76 | +                ('virtual', bool_), | 
|  | 77 | +                ('virtual_length', self.VIRTUAL_DTYPE)] | 
|  | 78 | +         | 
|  | 79 | +        # Array to simply merged | 
|  | 80 | +        self.nb_obs_by_tracks = None | 
|  | 81 | +        self.i_current_by_tracks = None | 
|  | 82 | +        self.nb_obs = 0 | 
|  | 83 | +        self.eddies = None | 
|  | 84 | + | 
|  | 85 | +    def swap_dataset(self, dataset): | 
|  | 86 | +        """ | 
|  | 87 | +        """ | 
|  | 88 | +        self.previous2_obs = self.previous_obs | 
|  | 89 | +        self.previous_obs = self.current_obs | 
|  | 90 | +        self.current_obs = EddiesObservations.load_from_netcdf(dataset) | 
|  | 91 | +         | 
|  | 92 | +    def store_correspondance(self, i_previous, i_current): | 
|  | 93 | +        correspondance = array( | 
|  | 94 | +            i_previous, | 
|  | 95 | +            dtype=self.correspondance_dtype) | 
|  | 96 | +        correspondance['out'] = i_current | 
|  | 97 | + | 
|  | 98 | +        if self.virtual: | 
|  | 99 | +            correspondance['virtual'] = i_previous >= nb_real_obs | 
|  | 100 | +     | 
|  | 101 | +    def id_generator(self, nb_id): | 
|  | 102 | +        """Generation id and incrementation | 
|  | 103 | +        """ | 
|  | 104 | +        values = arange(self.current_id, self.current_id + nb_id) | 
|  | 105 | +        self.current_id += nb_id | 
|  | 106 | +        return values | 
|  | 107 | +     | 
|  | 108 | +    def recense_dead_id_to_extend(self): | 
|  | 109 | +        """Recense dead id to extend in virtual observation | 
|  | 110 | +        """ | 
|  | 111 | +        nb_virtual_extend = 0 | 
|  | 112 | +        # List previous id which are not use in the next step | 
|  | 113 | +        dead_id = setdiff1d(self[-2]['id'], self[-1]['id']) | 
|  | 114 | +        nb_dead = dead_id.shape[0] | 
|  | 115 | +        logging.debug('%d death of real obs in this step', nb_dead) | 
|  | 116 | +        if not self.virtual: | 
|  | 117 | +            return | 
|  | 118 | +        # Creation of an virtual step for dead one | 
|  | 119 | +        self.virtual_obs = VirtualEddiesObservations(size=nb_dead + nb_virtual_extend) | 
|  | 120 | + | 
|  | 121 | +        # Find mask/index on previous correspondance to extrapolate | 
|  | 122 | +        # position | 
|  | 123 | +        list_previous_id = self[-2]['id'].tolist() | 
|  | 124 | +        i_dead_id = [list_previous_id.index(i) for i in dead_id] | 
|  | 125 | + | 
|  | 126 | +        # Selection of observations on N-2 and N-1 | 
|  | 127 | +        obs_a = self.previous2_obs.obs[self[-2][i_dead_id]['in']] | 
|  | 128 | +        obs_b = self.previous_obs.obs[self[-2][i_dead_id]['out']] | 
|  | 129 | +        # Position N-2 : A | 
|  | 130 | +        # Position N-1 : B | 
|  | 131 | +        # Virtual Position : C | 
|  | 132 | +        # New position C = B + AB | 
|  | 133 | +        self.virtual_obs['dlon'][:nb_dead] = obs_b['lon'] - obs_a['lon'] | 
|  | 134 | +        self.virtual_obs['dlat'][:nb_dead] = obs_b['lat'] - obs_a['lat'] | 
|  | 135 | +        self.virtual_obs['lon'][:nb_dead | 
|  | 136 | +            ] = obs_b['lon'] + self.virtual_obs['dlon'][:nb_dead] | 
|  | 137 | +        self.virtual_obs['lat'][:nb_dead | 
|  | 138 | +            ] = obs_b['lat'] + self.virtual_obs['dlat'][:nb_dead] | 
|  | 139 | +        # Id which are extended | 
|  | 140 | +        self.virtual_obs['track'][:nb_dead] = dead_id | 
|  | 141 | +        # Add previous virtual | 
|  | 142 | +        if nb_virtual_extend > 0: | 
|  | 143 | +            obs_to_extend = self.previous_virtual_obs.obs[i_virtual_dead_id][alive_virtual_obs] | 
|  | 144 | +            self.virtual_obs['lon'][nb_dead:] = obs_to_extend['lon'] + obs_to_extend['dlon'] | 
|  | 145 | +            self.virtual_obs['lat'][nb_dead:] = obs_to_extend['lat'] + obs_to_extend['dlat'] | 
|  | 146 | +            self.virtual_obs['track'][nb_dead:] = obs_to_extend['track'] | 
|  | 147 | +            self.virtual_obs['segment_size'][nb_dead:] = obs_to_extend['segment_size'] | 
|  | 148 | +        # Count | 
|  | 149 | +        self.virtual_obs['segment_size'][:] += 1 | 
|  | 150 | +     | 
|  | 151 | +    def track(self): | 
|  | 152 | +        """Run tracking | 
|  | 153 | +        """ | 
|  | 154 | +        START = True | 
|  | 155 | +        FLG_VIRTUAL = False | 
|  | 156 | +        # We begin with second file, first one is in previous | 
|  | 157 | +        for file_name in self.datasets[1:]: | 
|  | 158 | +            self.swap_dataset(file_name) | 
|  | 159 | +            logging.debug('%s match with previous state', file_name) | 
|  | 160 | +            logging.debug('%d obs to match', len(self.current_obs)) | 
|  | 161 | +         | 
|  | 162 | +            nb_real_obs = len(self.previous_obs) | 
|  | 163 | +            if FLG_VIRTUAL: | 
|  | 164 | +                logging.debug('%d virtual obs will be add to previous', | 
|  | 165 | +                              len(self.virtual_obs)) | 
|  | 166 | +                # If you comment this the virtual fonctionnality will be disable | 
|  | 167 | +                self.previous_obs = self.previous_obs.merge(self.virtual_obs) | 
|  | 168 | +         | 
|  | 169 | +            i_previous, i_current = self.previous_obs.tracking(self.current_obs) | 
|  | 170 | +            nb_match = i_previous.shape[0] | 
|  | 171 | + | 
|  | 172 | +            #~ self.store_correspondance(i_previous, i_current) | 
|  | 173 | +            correspondance = array(i_previous, dtype=self.correspondance_dtype) | 
|  | 174 | +            correspondance['out'] = i_current | 
|  | 175 | +             | 
|  | 176 | +            if self.virtual: | 
|  | 177 | +                correspondance['virtual'] = i_previous >= nb_real_obs | 
|  | 178 | +                 | 
|  | 179 | +            if START: | 
|  | 180 | +                START = False | 
|  | 181 | +                # Set an id for each match | 
|  | 182 | +                correspondance['id'] = self.id_generator(nb_match) | 
|  | 183 | +                self.append(correspondance) | 
|  | 184 | +                continue | 
|  | 185 | + | 
|  | 186 | +            # We set all id to UINT32_MAX | 
|  | 187 | +            id_previous = ones(len(self.previous_obs), dtype=self.ID_DTYPE) * self.UINT32_MAX | 
|  | 188 | +            # We get old id for previously eddies tracked | 
|  | 189 | +            previous_id = self[-1]['id'] | 
|  | 190 | +            id_previous[self[-1]['out']] = previous_id | 
|  | 191 | +            correspondance['id'] = id_previous[correspondance['in']] | 
|  | 192 | + | 
|  | 193 | +            if FLG_VIRTUAL: | 
|  | 194 | +                nb_rebirth = correspondance['virtual'].sum() | 
|  | 195 | +                if nb_rebirth != 0: | 
|  | 196 | +                    logging.debug('%d re-birth due to prolongation with' | 
|  | 197 | +                                  ' virtual observations', nb_rebirth) | 
|  | 198 | +                    # Set id for virtual | 
|  | 199 | +                    i_virtual = correspondance['in'][correspondance['virtual']] - nb_real_obs | 
|  | 200 | +                    correspondance['id'][correspondance['virtual']] = \ | 
|  | 201 | +                        self.virtual_obs['track'][i_virtual] | 
|  | 202 | +                    correspondance['virtual_length'][correspondance['virtual']] = \ | 
|  | 203 | +                        self.virtual_obs['segment_size'][i_virtual] | 
|  | 204 | +             | 
|  | 205 | +            # new_id is equal to UINT32_MAX we must add a new ones | 
|  | 206 | +            # we count the number of new | 
|  | 207 | +            mask_new_id = correspondance['id'] == UINT32_MAX | 
|  | 208 | +            nb_new_tracks = mask_new_id.sum() | 
|  | 209 | +            logging.debug('%d birth in this step', nb_new_tracks) | 
|  | 210 | +            # Set new id | 
|  | 211 | +            correspondance['id'][mask_new_id] = self.id_generator(nb_new_tracks) | 
|  | 212 | + | 
|  | 213 | +            self.append(correspondance) | 
|  | 214 | + | 
|  | 215 | +            # SECTION for virtual observation | 
|  | 216 | +            nb_virtual_prolongate = 0 | 
|  | 217 | +            if FLG_VIRTUAL: | 
|  | 218 | +                # Save previous state to count virtual obs | 
|  | 219 | +                self.previous_virtual_obs = self.virtual_obs | 
|  | 220 | +                virtual_dead_id = setdiff1d(self.virtual_obs['track'], | 
|  | 221 | +                                            correspondance['id']) | 
|  | 222 | +                list_previous_virtual_id = self.virtual_obs['track'].tolist() | 
|  | 223 | +                i_virtual_dead_id = [ | 
|  | 224 | +                    list_previous_virtual_id.index(i) for i in virtual_dead_id] | 
|  | 225 | +                # Virtual obs which can be prolongate | 
|  | 226 | +                alive_virtual_obs = self.virtual_obs['segment_size'][i_virtual_dead_id] < self.nb_virtual | 
|  | 227 | +                nb_virtual_prolongate = alive_virtual_obs.sum() | 
|  | 228 | +                logging.debug('%d virtual obs will be prolongate on the ' | 
|  | 229 | +                              'next step', nb_virtual_prolongate) | 
|  | 230 | + | 
|  | 231 | +            self.recense_dead_id_to_extend() | 
|  | 232 | + | 
|  | 233 | +            if self.virtual: | 
|  | 234 | +                FLG_VIRTUAL = True | 
|  | 235 | + | 
|  | 236 | +    def prepare_merging(self): | 
|  | 237 | +        # count obs by tracks (we add directly one, because correspondance | 
|  | 238 | +        # is an interval) | 
|  | 239 | +        self.nb_obs_by_tracks = zeros(self.current_id, dtype=self.N_DTYPE) + 1 | 
|  | 240 | +        for correspondance in self: | 
|  | 241 | +            self.nb_obs_by_tracks[correspondance['id']] += 1 | 
|  | 242 | +            if self.virtual: | 
|  | 243 | +                # When start is virtual, we don't have a previous correspondance | 
|  | 244 | +                self.nb_obs_by_tracks[correspondance['id'][correspondance['virtual']] | 
|  | 245 | +                                 ] += correspondance['virtual_length'][correspondance['virtual']] | 
|  | 246 | + | 
|  | 247 | +        # Compute index of each tracks | 
|  | 248 | +        self.i_current_by_tracks = self.nb_obs_by_tracks.cumsum() - self.nb_obs_by_tracks | 
|  | 249 | +        # Number of global obs | 
|  | 250 | +        self.nb_obs = nb_obs_by_tracks.sum() | 
|  | 251 | +        logging.info('%d tracks identified', self.current_id) | 
|  | 252 | +        logging.info('%d observations will be join', self.nb_obs) | 
|  | 253 | + | 
|  | 254 | +    def merge(self): | 
|  | 255 | +        # Start create netcdf to agglomerate all eddy | 
|  | 256 | +        self.eddies = TrackEddiesObservations(size=self.nb_obs) | 
|  | 257 | + | 
|  | 258 | +        # Calculate the index in each tracks, we compute in u4 and translate | 
|  | 259 | +        # in u2 (which are limited to 65535) | 
|  | 260 | +        logging.debug('Compute global index array (N)') | 
|  | 261 | +        n = arange(nb_obs, | 
|  | 262 | +                   dtype='u4') - self.i_current_by_tracks.repeat(self.nb_obs_by_tracks) | 
|  | 263 | +        self.eddies['n'][:] = uint16(n) | 
|  | 264 | +        logging.debug('Compute global track array') | 
|  | 265 | +        self.eddies['track'][:] = arange(self.current_id).repeat(self.nb_obs_by_tracks) | 
|  | 266 | + | 
|  | 267 | +        # Start loading identification again to save in the finals tracks | 
|  | 268 | +        # Load first file | 
|  | 269 | +        eddies_previous = EddiesObservations.load_from_netcdf(self.datasets[0]) | 
|  | 270 | +        # Set type of eddy with first file | 
|  | 271 | +        self.eddies.sign_type = eddies_previous.sign_type | 
|  | 272 | + | 
|  | 273 | +        # To know if the track start | 
|  | 274 | +        first_obs_save_in_tracks = zeros(i_current_by_tracks.shape, | 
|  | 275 | +                                            dtype=bool_) | 
|  | 276 | + | 
|  | 277 | +        for i, file_name in enumerate(FILENAMES[1:]): | 
|  | 278 | +            # Load current file (we begin with second one) | 
|  | 279 | +            self.current_obs = EddiesObservations.load_from_netcdf(file_name) | 
|  | 280 | +            # We select the list of id which are involve in the correspondance | 
|  | 281 | +            i_id = self[i]['id'] | 
|  | 282 | +            # Index where we will write in the final object | 
|  | 283 | +            index_final = i_current_by_tracks[i_id] | 
|  | 284 | + | 
|  | 285 | +            # First obs of eddies | 
|  | 286 | +            m_first_obs = -first_obs_save_in_tracks[i_id] | 
|  | 287 | +            if m_first_obs.any(): | 
|  | 288 | +                # Index in the current file | 
|  | 289 | +                index_in = self[i]['in'][m_first_obs] | 
|  | 290 | +                # Copy all variable | 
|  | 291 | +                for var, _ in eddies_current.obs.dtype.descr: | 
|  | 292 | +                    self.eddies[var][index_final[m_first_obs] | 
|  | 293 | +                        ] = eddies_previous[var][index_in] | 
|  | 294 | +                # Increment | 
|  | 295 | +                i_current_by_tracks[i_id[m_first_obs]] += 1 | 
|  | 296 | +                # Active this flag, we have only one first by tracks | 
|  | 297 | +                first_obs_save_in_tracks[i_id] = True | 
|  | 298 | +                index_final = i_current_by_tracks[i_id] | 
|  | 299 | + | 
|  | 300 | +            # Index in the current file | 
|  | 301 | +            index_current = self[i]['out'] | 
|  | 302 | +             | 
|  | 303 | +            if self.virtual: | 
|  | 304 | +                # If the flag virtual in correspondance is active, | 
|  | 305 | +                # the previous is virtual | 
|  | 306 | +                m_virtual = self[i]['virtual'] | 
|  | 307 | +                if m_virtual.any(): | 
|  | 308 | +                    index_virtual = index_final[m_virtual] | 
|  | 309 | +                    # Incrementing index | 
|  | 310 | +                    i_current_by_tracks[i_id[m_virtual] | 
|  | 311 | +                        ] += self[i]['virtual_length'][m_virtual] | 
|  | 312 | +                    # Get new index | 
|  | 313 | +                    index_final = i_current_by_tracks[i_id] | 
|  | 314 | + | 
|  | 315 | +            # Copy all variable | 
|  | 316 | +            for var, _ in eddies_current.obs.dtype.descr: | 
|  | 317 | +                self.eddies[var][index_final | 
|  | 318 | +                    ] = eddies_current[var][index_current] | 
|  | 319 | + | 
|  | 320 | +            # Add increment for each index used | 
|  | 321 | +            i_current_by_tracks[i_id] += 1 | 
|  | 322 | +            eddies_previous = eddies_current | 
|  | 323 | + | 
0 commit comments