source: issm/trunk-jpl/src/m/modeldata/interpMouginotAntTimeSeries1973to2018.m

Last change on this file was 27765, checked in by Mathieu Morlighem, 22 months ago

CHG: allow for any vector as time stamps
File size: 6.8 KB
Line 
1function [vxout vyout errxout erryout stdxout stdyout]= interpMouginotAntTimeSeries1973to2018(X,Y,T)
2%INTERPMOUGINOTANTTIMESERIES1973TO2018 - interpolate observed (time series) velocities
3%
4% Inputs
5% X,Y: spatial (scatter) coordinates
6% T: time (indexed by YEAR1 (YEAR2 is optional); see below)
7%
8% Outputs
9% vxout,vyout: interpolated velocities at X,Y, for each time requested in T
10%
11% Available time series:
12%
13% YEAR1 YEAR2
14% 1 1973 1975
15% 2 1973 1984
16% 3 1973 1988
17% 4 1984 1988
18% 5 1986 1988
19% 6 1988 1990
20% 7 1991 1992
21% 8 1995 1996
22% 9 2000 2001
23% 10 2002 2003
24% 11 2003 2004
25% 12 2005 2006
26% 13 2006 2007
27% 14 2007 2008
28% 15 2008 2009
29% 16 2009 2010
30% 17 2010 2011
31% 18 2011 2012
32% 19 2012 2013
33% 20 2013 2014
34% 21 2014 2015
35% 22 2015 2016
36% 23 2016 2017
37% 24 2017 2018
38%
39% Usage:
40% T refers to YEAR1, but the user can also use YEAR2 (e.g., the "1973" case in YEAR1).
41%
42% Then, these codes generate the same results:
43%
44% [vxout vyout]= interpMouginotAntTimeSeries1973to2018(md.mesh.x,md.mesh.y,[1986; 1991; 1995; 2000]);
45% [vxout vyout]= interpMouginotAntTimeSeries1973to2018(md.mesh.x,md.mesh.y,[1986 1988; 1991 1992; 1995 1996; 2000 2001]);
46%
47% Another examples:
48% [vxout vyout]= interpMouginotAntTimeSeries1973to2018(md.mesh.x,md.mesh.y,[1973 1975; 1973 1988; 1991 1992; 2011 2012]);
49% [vel]= interpMouginotAntTimeSeries1973to2018(md.mesh.x,md.mesh.y,[1986; 1991; 1995; 2000]);
50% [vxout vyout errxout erryout stdxout stdyout]= interpMouginotAntTimeSeries1973to2018(md.mesh.x,md.mesh.y,[1986; 1991; 1995; 2000]);
51
52%read data
53switch (oshostname()),
54 case {'ronne'}
55 nc = '/home/ModelData/Antarctica/MouginotVel/ASE_TimeSeries_1973-2018.nc';
56 case {'totten'}
57 nc = '/totten_1/ModelData/Antarctica/MouginotVel/ASE_TimeSeries_1973-2018.nc';
58 otherwise
59 error('hostname not supported yet');
60end
61
62xdata = double(ncread(nc,'x'));
63ydata = double(ncread(nc,'y'));
64year1 = ncread(nc,'YEAR1');
65year2 = ncread(nc,'YEAR2');
66
67% get the positions related to T
68if nargin==3
69 % initial checks %{{{
70 if ~isvector(T)
71 error('Size of input T not supported!');
72 end
73 T=T(:);
74 if size(T,2)==1 & any(T(:,1)==1973),
75 disp(' ');
76 disp(' Found year=1973 in T (array). Please, specify the data series using a second index.');
77 disp(' Data available for 1973:');
78 disp(' 1973 1975');
79 disp(' 1973 1984');
80 disp(' 1973 1988');
81 disp(' ');
82 disp(' Usage:');
83 disp(' [vxout vyout]= interpMouginotAntTimeSeries1973to2018(md.mesh.x,md.mesh.y,[1973 1975; 1973 1988; 1991 1992; 2011 2012])');
84 disp(' ');
85 error(' Change input T before continuing.');
86 end %}}}
87 pos = [];
88 for i=1:size(T,1),
89 flag = (T(i,1)==year1);
90 if size(T,2)==2, % ok, check both indexes (year1 and year2)
91 flag = (T(i,1)==year1).*(T(i,2)==year2);
92 end
93 pos = [pos; find(flag)];
94 end
95 % check again {{{
96 if length(pos)~=size(T,1) | length(unique(pos))~=length(pos),
97 disp(' ');
98 disp(' Time resquested does not exist in data set or is repeated!');
99 disp(' Data resquested:');
100 for i=1:length(T(:,1)),
101 str = [' ' int2str(T(i,1)) ' '];
102 if size(T,2)==2, % ok, check both indexes (year1 and year2)
103 str = [str int2str(T(i,2))];
104 end
105 disp(str);
106 end
107 disp(' ');
108 disp(' Data available (24 series):');
109 for i=1:length(year1),
110 str = [' ' int2str(year1(i)) ' ' int2str(year2(i))];
111 disp(str);
112 end
113 disp(' ');
114 disp(' Usage:');
115 disp(' [vxout vyout]= interpMouginotAntTimeSeries1973to2018(md.mesh.x,md.mesh.y,[1986; 1991; 1995; 2000])');
116 disp(' [vxout vyout]= interpMouginotAntTimeSeries1973to2018(md.mesh.x,md.mesh.y,[1986 1988; 1991 1992; 1995 1996; 2000 2001])');
117 disp(' [vxout vyout]= interpMouginotAntTimeSeries1973to2018(md.mesh.x,md.mesh.y,[1973 1975; 1973 1988; 1991 1992; 2011 2012])');
118 disp(' ');
119 error(' Change input T before continuing.');
120 end%}}}
121elseif nargin<3,
122 pos = 1:24; % all available data
123else
124 error('nargin not supported yet!');
125end
126if nargout~=1 & nargout~=2 & nargout~=6
127 error('nargout not supported!');
128end
129
130
131% get the spatial positions
132offset=2;
133
134xmin=min(X(:)); xmax=max(X(:));
135posx=find(xdata<=xmax);
136id1x=max(1,find(xdata>=xmin,1)-offset);
137id2x=min(numel(xdata),posx(end)+offset);
138
139ymin=min(Y(:)); ymax=max(Y(:));
140posy=find(ydata>=ymin);
141id1y=max(1,find(ydata<=ymax,1)-offset);
142id2y=min(numel(ydata),posy(end)+offset);
143
144disp([' -- Mouginot Time Series 1973 to 2018: loading velocities']);
145vxdata = [];
146vydata = [];
147if nargout==6 % it includes ERRX, ERRY, STDX and STDY
148 errxdata = [];
149 errydata = [];
150 stdxdata = [];
151 stdydata = [];
152end
153for i=1:length(pos),
154 disp([' step = ' int2str(i) '/' int2str(length(pos)) ', position = ' int2str(pos(i)) ', year = ' int2str(year1(pos(i))) ' - ' int2str(year2(pos(i)))]);
155 vx = double(ncread(nc,'VX',[id1x id1y pos(i)],[id2x-id1x+1 id2y-id1y+1 1],[1 1 1]));
156 vy = double(ncread(nc,'VY',[id1x id1y pos(i)],[id2x-id1x+1 id2y-id1y+1 1],[1 1 1]));
157 vxdata(:,:,i) = permute(vx,[2 1 3]);
158 vydata(:,:,i) = permute(vy,[2 1 3]);
159 if nargout==6 % it includes ERRX, ERRY, STDX and STDY
160 errx = double(ncread(nc,'ERRX',[id1x id1y pos(i)],[id2x-id1x+1 id2y-id1y+1 1],[1 1 1]));
161 erry = double(ncread(nc,'ERRY',[id1x id1y pos(i)],[id2x-id1x+1 id2y-id1y+1 1],[1 1 1]));
162 stdx = double(ncread(nc,'STDX',[id1x id1y pos(i)],[id2x-id1x+1 id2y-id1y+1 1],[1 1 1]));
163 stdy = double(ncread(nc,'STDY',[id1x id1y pos(i)],[id2x-id1x+1 id2y-id1y+1 1],[1 1 1]));
164 errxdata(:,:,i) = permute(errx,[2 1 3]);
165 errydata(:,:,i) = permute(erry,[2 1 3]);
166 stdxdata(:,:,i) = permute(stdx,[2 1 3]);
167 stdydata(:,:,i) = permute(stdy,[2 1 3]);
168 end
169end
170xdata=xdata(id1x:id2x);
171ydata=ydata(id1y:id2y);
172
173disp([' -- Mouginot Time Series 1973 to 2018: interpolating']);
174vxout = [];
175vyout = [];
176if nargout==6 % it includes ERRX, ERRY, STDX and STDY
177 errxout = [];
178 erryout = [];
179 stdxout = [];
180 stdyout = [];
181end
182for i=1:length(pos),
183 disp([' step = ' int2str(i) '/' int2str(length(pos)) ', position = ' int2str(pos(i)) ', year = ' int2str(year1(pos(i))) ' - ' int2str(year2(pos(i)))]);
184 vxout = [vxout InterpFromGrid(xdata,ydata,vxdata(:,:,i),double(X),double(Y))];
185 vyout = [vyout InterpFromGrid(xdata,ydata,vydata(:,:,i),double(X),double(Y))];
186 if nargout==6 % it includes ERRX, ERRY, STDX and STDY
187 errxout = [errxout InterpFromGrid(xdata,ydata,errxdata(:,:,i),double(X),double(Y))];
188 erryout = [erryout InterpFromGrid(xdata,ydata,errydata(:,:,i),double(X),double(Y))];
189 stdxout = [stdxout InterpFromGrid(xdata,ydata,stdxdata(:,:,i),double(X),double(Y))];
190 stdyout = [stdyout InterpFromGrid(xdata,ydata,stdydata(:,:,i),double(X),double(Y))];
191 end
192end
193
194%return vel if only one output is requested
195if nargout==1,
196 vxout = sqrt(vxout.^2+vyout.^2);
197end
Note: See TracBrowser for help on using the repository browser.