AntSimi / py-eddy-tracker

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

Collocated External Data. #24

Closed acapet closed 4 years ago

acapet commented 4 years ago

Hi !

I'm trying to load external data (SST) using the same data structure for collocation (eg. assess mean SST anomalies within contours).

I think I'm almost there, just miss the interpolation step.

datest = '20160721'
SLAfilename = "../l4_"+set1+"/dt_blacksea_allsat_phy_l4_"+datest+"_"+dat1+".nc"
SSTfilename = "../SST/"+datest+"000000-GOS-L4_GHRSST-SSTfnd-OISST_HR_REP-BLK-v02.0-fv01.0.nc4" 

g = RegularGridDataset(SLAfilename, x_name= 'longitude',y_name= 'latitude') 
t = RegularGridDataset(filename=SSTfilename, x_name="lon", y_name="lat") 
t.grid('analysed_sst')

So far, so good. Then I'd like to use :

g.add_grid('sst',ti)

with

ti=t.interp('analysed_sst',g.grid(lon_name),g.grid(lat_name))

but something is wrong.

Among what I noticed :

Thanks for insights !!

Art

AntSimi commented 4 years ago

Hi!

I think you should try something like this :

from numpy import meshgrid
# create mesh
lats, lons = meshgrid(g.y_c, g.x_c)
shape = lats.T.shape
# flat grid before interp
lons, lats = lons.reshape(-1), lats.reshape(-1)
# interp and reshape
ti = g.interp('sst', lons, lats).reshape(shape)

In my mind masked array or not must be the same, but when i read function ... if it's a problem send traceback

I wonder why you want interpolate sst grid on sla resolution before assess mean SST anomalies within contours? You could get sst pixel within contours directly .

Antoine

acapet commented 4 years ago

Thanks, that's exactly what I needed. For the record, here's how I got through :

from numpy import meshgrid
lons, lats = meshgrid(g.x_c, g.y_c)
shape = lats.shape

# flat grid before interp
lons, lats = lons.reshape(-1), lats.reshape(-1)

# interp and reshape
ti = t.interp('analysed_sst', lons, lats).reshape(shape).T
ti = ma.masked_invalid(ti)

g.add_grid('sst',ti)

I masked it for later use e.g. g.add_grid('sst_anom',ti-ti.mean()).

The reason I loaded SST everywhere, and not only within contours at this stage, is for a first visual inspection.

Thanks !!!

(PS: I'm preparing a show case with corresponding data files for the gallery)

acapet commented 4 years ago

It seems that add_grid doesn't complete g.variables_description. As a consequence, I can't for instance use

g.copy("sst", "sst_high")
g.bessel_high_filter('sst_high',100)

to higlight anomalies in SST.

Any hint ?

acapet commented 4 years ago

Ok, I went over this issue with

g.variables_description['sst'] = t.variables_description['analysed_sst']

but I'm unsure it would be entirely correct, since there has been an interpolation in between these two..

AntSimi commented 4 years ago

Good, maybe i will add a method to re-grid a grid on another one (with same code of interpolation) to take all informations in one command