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	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