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

Last change on this file since 27548 was 26683, checked in by Mathieu Morlighem, 3 years ago

CHG: more updates now that we are moving to totten
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 size(T,2)>2 | size(T,1)<1 | size(T,2)<1,
71 error('Size of input T not supported!');
72 end
73 if size(T,2)==1 & any(T(:,1)==1973),
74 disp(' ');
75 disp(' Found year=1973 in T (array). Please, specify the data series using a second index.');
76 disp(' Data available for 1973:');
77 disp(' 1973 1975');
78 disp(' 1973 1984');
79 disp(' 1973 1988');
80 disp(' ');
81 disp(' Usage:');
82 disp(' [vxout vyout]= interpMouginotAntTimeSeries1973to2018(md.mesh.x,md.mesh.y,[1973 1975; 1973 1988; 1991 1992; 2011 2012])');
83 disp(' ');
84 error(' Change input T before continuing.');
85 end %}}}
86 pos = [];
87 for i=1:size(T,1),
88 flag = (T(i,1)==year1);
89 if size(T,2)==2, % ok, check both indexes (year1 and year2)
90 flag = (T(i,1)==year1).*(T(i,2)==year2);
91 end
92 pos = [pos; find(flag)];
93 end
94 % check again {{{
95 if length(pos)~=size(T,1) | length(unique(pos))~=length(pos),
96 disp(' ');
97 disp(' Time resquested does not exist in data set or is repeated!');
98 disp(' Data resquested:');
99 for i=1:length(T(:,1)),
100 str = [' ' int2str(T(i,1)) ' '];
101 if size(T,2)==2, % ok, check both indexes (year1 and year2)
102 str = [str int2str(T(i,2))];
103 end
104 disp(str);
105 end
106 disp(' ');
107 disp(' Data available (24 series):');
108 for i=1:length(year1),
109 str = [' ' int2str(year1(i)) ' ' int2str(year2(i))];
110 disp(str);
111 end
112 disp(' ');
113 disp(' Usage:');
114 disp(' [vxout vyout]= interpMouginotAntTimeSeries1973to2018(md.mesh.x,md.mesh.y,[1986; 1991; 1995; 2000])');
115 disp(' [vxout vyout]= interpMouginotAntTimeSeries1973to2018(md.mesh.x,md.mesh.y,[1986 1988; 1991 1992; 1995 1996; 2000 2001])');
116 disp(' [vxout vyout]= interpMouginotAntTimeSeries1973to2018(md.mesh.x,md.mesh.y,[1973 1975; 1973 1988; 1991 1992; 2011 2012])');
117 disp(' ');
118 error(' Change input T before continuing.');
119 end%}}}
120elseif nargin<3,
121 pos = 1:24; % all available data
122else
123 error('nargin not supported yet!');
124end
125if nargout~=1 & nargout~=2 & nargout~=6
126 error('nargout not supported!');
127end
128
129
130% get the spatial positions
131offset=2;
132
133xmin=min(X(:)); xmax=max(X(:));
134posx=find(xdata<=xmax);
135id1x=max(1,find(xdata>=xmin,1)-offset);
136id2x=min(numel(xdata),posx(end)+offset);
137
138ymin=min(Y(:)); ymax=max(Y(:));
139posy=find(ydata>=ymin);
140id1y=max(1,find(ydata<=ymax,1)-offset);
141id2y=min(numel(ydata),posy(end)+offset);
142
143disp([' -- Mouginot Time Series 1973 to 2018: loading velocities']);
144vxdata = [];
145vydata = [];
146if nargout==6 % it includes ERRX, ERRY, STDX and STDY
147 errxdata = [];
148 errydata = [];
149 stdxdata = [];
150 stdydata = [];
151end
152for i=1:length(pos),
153 disp([' step = ' int2str(i) '/' int2str(length(pos)) ', position = ' int2str(pos(i)) ', year = ' int2str(year1(pos(i))) ' - ' int2str(year2(pos(i)))]);
154 vx = double(ncread(nc,'VX',[id1x id1y pos(i)],[id2x-id1x+1 id2y-id1y+1 1],[1 1 1]));
155 vy = double(ncread(nc,'VY',[id1x id1y pos(i)],[id2x-id1x+1 id2y-id1y+1 1],[1 1 1]));
156 vxdata(:,:,i) = permute(vx,[2 1 3]);
157 vydata(:,:,i) = permute(vy,[2 1 3]);
158 if nargout==6 % it includes ERRX, ERRY, STDX and STDY
159 errx = double(ncread(nc,'ERRX',[id1x id1y pos(i)],[id2x-id1x+1 id2y-id1y+1 1],[1 1 1]));
160 erry = double(ncread(nc,'ERRY',[id1x id1y pos(i)],[id2x-id1x+1 id2y-id1y+1 1],[1 1 1]));
161 stdx = double(ncread(nc,'STDX',[id1x id1y pos(i)],[id2x-id1x+1 id2y-id1y+1 1],[1 1 1]));
162 stdy = double(ncread(nc,'STDY',[id1x id1y pos(i)],[id2x-id1x+1 id2y-id1y+1 1],[1 1 1]));
163 errxdata(:,:,i) = permute(errx,[2 1 3]);
164 errydata(:,:,i) = permute(erry,[2 1 3]);
165 stdxdata(:,:,i) = permute(stdx,[2 1 3]);
166 stdydata(:,:,i) = permute(stdy,[2 1 3]);
167 end
168end
169xdata=xdata(id1x:id2x);
170ydata=ydata(id1y:id2y);
171
172disp([' -- Mouginot Time Series 1973 to 2018: interpolating']);
173vxout = [];
174vyout = [];
175if nargout==6 % it includes ERRX, ERRY, STDX and STDY
176 errxout = [];
177 erryout = [];
178 stdxout = [];
179 stdyout = [];
180end
181for i=1:length(pos),
182 disp([' step = ' int2str(i) '/' int2str(length(pos)) ', position = ' int2str(pos(i)) ', year = ' int2str(year1(pos(i))) ' - ' int2str(year2(pos(i)))]);
183 vxout = [vxout InterpFromGrid(xdata,ydata,vxdata(:,:,i),double(X),double(Y))];
184 vyout = [vyout InterpFromGrid(xdata,ydata,vydata(:,:,i),double(X),double(Y))];
185 if nargout==6 % it includes ERRX, ERRY, STDX and STDY
186 errxout = [errxout InterpFromGrid(xdata,ydata,errxdata(:,:,i),double(X),double(Y))];
187 erryout = [erryout InterpFromGrid(xdata,ydata,errydata(:,:,i),double(X),double(Y))];
188 stdxout = [stdxout InterpFromGrid(xdata,ydata,stdxdata(:,:,i),double(X),double(Y))];
189 stdyout = [stdyout InterpFromGrid(xdata,ydata,stdydata(:,:,i),double(X),double(Y))];
190 end
191end
192
193%return vel if only one output is requested
194if nargout==1,
195 vxout = sqrt(vxout.^2+vyout.^2);
196end
Note: See TracBrowser for help on using the repository browser.