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

问题:已有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);

  • 去除背景值0

空间分析工具——栅格计算器——SetNull("raster"==100,"raster")

  • 矢量转point后,IDW插值处理

  • 2
    点赞
  • 8
    收藏
    觉得还不错? 一键收藏
  • 打赏
    打赏
  • 1
    评论
好的,以下是一个Python脚本可以完成您的要求: ```python import gdal import ogr import numpy as np # 打开多波段栅格数据 ds = gdal.Open("input_filename.tif") # 获取栅格数据的行列数和波段数 cols = ds.RasterXSize rows = ds.RasterYSize bands = ds.RasterCount # 遍历波段,在每个波段中提取值,并将其转换为点并添加到shp文件中 for band in range(bands): band_num = band + 1 band_data = ds.GetRasterBand(band_num).ReadAsArray() # 获取栅格的地理参照信息 transform = ds.GetGeoTransform() x_start = transform[0] y_start = transform[3] pixel_width = transform[1] pixel_height = transform[5] # 创建一个shp文件来存储提取的点 shp_driver = ogr.GetDriverByName("ESRI Shapefile") shp_ds = shp_driver.CreateDataSource("output_filename.shp") shp_layer = shp_ds.CreateLayer("point_layer", geom_type=ogr.wkbPoint) shp_layer.CreateField(ogr.FieldDefn("value", ogr.OFTReal)) # 遍历栅格数据中的每个像素,将其转换为点并添加到shp文件中 for i in range(rows): for j in range(cols): x = x_start + pixel_width * j y = y_start + pixel_height * i value = band_data[i][j] point = ogr.Geometry(ogr.wkbPoint) point.AddPoint(x, y) feature = ogr.Feature(shp_layer.GetLayerDefn()) feature.SetGeometry(point) feature.SetField("value", value) shp_layer.CreateFeature(feature) feature = None print("Band " + str(band_num) + " processed and added to shapefile") # 将shp文件中的点和值导出为excal文件 excal_driver = ogr.GetDriverByName("XLSX") excal_ds = excal_driver.CreateDataSource("output_filename.xlsx") excal_layer = excal_ds.CreateLayer("data", geom_type=ogr.wkbNone) excal_layer.CreateField(ogr.FieldDefn("x", ogr.OFTReal)) excal_layer.CreateField(ogr.FieldDefn("y", ogr.OFTReal)) excal_layer.CreateField(ogr.FieldDefn("value", ogr.OFTReal)) for feature in shp_layer: excal_feature = ogr.Feature(excal_layer.GetLayerDefn()) excal_feature.SetField("x", feature.GetGeometryRef().GetX()) excal_feature.SetField("y", feature.GetGeometryRef().GetY()) excal_feature.SetField("value", feature.GetField("value")) excal_layer.CreateFeature(excal_feature) excal_feature = None excal_ds = None print("All bands processed and data exported to excal file") ``` 您需要将`input_filename.tif`替换为您的多波段栅格数据文件名,将`output_filename`替换为您想要的输出文件名。请注意,此代码假定多波段栅格数据是以与GDAL兼容的格式存储的。

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 1
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

Leon_124

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付:¥1
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值