点云生成栅格数据补充+带有地理坐标，2021-10-31

博主在18号的时候写下了点云格网化的代码，当时还没想到如何补充地理信息，所以栅格数据仅仅能借助ArcGIS完成，最近博主学习了geotiff的相关内容，于是做了一点补充，关于点云格网化的网址如下：

点云数据在Matlab生成格网数据-2021-10-18_~追风筝的猫的博客-CSDN博客

---

 代码和格网数据一样，只是空间参考R做出了改变。

%raster
% 栅格化点云数据
% 格网平均值、格网最大值、格网最小值
clear
a = randn(10000,3);
pcshow(a)


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

% 读取点云的xy极限值
max_x = max(a(:,1));
min_x = min(a(:,1));
max_y = max(a(:,2));
min_y = min(a(:,2));

% 栅格化点云的分辨率
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 = flipud(grid);

lat = [min_x,max_x];
lon = [min_y,max_y];
R = maprefcells(lat,lon ,size(grid));
% map开头为投影坐标，geo开头为地理坐标
% 32649为WSG84 49N带的坐标代码
geotiffwrite('test2.tif',grid,R,'CoordRefSysCode',32649)

结果就是把数据保存下来了。

该部分的代码只要稍加修改，就能成为function函数使用拉。

博主就不上传资源了。