% do_gom.m
%
% Program to create 1 deg by 1 deg squares of data. Squares will be pasted (do_paste.m)
% to produce a 3 sec bathy grid of the Gulf of Maine Region (GOM).
%
% Program calls: blockmedian, surface, and grdcut from Generic Mapping
% Tools.
% 
% Input: bathy data
%
% Output: 1 degree by 1 degree squares.
%
% Erin R. Twomey
% etwomey@usgs.gov
% 9/11/08
% 
% Revised by Rich P. Signell 10/26/11
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Create the 1 deg by 1 deg plus 10% squares.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Define the longitude and latitude vectors.

lonc = [-71.5:1:-63.5]; %cutting longitude min.
loncm = [-70.5:1:-63.5 -63.0]; %cutting longitude max.
latc = [39.5 40:1:45];
latcm = [40.0:1:46.0];
row = {'A','B', 'C', 'D', 'E', 'F', 'G'};

% Define non-changing variables.

%res = '3c';
res_sec = 3;
res = [num2str(res_sec) 'c'];
dx = res_sec/2;

tension = 0.35;
convergence = 0.1;
aspect_ratio = 0.75;

tic
%for i = 1:9 % Longitude: 1-9; 9 columns
for i =1:9 
    for j = 1:7 % Latitude: 1-7; 7 rows
           
i
j

%lonlat = sprintf('%.1f/%.1f/%.1f/%.1f',lon_min(i),lon_max(i),lat_min(j),lat_max(j))
cut = sprintf('%.11f/%.11f/%.11f/%.11f',lonc(i)+dx/3600,loncm(i)-dx/3600,latc(j)+dx/3600,latcm(j)-dx/3600)

lonlat = sprintf('%.11f/%.11f/%.11f/%.11f',lonc(i)-.1+dx/3600,loncm(i)+.1+dx/3600,latc(j)-.1+dx/3600,latcm(j)+.1+dx/3600)
%Create the blockmedian command. 

lon_min=num2str((lonc(i)-.1+dx/3600),11);
lon_max=num2str((loncm(i)+.1-dx/3600),11);
lat_min=num2str((latc(j)-.1+dx/3600),11);
lat_max=num2str((latcm(j)+.1-dx/3600),11);

%cmd = sprintf('!cat final_xyz/*.xyz | awk '' $1>=%s && $1<=%s && $2>=%s && $2<=%s {printf " %%.6f %%.6f %%.2f %%d %%s\\n",$1,$2,$3,$4,$5}'' > awk_%s%d.xyz',...
%    lon_min,lon_max,lat_min,lat_max,char(row(j)),i)
outfile_xyz=sprintf('awk_%s%d.xyz',char(row(j)),i)

!rm *.xyz
cmd = sprintf('!for file in final_xyz/*.xyz;do cat $file | awk '' $1>=%s && $1<=%s && $2>=%s && $2<=%s {printf " %%.6f %%.6f %%.2f %%d %%s\\n",$1,$2,$3,$4,$5}'' >> awk_%s%d.xyz; done',...
    lon_min,lon_max,lat_min,lat_max,char(row(j)),i)
eval(cmd)


blockmedian_in= outfile_xyz;
%blockmedian_in='awk_D1r.xyz';
blockmedian_out=sprintf('GOM_061311_%s%d_3sec.xyz',char(row(j)),i);
cmd1 = sprintf('!blockmedian %s -R%s -I%s -C -V > %s',blockmedian_in,lonlat,res,blockmedian_out)

surface_out=sprintf('GOM_061311_%s%d_3sec.grd',char(row(j)),i);

%Create the new surface command.

cmd2 = sprintf('!surface %s -R%s -G%s -C%d -I%s -A%d -T%d -V',blockmedian_out,lonlat,surface_out,convergence,res,aspect_ratio,tension)

% Cut out the degree square 
cut_out=sprintf('GOM_061311_%s%d_3sec_cut.grd',char(row(j)),i);

cmd3 = sprintf('!grdcut %s -R%s -G%s -V',surface_out,cut,cut_out)

% Convert gridline registration to pixel registration
cmd4 = sprintf('!grdedit -T %s',cut_out')

% evaluate the commands

eval(cmd1);
 eval(cmd2);
 eval(cmd3);
 eval(cmd4);
 
toc
%eval(cmd3)
% Move the blockmedianed xyz file to another directory.

%c = filename;
%d = sprintf('G:\\GOM\\GOM_blockfiles_091208');
%movefile(c,d);

%movefile(cmd(164:173),'G:/GOM/deg_square_files')
%movefile(blockmedian_out,'G:/GOM/GOM_112008')

movefile(outfile_xyz,'../GOM_20110613')
movefile(blockmedian_out,'../GOM_20110613')

    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Merge the cut grids together to get the final image.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%cmd4 = sprintf('!gdal_merge.bat -o GOM_final_3sec.grd *_cut.grd')
%eval(cmd4)
