| 1 |
function [madt]=aviso_get_madt(loc,grid_aviso) |
| 2 |
% function [madt]=aviso_get_madt(loc,grid) |
| 3 |
% extract madt value (mapped absolute dynamic topography) at 1 location |
| 4 |
% linear interpolation in space and time at each location. |
| 5 |
% |
| 6 |
% loc is a structure with three array of size (N,1) : |
| 7 |
% loc.lat : latitude [deg_north] |
| 8 |
% loc.lon : longitude [deg_east] |
| 9 |
% loc.time : date in Matlab julian date format |
| 10 |
% |
| 11 |
% grid is the structure returned by aviso_get_grid |
| 12 |
% |
| 13 |
% madt is the output array [unit: m], of size (N,1). |
| 14 |
% |
| 15 |
|
| 16 |
N=length(loc.lat); |
| 17 |
madt=NaN*zeros(N,1); |
| 18 |
|
| 19 |
loc.lon(loc.lon<0)=loc.lon(loc.lon<0)+360; |
| 20 |
|
| 21 |
Itime=interp1(grid_aviso.time,(1:grid_aviso.Ntime)'-1,loc.time); |
| 22 |
Ilon=interp1(grid_aviso.NbLongitudes,(1:grid_aviso.Nlon)'-1,loc.lon); |
| 23 |
Ilat=interp1(grid_aviso.NbLatitudes,(1:grid_aviso.Nlat)'-1,loc.lat); |
| 24 |
|
| 25 |
[X,Y,T]=meshgrid([0 1],[0 1],[0 1]); |
| 26 |
|
| 27 |
for kk=1:N, |
| 28 |
|
| 29 |
if isnan(Itime(kk)*Ilon(kk)*Ilat(kk)), continue, end |
| 30 |
it=floor(Itime(kk)); |
| 31 |
il=floor(Ilon(kk)); |
| 32 |
iL=floor(Ilat(kk)); |
| 33 |
URL=sprintf('http://opendap.aviso.oceanobs.com/thredds/dodsC/dataset-duacs-dt-ref-global-merged-madt-h?Grid_0001[%d:%d][%d:%d][%d:%d]',... |
| 34 |
it,it+1,il,il+1,iL,iL+1); |
| 35 |
a=loaddap(URL);data=a.Grid_0001; % data.Grid_0001 : lon x lat x time |
| 36 |
if any(data.Grid_0001(:)>1e10), continue, end |
| 37 |
madt(kk)=interp3(X,Y,T,data.Grid_0001/100,Ilon(kk)-il,Ilat(kk)-iL,Itime(kk)-it); |
| 38 |
|
| 39 |
end |