% GriddataILMag0305NEW.m  L Braile 3/17/05
% 
%
% Grid and contour Illinois magnetics data -- use interp1 (extrap) at edge to
% extrapolate the outer edge grid values
RandGridDataIllinoisMagnetics % obtain 200 random x-y
% locations and magnetics (z) values
dx = 6/19; % grid interval
[XI,YI] = meshgrid(0:dx:6); % define x-y grid locations
ZI = griddata(x,y,z,XI,YI,'cubic'); % grid (2-D interpolate)
%
% This section interpolates values at the edges of the map area 
% and within the map area where data points are sparse and the
% griddata code returns 'NaN' values.  This section assumes 
% that the grid area is slightly larger in each dimension than
% the x,y space of the input data.
%
% Find number of NaN values not on outer edge (NONE for this case)
[nx ny] = size(ZI); % nx is num of rows, ny is num of col
ny1 = ny-1;  nx1 = nx-1;
ic = 0;
for ix = 2:nx1
    for iy = 2:ny1
        iNaN = isnan(ZI(ix,iy));
        if iNaN == 1
            ic = ic +1;
        else
        end
    end
end
d = zeros(ic,1);
inx = zeros(ic,1);
iny = zeros(ic,1);
xm = mean(XI(:));
ym = mean(YI(:));
ic = 0;
for ix = 2:nx1
    for iy = 2:ny1
        iNaN = isnan(ZI(ix,iy));
        if iNaN == 1
            ic = ic +1;
            d(ic,1) = sqrt((XI(ix,iy) - xm)^2 + (YI(ix,iy) - ym)^2);
            inx(ic,1) = ix;
            iny(ic,1) = iy;
        else
        end
    end
end
A = [d inx iny];
% sort A so that NaN values are interpolated beginning with those
% closest to the center of the map and working outward
B = sortrows(A,1);
inx = B(:,2);
iny = B(:,3);
for i = 1:ic
    it = isnan(ZI(inx(i),:));
    xt = [];
    yt = [];
    zt = [];
    for iy = 2:ny1
        if it(iy) == 0
            xt = [xt XI(inx(i),iy)];
            zt = [zt ZI(inx(i),iy)];
        else
        end
    end
    ZI(inx(i),iny(i)) = interp1(xt,zt,XI(inx(i),iny(i)),'spline','extrap');
end
% Extrapolate to outer edge of map. Use 'nearest' for this case
% left and right edge
for ix = 2:nx1
    ZI(ix,1) = interp1(XI(ix,2:ny1),ZI(ix,2:ny1),XI(ix,1),'nearest','extrap');
    ZI(ix,ny) = interp1(XI(ix,2:ny1),ZI(ix,2:ny1),XI(ix,ny),'nearest','extrap');
end
% top and bottom and corners
for iy = 1:ny
    ZI(1,iy) = interp1(YI(2:nx1,iy),ZI(2:nx1,iy),YI(1,iy),'nearest','extrap');
    ZI(nx,iy) = interp1(YI(2:nx1,iy),ZI(2:nx1,iy),YI(ny,iy),'nearest','extrap');
end
%
figure
mesh(XI,YI,ZI);, hold on % Make a mesh of gridded surface
plot3(x,y,z,'o'), hold off % Show original x,y,z data values
set(gca,'fontsize',16,'linewidth',2)
xlabel('x','fontsize',16)
ylabel('y','fontsize',16)
title('200 random x,y data, Illinois Magnetics, Mesh','fontsize',16)
%
% Plot color image and overlay contours
figure
pcolor(XI,YI,ZI);
shading interp
colormap jet
colorbar
hold on
% plot the x-y locations of the original data
plot(x,y,'+k')
% Contour the gridded data
v = [-200:100:800];
[c,h] = contour(XI,YI,ZI,v,'w-');
clabel(c,h,'color','w','fontweight','bold','fontsize',15,'LabelSpacing',400)
axis square
set(gca,'fontsize',16,'linewidth',2)
xlabel('x','fontsize',16)
ylabel('y','fontsize',16)
title('200 random x,y data, Illinois Magnetics, coarse grid','fontsize',16)
hold off
%
% use interp2 to regrid for increased resolution and smoothness
figure
dx2 = 1/19;
[X2,Y2] = meshgrid(0:dx2:6); % set finer grid interval, grid = 0.05
Z2 = interp2(XI,YI,ZI,X2,Y2,'cubic');
% replot color display and contour map for interpolated data
pcolor(X2,Y2,Z2)
shading interp
colorbar
hold on
% Contour the gridded data
v = [-200:100:800];
[c,h] = contour(X2,Y2,Z2,v,'w-');
clabel(c,h,'color','w','fontweight','bold','fontsize',15,'LabelSpacing',400)
axis square
set(gca,'fontsize',16,'linewidth',2)
xlabel('x','fontsize',16)
ylabel('y','fontsize',16)
title('Illinois Magnetics, fine grid','fontsize',16)
hold off
%
% Read and plot the hand digitized data and the residual map
HandDigitizedGrid
figure
dx2 = 1/19;
[X2,Y2] = meshgrid(0:dx2:6); % set finer grid interval, grid = 0.05
Z3 = interp2(XI,YI,ZHD,X2,Y2,'cubic');
% replot color display and contour map for interpolated data
pcolor(X2,Y2,Z3)
shading interp
colorbar
hold on
% Contour the gridded data
v = [-200:100:800];
[c,h] = contour(X2,Y2,Z3,v,'w-');
clabel(c,h,'color','w','fontweight','bold','fontsize',15,'LabelSpacing',400)
axis square
set(gca,'fontsize',16,'linewidth',2)
xlabel('x','fontsize',16)
ylabel('y','fontsize',16)
title('Illinois Magnetics, Hand digitized, fine grid','fontsize',16)
hold off
% Residual surface (use fine grids)
Zdiff = Z2 - Z3;
% replot color display and contour map for interpolated data
pcolor(X2,Y2,Zdiff)
shading interp
colorbar
hold on
% Contour the gridded residual data
v = [-100 -50 -20 -10 10 20 50 100];
[c,h] = contour(X2,Y2,Zdiff,v,'k-');
clabel(c,h,'color','k','fontweight','bold','fontsize',15,'LabelSpacing',400)
axis square
% Plot data locations
plot(x,y,'+k')
set(gca,'fontsize',16,'linewidth',2)
xlabel('x','fontsize',16)
ylabel('y','fontsize',16)
title('Illinois Magnetics, Residual Surface, fine grid','fontsize',16)
hold off