-
Notifications
You must be signed in to change notification settings - Fork 48
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Using EddyTracking to compare observations and model outputs #7
Comments
Good try, but tracking was not think for that... In the example below i used a cost function based on contour overlap: With a ratio > 30 % Very little eddies in big ones will be not selected with this cost function import pylook
from matplotlib import pyplot as plt
from py_eddy_tracker.observations.observation import EddiesObservations
from numpy import linspace, ones, where
def reverse_index(index, nb):
m = ones(nb, dtype="bool")
m[index] = False
return where(m)[0]
a_filt = EddiesObservations.load_file("Anticyclonic_20200424.nc")
a_raw = EddiesObservations.load_file("Anticyclonic_20200424_raw.nc")
# Get indexs of close observation (based on overlap contour)
i, j, cost = a_filt.match(a_raw)
m = cost > 0.3
i, j = i[m], j[m]
# Get index not used
i_, j_ = reverse_index(i, len(a_filt)), reverse_index(j, len(a_raw))
# Subset
a_filt_junk = a_filt.index(i_)
a_raw_junk = a_raw.index(j_)
a_filt = a_filt.index(i)
a_raw = a_raw.index(j)
# Plot
ax = plt.subplot(111, projection="plat_carre")
ax.grid()
ax.set_xlim(0, 360)
ax.set_ylim(-80, 80)
kwargs_display = dict(lw=1, extern_only=True)
a_filt.display(ax, label="Anticyclonic unfiltered", color="r", **kwargs_display)
a_raw.display(ax, label="Anticyclonic", color="blue", **kwargs_display)
kwargs_display["lw"] = 0.25
a_filt_junk.display(
ax, label="Anticyclonic unfiltered not associate", color="r", **kwargs_display
)
a_raw_junk.display(ax, label="Anticyclonic not associate", color="b", **kwargs_display)
ax.legend(loc="upper right")
plt.show() If you want save observations in a file use "to_netcdf" method |
Thank you for your useful proposition :) However, since the eddies may be a bit displaced from obs to model, a criteria only on area is too restrictive for us (see figures). I worked a bit from your solution and the tracking function in EddiesObservations to develop another solution that I will copy/paste here. import cartopy.crs as ccrs
from matplotlib import pyplot as plt
from py_eddy_tracker.observations.observation import EddiesObservations
from py_eddy_tracker import generic
from numpy import linspace, ones, where, arange
root = 'data\\eddies-obs-model-tocompare\\'
AVISO_AC_file = root+'Anticyclonic_obs_20200424.nc'
GIOPS_AC_file = root+'Anticyclonic_20200424.nc'
def reverse_index(index, nb):
m = ones(nb, dtype="bool")
m[index] = False
return where(m)[0]
AVISO_AC = EddiesObservations.load_file(AVISO_AC_file)
GIOPS_AC = EddiesObservations.load_file(GIOPS_AC_file)
i_aviso, i_giops, cost_mat = AVISO_AC.tracking(GIOPS_AC)
# Get index not used
i_aviso_junk, i_giops_junk = reverse_index(i_aviso, len(AVISO_AC)), reverse_index(i_giops, len(GIOPS_AC))
# Subset
AVISO_AC_junk = AVISO_AC.index(i_aviso_junk)
GIOPS_AC_junk = GIOPS_AC.index(i_giops_junk)
AVISO_AC = AVISO_AC.index(i_aviso)
GIOPS_AC = GIOPS_AC.index(i_giops)
# Plot
plt.figure(figsize=(15,10))
data_crs=ccrs.PlateCarree()
ax = plt.axes(projection=data_crs)
ax.coastlines()
ax.set_extent([-80,-45,30,50])
ax.set_aspect('auto')
ax.grid()
kwargs_display = dict(lw=1.25, extern_only=True)
AVISO_AC.display(ax, label="Anticyclonic AVISO", color="black",transform=data_crs, **kwargs_display)
GIOPS_AC.display(ax, label="Anticyclonic GIOPS", color="red",transform=data_crs, **kwargs_display)
kwargs_display["lw"] = 1
AVISO_AC_junk.display(ax, label="Anticyclonic AVISO not associated", color="navy",transform=data_crs, **kwargs_display)
GIOPS_AC_junk.display(ax, label="Anticyclonic GIOPS not associated", color="darkorange",transform=data_crs, **kwargs_display)
ax.legend(loc="upper left")
ax.set(title='Association based on EddiesObservations\' Methods')
plt.show() |
Nice solution. |
Maybe you could use extern contour to increase match, you could also reduce accepted score. |
Indeed, some associations from the default tracking seems wrong in the above figure. (But the association also seems wrong if I reduce the ratio of overlap.) As for the match function, it would indeed increase the associations, and may be a good solution if we look at large radius eddies. Small eddies won't be properly associated as they may be displaced in GIOPS (their contour may not intersect). So I think I will work a bit more on the default tracking. Here is what I understand from it : I was thinking to put a threshold on the cost values to reduce the faulty association. |
Another things. I talked with a colleague today, and she suggested to also
search match between cyclonic from on set to anticyclonic of the other set.
If there are to much match, maybe dataset will be hard to compare!
Le mer. 27 mai 2020 à 17:19, annesophiefortin <notifications@github.com> a
écrit :
… When the cost < 2
[image: image]
<https://user-images.githubusercontent.com/64930132/83039147-c8abad00-a00b-11ea-83f4-515bb989ee5f.png>
—
You are receiving this because you were assigned.
Reply to this email directly, view it on GitHub
<#7 (comment)>,
or unsubscribe
<https://github.com/notifications/unsubscribe-auth/AIS7AZPYXVY5O63FM6K4UYTRTUVQBANCNFSM4NG5MSYQ>
.
|
Could you remember which field did you for aviso (h, u, v) and for giops, I just look quickly your figure, I am not sure it s a good things that 1/4 of observations are associate with opposite sign did you found same quantity when you do cyclonic aviso vs anticyclonic giops ? |
I am not sure I understand your comment. Are you asking about the sla/adt, ugos, vgos? As for the figure below, the number of eddies are about the same as above. The eddies that are match are about the same size and about at the same location such that the cost is low. |
Did you use same field in the two case, sla or adt(preferred) if it's not the case it could increase mismatch between cyclonic and anticyclonic |
I have no advice to set threshold value on this cost function, because i used rarely this one.
|
It is both sla (will change to adt for my analysis). I use the ugos, vgos of the AVISO file and the add_uv for GIOPS. As GIOPS is interpolated from its native grid to AVISO grid, I think some eddies may be deformed.
|
Hi :) AntSimi showed me your issue. So, my vision is oriented. Thus no intersection = absence in a dataset, an intersection with same rotation with small overlapping area = shift, same rotation and large overlapping area = good agreement, intersection with a different rotation = something is wrong I'm not convinced by the matching eddies which are not overlapping as with the default tracking, for me without overlap it's a different structure, geographical proximity is not enough. But it depends on what is important for your evaluation : the distance, the overlap, the common area, the similarity in characteristics even the eddies are (more or less) shifted... Did you filter the eddies in amplitude, radius or lifetime? It seems that you have few small eddies. |
Hi Dr. Pegliasco, I would like to thank you for your reply. |
Hi Dr. Delepoulle, |
You are right, you could get twice match or more with this method, but it's wanted. Because this case could be highlight interaction process. If you want unique match, you need to apply post processings. |
I would like to compare eddies from observations and models by using the EddyTracking since it gives useful informations like the cost_values for the paired eddies.
However I noticed that, in the tracking.yaml file, TRACK_DURATION_MIN must be at least 3 in order to have Anticyclonic.nc (Cyclonic.nc) and the Anticyclonic_track_too_short.nc (Cyclonic_track_too_short.nc) files. If the it is less than 3, I would obtain only the correspondances and the untracked files. But I would like to compare only 2 files at a times.
Do you think it could be feasible to compare only two files and get the cost_values of the paired eddies or should I develop my own tracking for that?
I don't know if it may be related, but it seems you had a similar issue: #3
Also, as mentionned in that issue, in the wiki (https://github.com/AntSimi/py-eddy-tracker/wiki#how-do-i-create-my-own-tracking-recipe),
from py_eddy_tracker.observations import EddiesObservations as Model
should be from py_eddy_tracker.observations.observation import EddiesObservations as Model
The text was updated successfully, but these errors were encountered: