AntSimi / py-eddy-tracker

Eddy identification and tracking
https://py-eddy-tracker.readthedocs.io/en/latest/
GNU General Public License v3.0
123 stars 51 forks source link

Tracks on land and longitude interval #76

Closed davcali closed 3 years ago

davcali commented 3 years ago

Hello,

I'm using your eddy tracker in the Mediterranean area over a period of 27 years, and in the North Atlantic area over a period of 10 year. I think I've run the detection and tracking steps correctly (or I hope so), but I have two questions:

1) why do I find tracks passing over the land? How can I avoid it?

Tracks_carto

tracks_longer20weeks

2) in the inputl netcdf file, longitudes are in the interval [-180°, 180°] , but the detection process generates Anticyclonic and Cyclonic netcdf files with longitudes in [0°, 360°]. How can I display the contours and the tracks in the original [-180°,180°] longitudes system? I'm pretty new to Python and maybe I'm missing something easy... For clarity I would like to report the eddies in this pic:

Eddies_not_centered_in_zero

in the area here:

Correct_area

For the first figure (not centered), I'm using the following code:

from matplotlib import pyplot as plt

from py_eddy_tracker import data
from py_eddy_tracker.observations.observation import EddiesObservations
from py_eddy_tracker.dataset.grid import RegularGridDataset

Load detection files

a = EddiesObservations.load_file("out_directory_NA/Anticyclonic_20100101.nc")
c = EddiesObservations.load_file("out_directory_NA/Cyclonic_20100101.nc")

Plot the speed and effective (dashed) contours

fig = plt.figure(figsize=(15, 8))
ax = fig.add_axes((0.05, 0.05, 0.9, 0.9))
ax.set_aspect("equal")
ax.set_xlim(0,360)
ax.set_ylim(20, 90)

g = RegularGridDataset(
    "CMEMS_DATA_NA_2010_2020_dailyzied/20100101.nc", "lon", "lat"
)
m = g.display(ax, "adt", vmin=-0.15, vmax=0.15, cmap="RdBu_r")
a.display(ax, label="Anticyclonic contour", color= "k", lw=1.5)
c.display(ax, label="Cyclonic contour",color= "k", ls = "dashed", lw=1.5)

ax.legend(loc="upper right")

plt.show()`

Thank you very much!

AntSimi commented 3 years ago

For longitude wrapping use 'ref' option documentation is here

AntSimi commented 3 years ago

I don't know exactly how your run tracking, but I guess that you use default tracker describe here which are clearly not adapted for med area and which didn't check land. You could use a custom tracker like areatracker which are more stable maybe with use of some virtual.

davcali commented 3 years ago

Hello Antoine,

thanks for the quick reply!

Did you use the area-tracker in the example on the website?

https://py-eddy-tracker.readthedocs.io/en/latest/python_module/08_tracking_manipulation/pet_display_track.html#sphx-glr-python-module-08-tracking-manipulation-pet-display-track-py

Does the defaul tracker not check for land in general? Do you think that for the North Atlantic area I could use the default tracker? In any case, I will surely try with the area-tracker.

Thanks for the link to the documentation for the "ref" option, I didn't really see it!

Thank you very much,

Davide

AntSimi commented 3 years ago

Sorry i don't remenber how i build those dataset. But you must consider that py-eddy-tracker is a toolbox, so default behaviour isn't appropriate in any case. My first advice is to use AreaTracker with virtual set from 2 to 5, but you could also customize your tracker : https://github.com/AntSimi/py-eddy-tracker/wiki

AntSimi commented 3 years ago

Here you could see a way how to manage track across land, but i think isn't usefull for AreaTracker.

davcali commented 3 years ago

Hello Antoine, thanks for your advices. I've tried with both AreaTracker and CheltonTracker. The problem with the land is persisting, and this seems a little bit weird to me. The Eddy Atlas provided on the AVISO website is produced with the Chelton Tracker, and it correctly gives the tracks avoiding the land. The link you gave in the last reply is related to the Chelton Tracker

https://github.com/AntSimi/py-eddy-tracker/blob/master/src/py_eddy_tracker/featured_tracking/old_tracker_reference.py#L46

Should it check for land, no? Is there something that I'm missing? In your opinion, in case, how and where could I introduce a check for the land in the Chelton or in the Area tracker? This is the only real issue I'm facing with your tracker (for the rest, it seems to work pretty fine)

Thanks!

AntSimi commented 3 years ago

So AVISO atlas is produce with py-eddy-tracker but with custom tracking function which are not disclosed. AVISO atlas doesn't check if there are land on eddy path but check the distance of the closest coast along eddy path and put an arbitrary threshold on this value which must greater than 20% of radius.

I am surprised, that you have land problem with AreaTracker(CheltonTracker also which have specific function across_ground), how did you setup tracking? Could you show or share some little example with AreaTracker which highlight problem? And share how did you run tracking?

In my opinion with AreaTracker there are no need of landmask, or maybe in specific region like maldive or aegan sea where altimetry product are defined over eddy island.

davcali commented 3 years ago

So AVISO atlas is produce with py-eddy-tracker but with custom tracking function which are not disclosed. AVISO atlas doesn't check if there are land on eddy path but check the distance of the closest coast along eddy path and put an arbitrary threshold on this value which must greater than 20% of radius.

Ah, ok, now it makes sense.

About the tracking, after the identification, I run a conf.yaml file like this:

PATHS:

  FILES_PATTERN: /home/davide/data/NA_CMEMS_2010_2020/out_directory_NA/Cyclonic_*.nc
  SAVE_DIR: /home/davide/data/NA_CMEMS_2010_2020

VIRTUAL_LENGTH_MAX: 3       
TRACK_DURATION_MIN: 5

CLASS:
     MODULE: py_eddy_tracker.featured_tracking.area_tracker
     CLASS: AreaTracker        #or CheltonTracker

What I obtain in this way (using AreaTracker or CheltonTracker):

Area_tracker Chelton_tracker_2

As you can see there are tracks over the land...

What do you think?

D

AntSimi commented 3 years ago

I guess you used filtered path? If it's case could you remove filtering step for my understanding. Could you also add some longitude/latitude label in order to estimate scale? Could you give informations about input grid specification?

davcali commented 3 years ago

Yes, correct, I used the filtering. Now it is removed, and I added some labels:

More_detail_1

About the input, these are CMEMS 1/4 degree observational data (is this the information you were asking about?)

AntSimi commented 3 years ago

Ok if it's AreaTracker => i am not surprise to see some track over island, because i am pretty sure those pixel are defined in altimetry maps. For CheltonTracker and across_ground method mask used is here, PyEddyTracker check if the direct path between 2 observations doesn't cross a land pixel . But the mask is define at 1/60, so it's pretty accurate.

It could be possible to overload AreaTracker to use this across_ground method.

davcali commented 3 years ago

Yes, the last image is made with Area Tracker, so I would agree with you. But it seems that also the Chelton Tracker has the same problem, as you can see in the pic I posted before (which is made with the Chelton Tracker)

Chelton_tracker_2

I think your line of reasoning is correct, but I don't undestand why, at this point, I get this problem also with the Chelton Tracker...

AntSimi commented 3 years ago

When you used CheltonTracker you change MODULE value i guess?

Maybe display mask for this area, it was computed with full resolution gshhs data. so i think it will be coherent with your coastline. Again could you remove filtered step for display.

davcali commented 3 years ago

Yes, I changed it, I used

MODULE: py_eddy_tracker.featured_tracking.old_tracker_reference

The mask seems to be coherent with my coastlines:

mask_1_60

No, I can't remove the filtering (my VM is busy with another process), but I think it doesn't make a great difference...

AntSimi commented 3 years ago

So could you share just few days of detection maybe (5-10 days) and yaml from tracking to understand this case. Less data is the best, if you could show the problem on small area and only on few days ...

davcali commented 3 years ago

Could be this the problem? (days 1, 2 and 12)

Step1 Step2 Step12

AntSimi commented 3 years ago

Here you have detection over land but eddy path which is built with center will not cross land, that why aviso atlas use a grid which contains distance with closest coast to be able to see coast on side of eddy path.

AntSimi commented 3 years ago

In my opinion this case is OK because eddy may exists but maybe a little bit smaller. In order to have a better understanding you could subset your atlas on a shorter period and area with Eddysubsetter and animate with EddyAnim