clear all
close all
%
%----------------------------------------------------------
ltoponame = 'Topo_MIKE21.txt';
gtoponame = '/Users/pmarches/Roms_tools/Topo/GEBCO_2014_SCS.nc';
otoponame = 'Topo_Mekong.nc';

lonmin=103;
lonmax=109;
latmin=7;
latmax=12;
%----------------------------------------------------------
%
addpath('/Users/pmarches/Roms_tools/Run_mekong/geodetic')
%
%  Read local topo (UTM)
%
BATHY=load(ltoponame);
E=BATHY(:,1); 
N=BATHY(:,2);
H=-BATHY(:,3);
%
% make a plot in UTM coordinates
%
%h=max(abs(H),0.1);
%figure
%scatter(E,N,h,h,'filled');
%colorbar
%
% Convert UTM to LAT,LON
%
Zone=48;
rad2deg=180/pi;
[LAT,LON]=utm2ell(N,E,Zone);
LAT=rad2deg*LAT;
LON=rad2deg*LON;
%
%  Open global topo file
%
nc=netcdf(gtoponame,'r');
tlon=nc{'lon'}(:);
tlat=nc{'lat'}(:);
dl=tlon(2,1)-tlon(1,1);
%
%  Find domain limits
%
%lonmin=min(LON)-dl;
%lonmax=max(LON)+dl;
%latmin=min(LAT)-dl;
%latmax=max(LAT)+dl;
%
%  get a subgrid
%
j=find(tlat>=latmin & tlat<=latmax);
i1=find(tlon-360>=lonmin & tlon-360<=lonmax);
i2=find(tlon>=lonmin & tlon<=lonmax);
i3=find(tlon+360>=lonmin & tlon+360<=lonmax);
x=cat(1,tlon(i1)-360,tlon(i2),tlon(i3)+360);
y=tlat(j);
[lon,lat]=meshgrid(x,y);
%
%  Read data
%
if ~isempty(i2)
  topo=-nc{'topo'}(j,i2);
else
  topo=[];
end
if ~isempty(i1)
  topo=cat(2,-nc{'topo'}(j,i1),topo);
end
if ~isempty(i3)
  topo=cat(2,topo,-nc{'topo'}(j,i3));
end
close(nc);
%
hg=topo; clear topo;
hg(hg>10)=hg(hg>10)+5;
%
%  In global topo, find land points and 
%  remove areas where local topo is available
%  Then, add local and global topo data
%
[M L]=size(lon);
HL  =999*ones(size(lon));
LONL=999*ones(size(lon));
LATL=999*ones(size(lon));
dlo2=dl/2;
% land area uncovered by local data
for j=2:M-1; for i=2:L-1;
 D=find(LON<lon(j,i+1)-dlo2 & LON>lon(j,i-1)+dlo2 & ...
        LAT<lat(j+1,i)-dlo2 & LAT>lat(j-1,i)+dlo2);
 if isempty(D) & hg(j,i)<1, 
  HL(j,i)=NaN;
  LONL(j,i)=lon(j,i);
  LATL(j,i)=lat(j,i);
 end;
end; end;
% sea area covered only by global data
for j=6:M-5; for i=6:L-5;
 D=find(LON<lon(j,i+5) & LON>lon(j,i-5) & ...
        LAT<lat(j+5,i) & LAT>lat(j-5,i));
 if isempty(D) & hg(j,i)>1, 
  HL(j,i)=hg(j,i);
  LONL(j,i)=lon(j,i);
  LATL(j,i)=lat(j,i);
 end;
end; end;
% remove points in common local and global areas
HL(HL==999)=[];
LONL(LONL==999)=[];
LATL(LATL==999)=[];
% concatenate local, land and global points
H=[H;HL'];
LON=[LON;LONL'];
LAT=[LAT;LATL'];
%
%  Interpolate composite topo 
%  on a high-resolution regular grid
%
X=[lonmin:1/1000:lonmax];
Y=[latmin:1/1000:latmax];
[lon,lat]=meshgrid(X,Y);
h=griddata(LON,LAT,H,lon,lat,'linear');
h(isnan(h))=-1;
%
%  Create and fill composite topo file
%
[M L]=size(lon);
nw = netcdf(otoponame, 'clobber');
nw('lat') = M;
nw('lon') = L;
nw{'lat'} = ncdouble('lat');
nw{'lat'}.long_name = ncchar('latitude');
nw{'lat'}.long_name = 'latitude';
nw{'lat'}.units = ncchar('degrees_north');
nw{'lat'}.units = 'degrees_north';
nw{'lon'} = ncdouble('lon');
nw{'lon'}.long_name = ncchar('longitude');
nw{'lon'}.long_name = 'longitude';
nw{'lon'}.units = ncchar('degrees_east');
nw{'lon'}.units = 'degrees_east';
nw{'topo'} = ncdouble('lat', 'lon');
nw{'topo'}.long_name = ncchar('height_above_reference_ellipsoid');
nw{'topo'}.long_name = 'height_above_reference_ellipsoid';
nw{'topo'}.units = ncchar('Meters');
nw{'topo'}.units = 'Meters';
nw.title = ncchar('Composite of local and global topo data');
nw.title = 'Composite of local and global topo data';
nw.date = ncchar(date);
nw.date = date;
nw{'lat'}(:)=Y;
nw{'lon'}(:)=X;
nw{'topo'}(:)=-h; % positive on land
close(nw)
%
%  Make a plot
%
domaxis=[min(min(lon)) max(max(lon))  ...
         min(min(lat)) max(max(lat))];
colaxis=[min(min(h)) max(max(h))];
figure
m_proj('mercator',...
       'lon',[domaxis(1) domaxis(2)],...
       'lat',[domaxis(3) domaxis(4)]);
m_pcolor(lon,lat,h);
shading flat
colorbar
m_gshhs_f('color','r')
m_grid('box','fancy',...
       'xtick',5,'ytick',5,'tickdir','out',...
       'fontsize',7);
figureHandle = gcf;
set(findall(figureHandle,'type','text'),'fontSize',14, ...
                                        'fontWeight','bold')