点云数据在Matlab生成格网数据-2021-10-18

博主10-31日的更新链接如下：就是原模原样的代码，但是把R空间参考改了，可以储存地理坐标或者投影坐标了

点云生成栅格数据补充+带有地理坐标，2021-10-31_~追风筝的猫的博客-CSDN博客

---

点云栅格化，即把点数据划分为若干个大小相等的格网单元，并将一定的属性值填充进入格网单元内，如高程值、强度值等，最终形成一个带有坐标系的影像。

---

按照点云中的xyz坐标值，直接判断放到相应的栅格中，即可实现点云栅格化。

博主很长一段时间研究格网中怎么放入坐标系，怎么储存带有坐标系的地理影像，但是都以失败告终，所以在这个博客中，主要是存放一段Matlab结合GIS生成栅格数据的代码。

---

值得说明的是：点本身是不具备物理意义的，因此博主在写的过程中，借助了R语言的参考，将原本的每个点，转换成了圆盘，这样每个圆盘就有实际的位置，有了物理意义，并且，当点稀疏且分辨率过高时，圆盘可以有效减少栅格化过程中产生的凹坑和空白线条。

---

%raster
% 栅格化点云数据
% 格网平均值、格网最大值、格网最小值
clear
[X,Y,Z] = sphere(100);
a = [X(:),Y(:),Z(:)];
pcshow(a)

% 提取xyz三维坐标
x = a(:,1);
y = a(:,2);
z = a(:,3);

% 圆盘大小为10cm
r = 0.1;
% 8个点的间距
p = 45;
% 创建新的集合
b = [];

% 读取点云的xy极限值
max_x = max(x);
min_x = min(x);
max_y = max(y);
min_y = min(y);

% 栅格化点云的分辨率
cells = 0.1;

% 计算格网的长度
lenx = ceil((max_x - min_x)/cells);
leny = ceil((max_y - min_y)/cells);
grid = zeros(leny,lenx);

% 制作圆盘
tic
temp1 = [a(:,1)+r*cosd(67.5),a(:,2)+r*sind(67.5),a(:,3)];
temp2 = [a(:,1)+r*cosd(22.5),a(:,2)+r*sind(22.5),a(:,3)];
temp3 = [a(:,1)-r*cosd(67.5),a(:,2)+r*sind(67.5),a(:,3)];
temp4 = [a(:,1)-r*cosd(22.5),a(:,2)+r*sind(22.5),a(:,3)];
temp5 = [a(:,1)+r*cosd(22.5),a(:,2)-r*sind(22.5),a(:,3)];
temp6 = [a(:,1)+r*cosd(67.5),a(:,2)-r*sind(67.5),a(:,3)];
temp7 = [a(:,1)-r*cosd(22.5),a(:,2)-r*sind(22.5),a(:,3)];
temp8 = [a(:,1)-r*cosd(67.5),a(:,2)-r*sind(67.5),a(:,3)];
b = [b;temp1;temp2;temp3;temp4;temp5;temp6;temp7;temp8;];
toc
clear temp1 temp2 temp3 temp4 temp5 temp6 temp7 temp8

tic
for j = 1:leny
    limity1 = (j-1) * cells + min_y;
    limity2 = j * cells + min_y;
    c = b(b(:,2)>=limity1 & b(:,2)<limity2,:);
    for i = 1:lenx
        limitx1 = (i-1) * cells + min_x;
        limitx2 = i * cells + min_x;
        d = c(c(:,1)>=limitx1 & c(:,1)<limitx2,:);
        if isempty(d) == 1
             grid(j,i) = NaN;
        else
            % 以最大值作为每个格网的值
            grid(j,i) = max(d(:,3));
        end
    end
end
toc
% grid = cell2mat(grid);
% 翻转格网
grid = flipud(grid);
imshow(grid)

% R是参数，reference，在arcgis中用ascii转栅格中使用
R = {
strcat('NCOLS',32,num2str(leny));
strcat('NROWS',32,num2str(lenx));
strcat('XLLCORNER',32,num2str(min_x));
strcat('YLLCORNER',32,num2str(min_y));
strcat('CELLSIZE',32,num2str(cells));
strcat('NODATA_VALUE',32,'NULL')
}

% 接下来使用save函数保存grid

产生数据如下：

栅格化结果：

---

打开保存的txt栅格数据，然后再排头加上R参数。比如下图，然后储存。

打开ArcGIS，搜索找到ASCII转栅格工具，输入数据为导出的txt数据，即可转换。

这里博主的惰性犯了，下面是GIS的帮助文档可看。

栅格转 ASCII—帮助 | ArcGIS for Desktop

不过这个转换好像有点玄学，博主试了很多次才成功。

---

如果一直是这样操作会很麻烦把，但是博主暂时没有更好的办法了，如果各位有建议，或者办法，可以给博主留言，一起讨论。