kriging.rar_kriging_kriging matlab_插值 面 matlab_球面

function [z, varargout] = kriging(xi, yi, zi, x, y, pl) % Perform interpolation using Ordinary kriging. % [z, z_var] = kriging(xi, yi, zi, x, y) % % [x, y, z] = kriging(...) % [x, y, z, z_var] = kriging(...) % [x, y, z, z_var, time] = kriging(...) % % Matt Foster <ee1mpf@bath.ac.uk> % Code heavily based on recipe from MATLAB recipes for Earth sciences. Martin % H. Trauth, Springer 2006 warning off all; tic; if nargin < 6 pl = 0; end % Create Variogram -- warning.. even this is slow! [X1,X2] = meshgrid(xi); [Y1,Y2] = meshgrid(yi); [Z1,Z2] = meshgrid(zi); D = sqrt((X1 - X2).^2 + (Y1 - Y2).^2); G = 0.5*(Z1 - Z2).^2; D2 = D.*(diag(xi*NaN)+1); lag = mean(min(D2)); hmd = max(D(:))/2; max_lags = floor(hmd/lag); LAGS = ceil(D/lag); for i = 1:max_lags SEL = (LAGS == i); DE(i) = mean(mean(D(SEL))); PN(i) = sum(sum(SEL == 1))/2; GE(i) = mean(mean(G(SEL))); end lags=0:max(DE); mod_func = fit_spherical(DE, GE, var(zi(:))); model = mod_func(D); if pl plot(DE, GE, '.'); hold on plot(DE, mod_func(DE)); hold off dlmwrite('variogram.dat', [DE', GE', mod_func(DE)'], '\t'); end % add zeros as the lat for and col of the variogram model. n = length(xi); model(:,n+1) = 1; model(n+1,:) = 1; % add a zero to the bottom right corner model(n+1,n+1) = 0; % Invert the variogram model model_inv = inv(model); Zg = nan(size(x)); s2_k = nan(size(x)); % Perform the Kriging for k = 1:length(Zg(:)); DOR = ((xi - x(k)).^2+(yi - y(k)).^2).^0.5; G_R = mod_func(DOR); G_R(n+1) = 1; E = model_inv*G_R; Zg(k) = sum(E(1:n,1).*zi); s2_k(k) = sum(E(1:n,1).*G_R(1:n,1))+E(n+1,1); end whos z z = Zg; z_var = s2_k; time = toc; if nargout == 2 || nargout == 3 varargout{1} = z_var; end if nargout == 3 varargout{2} = time; end if nargout >= 4 tmp = z; z = xx; varargout{1} = yy; varargout{2} = tmp; varargout{3} = z_var; end if nargout == 5 varargout{4} = time; end warning on all; end function mod_func = fit_spherical(DE, GE, var_z) % Only fit to the section below 90% of the variance ind = 1:min(find(GE > var_z * 0.9)); coef = [DE(ind);DE(ind).^3]' \ GE(ind)'; c = sqrt((-0.5 .* (1.5 ./ coef(1)).^-3)./coef(2)); a = 1.5 .* c / coef(1); %---------------------------------------------------- % lag_var = DE(min(find(GE > var_z))); % if a > lag_var; % a = lag_var; % end % % if c > var_z % c = mean(GE(min(find(GE > var_z)):end)); % end %---------------------------------------------------- if isreal(a) % Return an anonymous function that models the variogram mod_func = @(h) (coef(1) .* h + coef(2) .* h.^3) .* (h <= a) ... + c.*ones(size(h)) .* (h > a); else mod_func = @(h) (coef(1) .* h + coef(2) .* h.^3); end end