|
| 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