mlosch/interp_llc/extract_s2004.m

function [H,x,y]=extract_s2004(lon_lo,lon_hi,lat_lo,lat_hi)
%function [H,x,y]=extract_s2004(lon_lo,lon_hi,lat_lo,lat_hi)
%
% extract topography data from GEBCO 1min data set 
%
  
  ncquiet;
  
  load DATAPATH datapath

% $$$   fname = fullfile(datapath,'gebcosubmap.mat');
% $$$   load(fname)
  fname = fullfile(datapath,'s2004_2min.grd');
  ncload(fname)
  dmap = z;
  clear z
% $$$   nc = netcdf(fname,'r');
% $$$   if isempty(nc);
% $$$     error([fname ' not found'])
% $$$   end
  lon = x;
  lat = y;
  
  if lon(1)==lon(end)-360
    % remove last longitude
    lon(end) = [];
    dmap(:,end) = [];
  end
  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 = dmap(jj,ii)';
  
  return
1	mlosch	1.1	function [H,x,y]=extract_s2004(lon_lo,lon_hi,lat_lo,lat_hi)
2			%function [H,x,y]=extract_s2004(lon_lo,lon_hi,lat_lo,lat_hi)
3			%
4			% extract topography data from GEBCO 1min data set
5			%
6
7			ncquiet;
8
9			load DATAPATH datapath
10
11			% $$$ fname = fullfile(datapath,'gebcosubmap.mat');
12			% $$$ load(fname)
13			fname = fullfile(datapath,'s2004_2min.grd');
14			ncload(fname)
15			dmap = z;
16			clear z
17			% $$$ nc = netcdf(fname,'r');
18			% $$$ if isempty(nc);
19			% $$$ error([fname ' not found'])
20			% $$$ end
21			lon = x;
22			lat = y;
23
24			if lon(1)==lon(end)-360
25			% remove last longitude
26			lon(end) = [];
27			dmap(:,end) = [];
28			end
29			nx = length(lon);
30			ny = length(lat);
31			if size(lon,1) > 1; lon = lon'; end
32			if size(lat,1) > 1; lat = lat'; end
33
34			if lon_lo<0
35			i0=find( lon(1:end-1)<=lon_lo+360 & lon(2:end)>lon_lo+360);
36			else
37			i0=find( lon(1:end-1)<=lon_lo & lon(2:end)>lon_lo);
38			end
39			if lon_hi<0
40			i1=find( lon<lon_hi+360 & [lon(2:end) 360]>=lon_hi+360)+1;
41			else
42			i1=find( lon<lon_hi & [lon(2:end) 360]>=lon_hi)+1;
43			end
44			j0=find( [-90 lat(1:end-1)]<=lat_lo & lat>lat_lo)-1; j0=max(1,j0);
45			j1=find( lat<lat_hi & [lat(2:end) 90]>=lat_hi)+1; j1=min(ny,j1);
46			if i0<i1
47			ii=mod((i0:i1)-1,nx)+1;
48			else
49			ii=mod([i0:nx 1:i1]-1,nx)+1;
50			end
51			jj=j0:j1;
52			x = lon(ii);
53			y = lat(jj);
54			H = dmap(jj,ii)';
55
56			return