mlosch/interp_llc/extract_brios.m

function [H,Z,x,y]=extract_brios(lon_lo,lon_hi,lat_lo,lat_hi)
%function [H,Z,x,y]=extract_brios(lon_lo,lon_hi,lat_lo,lat_hi)
%
% extract topography data from BRIOS data set
%
  
  ncquiet;
  
  load DATAPATH datapath

  fname = fullfile(datapath,'shelfice.mat');
  theResult = load(fname);
  if isempty(theResult);
    error([fname ' not found'])
  end
% $$$   lon = theResult.x+180;
% $$$   lat = theResult.y;
% $$$   nx = length(lon);
% $$$   ny = length(lat);
% $$$   if size(lon,1) > 1; lon = lon'; end
% $$$   if size(lat,1) > 1; lat = lat'; end
% $$$   
% $$$   if lon_lo<0
% $$$     i0=find( lon(1:end-1)<=lon_lo+360 & lon(2:end)>lon_lo+360);
% $$$   else
% $$$     i0=find( lon(1:end-1)<=lon_lo & lon(2:end)>lon_lo);
% $$$   end
% $$$   if lon_hi<0
% $$$     i1=find( lon<lon_hi+360 & [lon(2:end) 360]>=lon_hi+360)+1;
% $$$   else
% $$$     i1=find( lon<lon_hi & [lon(2:end) 360]>=lon_hi)+1;
% $$$   end
% $$$   j0=find( [-90 lat(1:end-1)]<=lat_lo & lat>lat_lo)-1; j0=max(1,j0);
% $$$   j1=find( lat<lat_hi & [lat(2:end) 90]>=lat_hi)+1; j1=min(ny,j1);
% $$$   if i0<i1
% $$$     ii=mod((i0:i1)-1,nx)+1;
% $$$   else
% $$$     ii=mod([i0:nx 1:i1]-1,nx)+1;
% $$$   end
% $$$   jj=j0:j1;
% $$$   
% $$$   x = lon(ii);
% $$$   y = lat(jj);
% $$$   H = -theResult.h(ii,jj);
% $$$   Z = -theResult.z(ii,jj);

  x = theResult.x(:)+180;
  y = theResult.y(:);
  H = -theResult.h;
  Z = -theResult.z;
  
  return
1	function [H,Z,x,y]=extract_brios(lon_lo,lon_hi,lat_lo,lat_hi)
2	%function [H,Z,x,y]=extract_brios(lon_lo,lon_hi,lat_lo,lat_hi)
3	%
4	% extract topography data from BRIOS data set
5	%
6
7	ncquiet;
8
9	load DATAPATH datapath
10
11	fname = fullfile(datapath,'shelfice.mat');
12	theResult = load(fname);
13	if isempty(theResult);
14	error([fname ' not found'])
15	end
16	% $$$ lon = theResult.x+180;
17	% $$$ lat = theResult.y;
18	% $$$ nx = length(lon);
19	% $$$ ny = length(lat);
20	% $$$ if size(lon,1) > 1; lon = lon'; end
21	% $$$ if size(lat,1) > 1; lat = lat'; end
22	% $$$
23	% $$$ if lon_lo<0
24	% $$$ i0=find( lon(1:end-1)<=lon_lo+360 & lon(2:end)>lon_lo+360);
25	% $$$ else
26	% $$$ i0=find( lon(1:end-1)<=lon_lo & lon(2:end)>lon_lo);
27	% $$$ end
28	% $$$ if lon_hi<0
29	% $$$ i1=find( lon<lon_hi+360 & [lon(2:end) 360]>=lon_hi+360)+1;
30	% $$$ else
31	% $$$ i1=find( lon<lon_hi & [lon(2:end) 360]>=lon_hi)+1;
32	% $$$ end
33	% $$$ j0=find( [-90 lat(1:end-1)]<=lat_lo & lat>lat_lo)-1; j0=max(1,j0);
34	% $$$ j1=find( lat<lat_hi & [lat(2:end) 90]>=lat_hi)+1; j1=min(ny,j1);
35	% $$$ if i0<i1
36	% $$$ ii=mod((i0:i1)-1,nx)+1;
37	% $$$ else
38	% $$$ ii=mod([i0:nx 1:i1]-1,nx)+1;
39	% $$$ end
40	% $$$ jj=j0:j1;
41	% $$$
42	% $$$ x = lon(ii);
43	% $$$ y = lat(jj);
44	% $$$ H = -theResult.h(ii,jj);
45	% $$$ Z = -theResult.z(ii,jj);
46
47	x = theResult.x(:)+180;
48	y = theResult.y(:);
49	H = -theResult.h;
50	Z = -theResult.z;
51
52	return