2011数模国赛A题部分Matlab代码

close all x=gd(:,1); y=gd(:,2); z=gd(:,3); scatter(x,y,5,z)%散点图 %figure [X,Y,Z]=griddata(x,y,z,linspace(1,30000)',linspace(1,20000),'v4');%插值 %pcolor(X,Y,Z);shading interp%伪彩色图 figure [C,h]=contourf(X,Y,Z) %等高线图; a=0:pi/36:2*pi; mymaker=['p';'o';'*';'s';'d'] hold on for i=1:size(gd,1) plot(gd(i,1),gd(i,2),mymaker(gd(i,4))) end colormap gray %plot(x,y,'*') hold on for i=1:size(gd,1) plot(50*gd(i,5)*cos(a)+gd(i,1),50*gd(i,5)*sin(a)+gd(i,2),'r-'); end colorbar

b={'As','Cd','Cr','Cu','Hg','Ni','Pb','Zn'};
load nd
for i=1:1
    Y=fix(nd(:,i));
    theta = [10 10]; lob = [1e-1 1e-1]; upb = [20 20];
    %[dmodel, perf] = dacefit(S, Y, @regpoly0, @corrspherical, theta, lob, upb) %球状模型
    [dmodel, perf] = dacefit(S, Y, @regpoly0, @correxp, theta, lob, upb)      %指数模型
    %>>>>>>>>>>>>>>predict<<<<<<<<<<<<<<<<<<<<<<<<<
    m=100;
    X = gridsamp([0 0;30000 20000], m);
    [YX MSE] = predictor(X, dmodel);
    X1 = reshape(X(:,1),m,m); X2 = reshape(X(:,2),m,m);
    YX = reshape(YX, size(X1));
    %figure(i), mesh(X1, X2, YX)
    figure(i), mesh(X1, X2, YX)
    hold on,
    plot3(S(:,1),S(:,2),Y,'.k', 'MarkerSize',10)
    %plot3(posi(:,1),posi(:,2),'*')    
    hold off
    xlabel('x/m')
    ylabel('y/m')
    zlabel('浓度')
    title([b(i)])
    colormap gray
    %print(gcf,'-dpng',['E:\360data\重要数据\桌面\空间分布\kriging插值图\' cell2mat(b(i)) '.png'])
    %print(gcf,'-djpeg',['E:\360data\重要数据\桌面\空间分布\等值线\' cell2mat(b(i)) '.jpg'])
end
%------------------------------------

% >>>>>>>>>>>>>>>> Kriging_fit.m ====================%S=data(:,1:2);Y=data(:,8);theta = [10 10]; lob = [1e-1 1e-1]; upb = [20 20];%[dmodel, perf] = dacefit(S, Y, @regpoly0, @corrspherical, theta, lob, upb)[dmodel, perf] = dacefit(S, Y, @regpoly0, @correxp, theta, lob, upb)%>>>>>>>>>>>>>>predict<<<<<<<<<<<<<<<<<<<<<<<<<m=80;X = gridsamp([0 0;30000 20000], m);[YX MSE] = predictor(X, dmodel);X1 = reshape(X(:,1),m,m); X2 = reshape(X(:,2),m,m);YX = reshape(YX, size(X1));figure(1), contourf(X1, X2, YX)hold on,plot3(S(:,1),S(:,2),Y,'.k', 'MarkerSize',10)hold offxlabel('x')ylabel('y')zlabel('浓度')title('Zn浓度空间分布')%===========

b={'生活区','工业区','山区','交通区','公园绿地区'} for i=1:5 figure(i) a=eval(['a' num2str(i)]); a=a'; bar(a); legend('背景值','观测值均值','国家二级标准') set(gca,'XTickLabel',['As';'Cd';'Cr';'Cu';'Hg';'Ni';'Pb';'Zn']) title([b(i)]) xlabel('x') ylabel('y') colormap gray set(gca,'fontsize',16) print(gcf,'-dpng',['E:\360data\重要数据\桌面\城市区域\' num2str(i) '.png']) end close all
===============================================