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	mlosch	1.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