????????以2021年美赛MCM赛题B题为例,对于这个题目分析的后续算法暂且不论,我们需要将问题可视化,即将澳大利亚各个地点的热量信息化成一张分布图。数据可以参考以下网址:????????Fires from Space: Australia | Kagglehttps://www.kaggle.com/carlosparadis/fires-from-space-australia-and-new-zeland????????以下是我们得到的信息状态:
?
? ? ? ? ?很明显信息状态分为四个值:经度、纬度、亮度、置信度。我们需要对这些信息进行处理,经纬度转化成x-y(以1km为度量)坐标系,亮度与置信度想乘表示目标点的值。
? ? ? ? 处理过程我也并不是处理的很好,如下:
clear;clc%%清空工作区
load modis_20200104%%导入澳大利亚2020年1月4日高程数据
x=modis_20200104(:,1);%纬度
y=modis_20200104(:,2);%经度
h=modis_20200104(:,3).*modis_20200104(:,4)/100;%热度
xx=(y-min(y(:)))*111;%%经纬度转km
yy=(x-min(x(:)))*111;
X=int32(xx)+1;%%保留到整数
Y=int32(yy)+1;
? ? ? ? 关于经纬度的算法比较复杂,大家可以自行搜索相关资料,以下展示分布图的绘制方法:
%%绘制这张图的核心在于离散的点到一格一格的像素之间的转化,比如可能两个点非常接近,我们就需要将他们进行一个平均处理,使得最终图片的整体呈现效果较好即可
[lx,ly]=meshgrid(0:1:max(X(:)),0:1:max(Y(:)));%%生成网格用于装我们要绘制的图形的值
Z=zeros(max(Y(:))+1,max(X(:))+1);%%Z用于存储每个网格中的值
rep=zeros(max(Y(:))+1,max(X(:))+1);%%rep用于存储每个网格中点的个数,我们最后通过Z./rep得到每个网格中的热度平均值
for i=1:size(X,1)
rep(Y(i),X(i))=rep(Y(i),X(i))+1;
Z(Y(i),X(i))=Z(Y(i),X(i))+h(i);
end
rep=rep+eps;%%eps表示一个特别特别小的数,防止某个网格没有对应的点,0/0无法计算
Z=Z./rep;
[row,col]=find(Z>0);%%查找大于0的点(即有值的点),row存储行数组,col存储列数组,row(i)与col(i)一一对应
for i=1:size(row)
high(i,1)=Z(row(i),col(i)) ;%%赋值准备绘图
end
figure(1);
surf(lx,ly,Z);%%绘制surf图像
axis equal%%让坐标系刻度均匀
shading interp%%光滑
title('Fire Distribution(km)')%%标题
hold on
????????但是这个代码仍有一点不足,在于其数据的点与点之间实在是相差太远(题目要分析维多利亚州东部,但是数据包含了整个澳大利亚),实际呈现效果如下(需要点一下上面那个立方体旋转的图标,然后右键图像,点击转至X-Y视图):
?????????很明显,只有右下角那一部分使我们需要的,那么我们就需要对我们的代码进行改进,以让其只显示我们制定的部分,改进部分代码比较简单,如下:
x_start=3380;%%表示x轴搜索的起点
x_len=350;%%表示x轴长度
y_start=400;%%表示y轴搜索的起点
y_len=250;%%表示y轴的长度
[lx,ly]=meshgrid(0:1:x_len,0:1:y_len);%%网格
Z=zeros(y_len+1,x_len+1);%%值
rep=zeros(y_len+1,x_len+1);%%点数
for i=1:size(X,1)%%对每个点,只要要求范围内的
if(X(i)>x_start&X(i)<(x_start+x_len)&Y(i)>y_start&Y(i)<(y_start+y_len))
rep(Y(i)-y_start,X(i)-x_start)=rep(Y(i)-y_start,X(i)-x_start)+1;
Z(Y(i)-y_start,X(i)-x_start)=Z(Y(i)-y_start,X(i)-x_start)+h(i);
end
end
rep=rep+eps;%%加很小的数
Z=Z./rep;
[row,col]=find(Z>0);%%查找
for i=1:size(row)
high(i,1)=Z(row(i),col(i)) ;%%赋值
end
figure(1);
surf(lx,ly,Z);%%绘图
axis equal%%坐标系同刻度
shading interp%%平滑
title('Fire Distribution(km)')
hold on
????????得到的图形如下(需要点一下上面那个立方体旋转的图标,然后右键图像,点击转至X-Y视图):
????????这就很有内味了。
? ? ? ? 对数据的处理仍有很大不足,这个例子只是为了让我们知道绘制此类图像的方法。
|