ecco-godae/climatology/gen_bathy.m

%
%  Ed Hill
%
%  Generate approximate bathymetry for the llc grid
%

%  clear all ; close all

do_print = 0;
dbug = 0;

%  Get the ETOPO2 data
disp('Reading the ETOPO2 data...');
nlat =  5400;
nlon = 10800;
lons = linspace(-180, (180 - 2/60), nlon);
lats = linspace(90, -(90 - 2/60), nlat);
fid = fopen('ETOPO2.raw.bin', 'r', 'ieee-be');
et2 =  reshape( fread(fid, nlat*nlon, 'int16'), [ nlon nlat ]);
fid = fclose(fid);

if dbug > 10
  surf( et2(1:500,1:500) ), view(2), shading flat
  ilon = [  1:10:5000 ];
  ilat = [ 1:10:5400 ];
  surf(lons(ilon), lats(ilat), et2(ilon,ilat)'), view(2), ...
      shading flat, axis equal
end

%  Get the llc grid information
disp('Reading the grid information...');
fnc_in = 'llc_p90_%d.nc';
vnms = { 'XC' 'YC' 'XG' 'YG' };
d2r = pi/180.0;
ginfo = {};
for k = 1:5
  nc = netcdf(sprintf(fnc_in,k),'nowrite');
  for j = 1:length(vnms)
    eval(sprintf('ginfo(k).%s = nc{''%s''}(:);',vnms{j},vnms{j}));
  end
  nc = close(nc);
  
  %  compute 3D coords of the C and G points
  cor = zeros( [ size(ginfo(k).XG) 3 ]);
  [cor(:,:,1),cor(:,:,2),cor(:,:,3)] = sph2cart( ginfo(k).XG*d2r, ...
                                                 ginfo(k).YG*d2r,1 );
  cen = zeros( [ size(ginfo(k).XC) 3 ]);
  [cen(:,:,1),cen(:,:,2),cen(:,:,3)] = sph2cart( ginfo(k).XC*d2r, ...
                                                 ginfo(k).YC*d2r,1 );
  ginfo(k).cor = cor;
  ginfo(k).cen = cen;
  
  %  empty bathy
  ginfo(k).bathy = zeros(size( size(ginfo(k).XC) ));
end

%  Do a quick-and-dirty regrid
disp('Performing the regrid...');
%ph(
%%%s_lon = 0.0;
s_lon = -37.0;
%ph)
s_lat = 0.0;
for k = 1:5
  disp(sprintf('  k = %d',k));
  %
  %  /  i,j+1 i+1,j+1 \
  %  \  i,j   i+i,j   /
  for j = 1:size(ginfo(k).XC,2)
    jr = [ j j+1 ];
    for i = 1:size(ginfo(k).XC,1)
      ir = [ i i+1 ];
      lon_min = min(min( ginfo(k).XG(ir,jr) + s_lon ));
      lon_max = max(max( ginfo(k).XG(ir,jr) + s_lon ));
      lat_min = min(min( ginfo(k).YG(ir,jr) + s_lat ));
      lat_max = max(max( ginfo(k).YG(ir,jr) + s_lat ));
      % [ lon_min lon_max lat_min lat_max ]
      lon_mm = [ mod(lon_min + 180,360)-180 ...
                 mod(lon_max + 180,360)-180 ];
      lon_min = min(lon_mm);
      lon_max = max(lon_mm);

      d_lon = lon_max - lon_min;
      if (abs(d_lon) > 45 && lat_min < 88)
        i_lon = find( lons >= lon_max  ...
                      |  lons <= lon_min );
      else
        i_lon = find(lon_min <= lons  &  lons <= lon_max);
      end
      i_lat = find(lat_min <= lats & lats <= lat_max);
      
      %  average the bathy
      bv = et2(i_lon,i_lat);
      if length(bv(:)) > 1
        ginfo(k).bathy(i,j) = mean(bv(:));
      else
        %  If not averaging, then interpolate
        %  ginfo(k).bathy(i,j) = NaN;
        i_lon = interp1(lons,[1:length(lons)],lon_min,'nearest');
        i_lat = interp1(lats,[1:length(lats)],lat_min,'nearest');
        ginfo(k).bathy(i,j) = et2(i_lon,i_lat);
      end
    end
  end
end

%  Plot the resulting bathymetry
k = 3;
for k = 1:5
  disp(sprintf('  k = %d',k));
  ginfo(k).b_oce = ginfo(k).bathy;
  %  ginfo(k).b_oce(find(ginfo(k).b_oce >= 0.0)) = NaN;
  ginfo(k).b_oce(find(ginfo(k).b_oce >= 0.0)) = 1000;
  if k == 2
    hold on
  end
  %  figure(k)
  surf( ginfo(k).cor(:,:,1), ...
        ginfo(k).cor(:,:,2), ...
        ginfo(k).cor(:,:,3), ginfo(k).b_oce )
  if k == 5
    hold off
  end
end
axis equal, view(2)


%  Plot the bathymetry at 1/2 resolution
k = 3;
for k = 1:5
  disp(sprintf('  k = %d',k));
  ginfo(k).c_half = ginfo(k).cor(1:2:end,1:2:end,:);
  ginfo(k).b_half = ginfo(k).bathy(1:2:end,1:2:end);
  ginfo(k).b_half(find(ginfo(k).b_half >= 0.0)) = 1000;
  if k == 2
    hold on
  end
  %  figure(k)
  surf( ginfo(k).c_half(:,:,1), ...
        ginfo(k).c_half(:,:,2), ...
        ginfo(k).c_half(:,:,3), ginfo(k).b_half )
  if k == 5
    hold off
  end
end
axis equal, view(2)

%--- BREAK POINT HERE
%--- do e.g. view(80,45) to change view angle


if do_print > 0
  print -depsc llc_bathy.eps
end


%  Write the bathymetry to netCDF compatible with MNC
id = 'bathy';
units = 'm';
for k = 1:5
  disp(sprintf('  k = %d',k));
  ginfo(k).b_out = ginfo(k).bathy;
  ginfo(k).b_out(find(ginfo(k).bathy >= 0.0)) = 0.0;
  bfile = sprintf('llc_p90_bathy.f%03d.nc',k);
  
  nc = netcdf(bfile, 'clobber');
  nc.reference = 'bathymetry';
  nc.author = 'Ed Hill <eh3@mit.edu>';
  nc.date = eval('date');
  nc('X') = size( ginfo(k).XC, 1 );
  nc('Y') = size( ginfo(k).XC, 2 );
  nc{ id } = { 'Y', 'X' };
  nc{ id }.units = units;
  nc{ id }(:) = ginfo(k).b_out';
  nc = close(nc);
end

%  Write the bathymetry to a set of 90-by-90 per-tile MDSIO-compatible
%  files
%       face  is  js
iti = [   1    1   1  ; ...
          1   91   1  ; ...
          1  181   1  ; ...
          1  271   1  ; ...
          2    1   1  ; ...
          2   91   1  ; ...
          2  181   1  ; ...
          2  271   1  ; ...
          3    1   1  ; ...
          4    1   1  ; ...
          4    1  91  ; ...
          4    1 181  ; ...
          4    1 271  ; ...
          5    1   1  ; ...
          5    1  91  ; ...
          5    1 181  ; ...
          5    1 271  ];
for k = 1:size(iti,1)
  disp(sprintf('  k = %d',k));
  fname = sprintf('bathy.%03d.001.data',k);
  disp( fname )
  fid = fopen(fname, 'w', 'ieee-be');
  is = iti(k,3);  ie = is + 90 - 1;
  js = iti(k,2);  je = js + 90 - 1;
  tmp = ginfo(iti(k,1)).b_out(is:ie,js:je);
  fwrite(fid,tmp,'real*8',0,'ieee-be');
  fclose(fid);
  surf(tmp), view(2), axis equal, shading flat, pause(1)
end

1	edhill	1.1	%
2			% Ed Hill
3			%
4			% Generate approximate bathymetry for the llc grid
5			%
6
7			% clear all ; close all
8
9			do_print = 0;
10			dbug = 0;
11
12			% Get the ETOPO2 data
13			disp('Reading the ETOPO2 data...');
14			nlat = 5400;
15			nlon = 10800;
16			lons = linspace(-180, (180 - 2/60), nlon);
17			lats = linspace(90, -(90 - 2/60), nlat);
18			fid = fopen('ETOPO2.raw.bin', 'r', 'ieee-be');
19			et2 = reshape( fread(fid, nlat*nlon, 'int16'), [ nlon nlat ]);
20			fid = fclose(fid);
21
22			if dbug > 10
23			surf( et2(1:500,1:500) ), view(2), shading flat
24			ilon = [ 1:10:5000 ];
25			ilat = [ 1:10:5400 ];
26			surf(lons(ilon), lats(ilat), et2(ilon,ilat)'), view(2), ...
27			shading flat, axis equal
28			end
29
30			% Get the llc grid information
31			disp('Reading the grid information...');
32			fnc_in = 'llc_p90_%d.nc';
33			vnms = { 'XC' 'YC' 'XG' 'YG' };
34			d2r = pi/180.0;
35			ginfo = {};
36			for k = 1:5
37			nc = netcdf(sprintf(fnc_in,k),'nowrite');
38			for j = 1:length(vnms)
39			eval(sprintf('ginfo(k).%s = nc{''%s''}(:);',vnms{j},vnms{j}));
40			end
41			nc = close(nc);
42
43			% compute 3D coords of the C and G points
44			cor = zeros( [ size(ginfo(k).XG) 3 ]);
45			[cor(:,:,1),cor(:,:,2),cor(:,:,3)] = sph2cart( ginfo(k).XG*d2r, ...
46			ginfo(k).YG*d2r,1 );
47			cen = zeros( [ size(ginfo(k).XC) 3 ]);
48			[cen(:,:,1),cen(:,:,2),cen(:,:,3)] = sph2cart( ginfo(k).XC*d2r, ...
49			ginfo(k).YC*d2r,1 );
50			ginfo(k).cor = cor;
51			ginfo(k).cen = cen;
52
53			% empty bathy
54			ginfo(k).bathy = zeros(size( size(ginfo(k).XC) ));
55			end
56
57			% Do a quick-and-dirty regrid
58			disp('Performing the regrid...');
59	heimbach	1.5	%ph(
60			%%%s_lon = 0.0;
61			s_lon = -37.0;
62			%ph)
63			s_lat = 0.0;
64	edhill	1.1	for k = 1:5
65			disp(sprintf(' k = %d',k));
66			%
67			% / i,j+1 i+1,j+1 \
68			% \ i,j i+i,j /
69			for j = 1:size(ginfo(k).XC,2)
70			jr = [ j j+1 ];
71			for i = 1:size(ginfo(k).XC,1)
72			ir = [ i i+1 ];
73			lon_min = min(min( ginfo(k).XG(ir,jr) + s_lon ));
74			lon_max = max(max( ginfo(k).XG(ir,jr) + s_lon ));
75			lat_min = min(min( ginfo(k).YG(ir,jr) + s_lat ));
76			lat_max = max(max( ginfo(k).YG(ir,jr) + s_lat ));
77			% [ lon_min lon_max lat_min lat_max ]
78			lon_mm = [ mod(lon_min + 180,360)-180 ...
79			mod(lon_max + 180,360)-180 ];
80			lon_min = min(lon_mm);
81			lon_max = max(lon_mm);
82
83			d_lon = lon_max - lon_min;
84			if (abs(d_lon) > 45 && lat_min < 88)
85			i_lon = find( lons >= lon_max ...
86			\| lons <= lon_min );
87			else
88			i_lon = find(lon_min <= lons & lons <= lon_max);
89			end
90			i_lat = find(lat_min <= lats & lats <= lat_max);
91
92			% average the bathy
93			bv = et2(i_lon,i_lat);
94			if length(bv(:)) > 1
95			ginfo(k).bathy(i,j) = mean(bv(:));
96			else
97			% If not averaging, then interpolate
98			% ginfo(k).bathy(i,j) = NaN;
99			i_lon = interp1(lons,[1:length(lons)],lon_min,'nearest');
100			i_lat = interp1(lats,[1:length(lats)],lat_min,'nearest');
101			ginfo(k).bathy(i,j) = et2(i_lon,i_lat);
102			end
103			end
104			end
105			end
106
107			% Plot the resulting bathymetry
108			k = 3;
109			for k = 1:5
110			disp(sprintf(' k = %d',k));
111			ginfo(k).b_oce = ginfo(k).bathy;
112	edhill	1.2	% ginfo(k).b_oce(find(ginfo(k).b_oce >= 0.0)) = NaN;
113			ginfo(k).b_oce(find(ginfo(k).b_oce >= 0.0)) = 1000;
114			if k == 2
115			hold on
116			end
117			% figure(k)
118	edhill	1.1	surf( ginfo(k).cor(:,:,1), ...
119			ginfo(k).cor(:,:,2), ...
120			ginfo(k).cor(:,:,3), ginfo(k).b_oce )
121	edhill	1.2	if k == 5
122			hold off
123			end
124	edhill	1.1	end
125	edhill	1.2	axis equal, view(2)
126	edhill	1.1
127
128	edhill	1.3	% Plot the bathymetry at 1/2 resolution
129			k = 3;
130			for k = 1:5
131			disp(sprintf(' k = %d',k));
132			ginfo(k).c_half = ginfo(k).cor(1:2:end,1:2:end,:);
133			ginfo(k).b_half = ginfo(k).bathy(1:2:end,1:2:end);
134			ginfo(k).b_half(find(ginfo(k).b_half >= 0.0)) = 1000;
135			if k == 2
136			hold on
137			end
138			% figure(k)
139			surf( ginfo(k).c_half(:,:,1), ...
140			ginfo(k).c_half(:,:,2), ...
141			ginfo(k).c_half(:,:,3), ginfo(k).b_half )
142			if k == 5
143			hold off
144			end
145			end
146			axis equal, view(2)
147
148	heimbach	1.4	%--- BREAK POINT HERE
149			%--- do e.g. view(80,45) to change view angle
150
151
152	edhill	1.3
153	edhill	1.1	if do_print > 0
154			print -depsc llc_bathy.eps
155			end
156
157
158	edhill	1.2	% Write the bathymetry to netCDF compatible with MNC
159			id = 'bathy';
160			units = 'm';
161			for k = 1:5
162			disp(sprintf(' k = %d',k));
163			ginfo(k).b_out = ginfo(k).bathy;
164			ginfo(k).b_out(find(ginfo(k).bathy >= 0.0)) = 0.0;
165			bfile = sprintf('llc_p90_bathy.f%03d.nc',k);
166
167			nc = netcdf(bfile, 'clobber');
168			nc.reference = 'bathymetry';
169			nc.author = 'Ed Hill <eh3@mit.edu>';
170			nc.date = eval('date');
171			nc('X') = size( ginfo(k).XC, 1 );
172			nc('Y') = size( ginfo(k).XC, 2 );
173			nc{ id } = { 'Y', 'X' };
174			nc{ id }.units = units;
175			nc{ id }(:) = ginfo(k).b_out';
176			nc = close(nc);
177			end
178
179			% Write the bathymetry to a set of 90-by-90 per-tile MDSIO-compatible
180			% files
181			% face is js
182			iti = [ 1 1 1 ; ...
183			1 91 1 ; ...
184			1 181 1 ; ...
185			1 271 1 ; ...
186			2 1 1 ; ...
187			2 91 1 ; ...
188			2 181 1 ; ...
189			2 271 1 ; ...
190			3 1 1 ; ...
191			4 1 1 ; ...
192			4 1 91 ; ...
193			4 1 181 ; ...
194			4 1 271 ; ...
195			5 1 1 ; ...
196			5 1 91 ; ...
197			5 1 181 ; ...
198			5 1 271 ];
199			for k = 1:size(iti,1)
200			disp(sprintf(' k = %d',k));
201			fname = sprintf('bathy.%03d.001.data',k);
202	heimbach	1.4	disp( fname )
203	edhill	1.2	fid = fopen(fname, 'w', 'ieee-be');
204			is = iti(k,3); ie = is + 90 - 1;
205			js = iti(k,2); je = js + 90 - 1;
206	heimbach	1.4	tmp = ginfo(iti(k,1)).b_out(is:ie,js:je);
207	edhill	1.2	fwrite(fid,tmp,'real*8',0,'ieee-be');
208			fclose(fid);
209	heimbach	1.4	surf(tmp), view(2), axis equal, shading flat, pause(1)
210	edhill	1.2	end
211	edhill	1.1