Collocate In-situ observations

[acapet]

Hello,

I have a data frame (df) for in-situ profiles with date, lon, lat, and other information.

My current objective is to add columns to this data frame, on the basis of a set of eddy tracks.

• Is the float within an eddy contour? YES/NO
• Lon, Lat for the eddy center,
• Radius and Age (lifetime) of this eddy.

I'm not entirely sure how to address that.

My first approach was to loop on days, subset both the in situ data frame and the eddies collection, and then use the inside function. With the following code, I can determine if a given profile is within any eddy identified for that day.

ta1 = TrackEddiesObservations.load_file("../out_"+set1+"/Tracks_Overlap/Anticyclonic.nc")
dmin=ta1.time.min()
dmax=ta1.time.max()

for d in range(dmin,dmax):
    a1=ta1.extract_with_period((d,d))
    dfsub = df[df["time"]==d]
    for dl in dfsub.iterrows():
        isin = a1.inside(x=np.array([dl[1].longitude], dtype=np.float64), y=np.array([dl[1].latitude], dtype=np.float64))
        if isin:
                ....

However I have difficulties to get the index of the eddy for wich the condition is true, and then access the relevant info.

Should I loop over all eddies in ´a1´ ?
How can this be done?

I guess I just need a few general guidance, in particular regarding indexing through the TrackEddiesObservations objects.

PS: I guess looping on days isn't strictly necessary here, but it was also done for visualization reasons.

Thanks !

[AntSimi]

Inside method return for each given position if there are in one eddy of the list (you get an array of boolean).
So you could imagine use like this:

mask = a1.inside(x=df.longitude.values.astype('f8'), y=df.latitude.values.astype('f8'))
obs_inside = df[mask]     

To get parameters of eddy which enclosed Insitu observation, use contains method which will return index of a1 observations (-1 in case of insitu obs are not enclosed by eddy) instead of a boolean:

index_eddy = a1.contains(x=df.longitude.values.astype('f8'), y=df.latitude.values.astype('f8'))
mask = index_eddy != -1
insitu_inside = df[mask] 
eddy_lat_which_enclosed_insitu = a1.lat[index_eddy[mask]]

[acapet]

Thanks! Somehow, the last line in the following induces a segmentation fault, and I can't find the dumped core to give additional information..

for d in range(dmin,dmax):
    a1=ta1.extract_with_period((d,d))
    dfsub = df[df["time"]==d]
    mask = a1.inside(x=dfsub.longitude.values.astype('f8'), y=dfsub.latitude.values.astype('f8'))

    if any(mask):
        print(d)
        index_eddy = a1.contains(x=dfsub.longitude.values.astype('f8'), y=dfsub.latitude.values.astype('f8'))

[AntSimi]

All your longitude and latitude are valid? Could you share min/max? You could also skip mask computation and use directly contains to avoid double computation.

[acapet]

They are valid, but generally contain only 1 element.

    print(dfsub.longitude.values.astype('f8'))
    print(dfsub.latitude.values.astype('f8'))

works fine below the if part above

[AntSimi]

Which version of py eddy tracker did you use?

[acapet]

I updated everything (pull, re-install) two days ago, and used the overlap tracker (but that was before, though, the tracking files were produced some months ago)

[CoriPegliasco]

Hi,
I also have an issue with contains (on the latest version)

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-33-d13a49993c9c> in <module>
      7                             dissociate_network=True,
      8                             correct_close_events=45,
----> 9                             remove_dead_end=14
     10                             )

--/dev/py-eddy-tracker/src/py_eddy_tracker/observations/network.py in analysis_coherence(self, date_function, uv_params, advection_mode, dt_advect, step_mesh, output_name, dissociate_network, correct_close_events, remove_dead_end)
   1368             output_name=output_name,
   1369             dt_advect=dt_advect,
-> 1370             step_mesh=step_mesh,
   1371         )
   1372 

--/dev/py-eddy-tracker/src/py_eddy_tracker/observations/network.py in segment_coherence(self, date_function, uv_params, advection_mode, dt_advect, step_mesh, output_name)
   1496                         i_target_b,
   1497                         pct_target_b,
-> 1498                         n_days=-dt_advect,
   1499                     )
   1500 

--/dev/py-eddy-tracker/src/py_eddy_tracker/observations/groups.py in particle_candidate(x, y, c, eddies, t_start, i_target, pct, **kwargs)
    121     translate_end = where(m_end)[0]
    122     # Id eddies for each alive particle (in core and extern)
--> 123     i_end = e_end.contains(x, y)
    124     # compute matrix and fill target array
    125     get_matrix(i_start, i_end, translate_start, translate_end, i_target, pct)

--/dev/py-eddy-tracker/src/py_eddy_tracker/observations/observation.py in contains(self, x, y, intern)
   2057         m = ~(isnan(x) + isnan(y))
   2058         i = -ones(x.shape, dtype="i4")
-> 2059         i[m] = poly_indexs(x[m], y[m], self[xname], self[yname])
   2060         return i
   2061 

~/.conda/envs/eddy37/lib/python3.7/site-packages/numba/np/arraymath.py in array_min_impl()
    533         def array_min_impl(arry):
    534             if arry.size == 0:
--> 535                 raise ValueError(MSG)
    536 
    537             it = np.nditer(arry)

ValueError: zero-size array to reduction operation minimum which has no identity

Could it be because I have longitudes between [-180°, 180°] and poly_indexs is waiting for [0°, 360°]?
Thus the validity for longitudes is ok if >0 ?

[AntSimi]

ok @acapet, i do some update to speed up poly_index few month ago, maybe it's a side effect ..., to understand the problem i will need a Minimal Working Example, if you are ok, you could create like that:

import pickle
a1.write_file(filename="eddy.nc")
x, y = dfsub.longitude.values.astype('f8'), dfsub.latitude.values.astype('f8')
with open("insitu_position.pkl", "wb") as pkl:
    pickle.dump((x,y), pkl)
a1.contains(x=x, y=y)

[CoriPegliasco]

The error and associated files

error_contains.tar.gz

from py_eddy_tracker.observations.network import NetworkObservations
import pickle
n = NetworkObservations.load_file("obs_error.nc")
with open("test_particles", "rb") as pkl:
           x,y= pickle.load(pkl)
i_end = n.contains(x, y)

[AntSimi]

@CoriPegliasco in your case x & y are empty array...

[acapet]

Here are the files for the minimal example :

ProfileCollocationFail.tar.gz

I have to mention that the fail is not systematic, which makes it even stranger.
Executing the following goes fine sometimes, but fails at some point if I repeat ~5-10 times.

from py_eddy_tracker.observations.network import NetworkObservations
import pickle
n = NetworkObservations.load_file("./ProfileCollocationFail/eddy.nc")
with open("./ProfileCollocationFail/insitu_position.pkl", "rb") as pkl:
       x,y= pickle.load(pkl)
i_end = n.contains(x, y)

[AntSimi]

So, it's due to a specific case, which give a bad index and sometime this index is in a forbidden area...
I didn't notice this wrong indexing/code due too a numba silence (numba/numba#7045).
I will release a patch soon.

[acapet]

Ok. If I got this right, the problem is due to applying the method when applied on small sets (with one eddy, or one profile).
This arises from the fact that I loop on individual days, which is necessary for matching eddies and external data since "contains" only considers lon and lat, and not time.

I guess the patch will totally do for me. If that goes further, one could of a version of "contain" including the time dimension.

[AntSimi]

This arises from the fact that I loop on individual days, which is necessary for matching eddies and external data since "contains" only considers lon and lat, and not time

It come when you call contains with one or two observations to compare with eddies.

Root problem come from a bad code which was silent. So i will re-wrote it.
And this root problem is in a method which are use in several other method...

If that goes further, one could of a version of "contain" including the time dimension.

I am not pro to do that, because you must manage time sampling which could be different between your eddies and other observations. My opinion is : it's best that user manage explicitly time coherence.
Maybe, it could be build with iter_on method.