Matlab读取多波段遥感数据，依据多波段信息建立额外栅格图层，处理后导出

Leon_124

已于 2024-01-29 21:45:12 修改

阅读量1.4k

点赞数 2

分类专栏： Matlab 文章标签： matlab

于 2022-09-01 22:31:55 首次发布

本文链接：https://blog.csdn.net/sdau_ly124/article/details/126652640

版权

Matlab 专栏收录该内容

14 篇文章 3 订阅

订阅专栏

问题：已有RGB三波段的DEM图，要求创造第四图层专门存放海拔信息

Matlab创造带地理信息图层

clear all
clc;
%% 读多层波段
filepath = 'D:\桌面\0825水位处理\地形地貌图_A3打印1.tif';
Info = imfinfo(filepath);

tif='tif';
format=Info.Format;
if  (strcmp(format ,tif)==0)
    disp('载入的不是tif图像，请确认载入的数据');                %%确保载入的图像是tiff图像
end
Slice=size(Info,1);                                          %%获取图片z向帧数
Width=Info.Width;
Height=Info.Height;
 
Image=zeros(Height,Width,Slice*3);
 
for i=1:Slice
    Image(:,:,(i-1)*3+1:i*3)=imread(filepath,i);                 %%一层一层的读入彩色图像
end

%% 单一波段
filepath2 = 'D:\0825水位处理\DEM提取\单一boduan.tif';
Info2 = imfinfo(filepath2);

tif='tif';
format2=Info2.Format;
if  (strcmp(format2 ,tif)==0)
    disp('载入的不是tif图像，请确认载入的数据');                %%确保载入的图像是tiff图像
end
Slice=size(Info2,1);                                          %%获取图片z向帧数
Width2=Info2.Width;
Height2=Info2.Height;
 
Image2=zeros(Height2,Width2,Slice*1);
 
for i=1:Slice
    Image2(:,:,(i-1)*1+1:i*1)=imread(filepath2,i);                 %%一层一层的读入彩色图像
end

%% 分离多层波段中每一波段
R_BODUAN = Image(:,:,1); %% R波段
G_BODUAN = Image(:,:,2); %% G波段
B_BODUAN = Image(:,:,3); %% B波段

%创建新矩阵
res = Image2;
res(:,:) = 0;
%根据指定规则给新矩阵每一点赋值 
res((R_BODUAN== 0) & (G_BODUAN == 0)&(B_BODUAN== 254)) = 0.5; %R=0 AND G=0 AND B=254 设置海拔为0.5
res((R_BODUAN==56)&(G_BODUAN==106)&(B_BODUAN==255))=1.25;
res((R_BODUAN==52)&(G_BODUAN==191)&(B_BODUAN==255))=1.75;
res((R_BODUAN==88)&(G_BODUAN==255)&(B_BODUAN==255))=2.25;
res((R_BODUAN==123)&(G_BODUAN==230)&(B_BODUAN==0))=2.75;
res((R_BODUAN==237)&(G_BODUAN==255)&(B_BODUAN==57))=3.25;
res((R_BODUAN==255)&(G_BODUAN==200)&(B_BODUAN==1))=3.75;
res((R_BODUAN==255)&(G_BODUAN==119)&(B_BODUAN==1))=4.5;
res((R_BODUAN==254)&(G_BODUAN==0)&(B_BODUAN==0))=5;
%% 写出图像
% 读取tif及空间坐标信息
%% 写出图像
% 读取tif及空间坐标信息
[A,R] = geotiffread(filepath2); %一定要把A写上，虽然后面用不到，但是不写会报错。原因就是错吧R认为A
info3 = geotiffinfo(filepath2);

% 保存tif并使用原始数据空间坐标信息
subR = R;
filename = 'Lixiahe.tif';
geotiffwrite(filename, res, subR,  ...
       'GeoKeyDirectoryTag', info3.GeoTIFFTags.GeoKeyDirectoryTag);
figure
mapshow(filename);