ENVI文件在matlab中的读取

最新推荐文章于 2025-02-21 15:57:17 发布

xiaowei9696

最新推荐文章于 2025-02-21 15:57:17 发布

阅读量1.1w

点赞数 7

分类专栏： MATLAB 遥感图像处理 IDL/ENVI

本文链接：https://blog.youkuaiyun.com/xiaowei9696/article/details/81146293

版权

遥感图像处理同时被 3 个专栏收录

6 篇文章

订阅专栏

IDL/ENVI

3 篇文章

订阅专栏

MATLAB

2 篇文章

订阅专栏

2018.07.21

皇天在上，终于找到读取的方法了，研究可以继续深入往下走了。

（1）tiff文件的读取。使用的MATLAB版本为2015b

[A,R] = geotiffread(filename)

http://ww2.mathworks.cn/help/map/ref/geotiffread.html?searchHighlight=geotiffread&s_tid=doc_srchtitle

A的格式与tiff原文件中的格式相同
R为reference，格式为MapCellsReference，解释见http://ww2.mathworks.cn/help/map/ref/map.rasterref.mapcellsreference.html?s_tid=doc_ta

ref = 

  MapCellsReference with properties:

            XWorldLimits: [659685 885615]   %用一个正方形代表像素，指的是像素顶点上的坐标范围
            YWorldLimits: [3240285 3471015] %与x同义
              RasterSize: [7691 7531]  %影像像素数
    RasterInterpretation: 'cells'   %即表示像素的值在cells中间而不在边缘上
        ColumnsStartFrom: 'north'
           RowsStartFrom: 'west'
      CellExtentInWorldX: 30    %像元在实际世界中的大小，单位为meter
      CellExtentInWorldY: 30
    RasterExtentInWorldX: 225930    %整个影像在实际世界中的大小，单位为meter
    RasterExtentInWorldY: 230730
        XIntrinsicLimits: [0.5 7531.5] %感觉是指代像元中心的意思，文档这部分读不太懂
        YIntrinsicLimits: [0.5 7691.5]
      TransformationType: 'rectilinear'
    CoordinateSystemType: 'planar'

info = geotiffinfo(filename)

http://ww2.mathworks.cn/help/map/ref/geotiffread.html?searchHighlight=geotiffread&s_tid=doc_srchtitle

info是一个struct，里面包含的信息量很多

info=geotiffinfo(fn)

info = 

             Filename: 'E:\Users\DATA\2015_Landsat7&Landsat8_SR\LC08…'
          FileModDate: '05-7月-2018 19:48:49'
             FileSize: 116057771
               Format: 'tif'
        FormatVersion: []
               Height: 7691
                Width: 7541
             BitDepth: 16
            ColorType: 'grayscale'
            ModelType: 'ModelTypeProjected'
                  PCS: 'WGS 84 / UTM zone 49N'
           Projection: 'UTM zone 49N'
               MapSys: 'UTM_NORTH'
                 Zone: 49
         CTProjection: 'CT_TransverseMercator'
             ProjParm: [7x1 double]
           ProjParmId: {7x1 cell}
                  GCS: 'WGS 84'
                Datum: 'World Geodetic System 1984'
            Ellipsoid: 'WGS 84'
            SemiMajor: 6378137
            SemiMinor: 6.3568e+06
                   PM: 'Greenwich'
    PMLongToGreenwich: 0
            UOMLength: 'metre'
    UOMLengthInMeters: 1
             UOMAngle: 'degree'
    UOMAngleInDegrees: 1
            TiePoints: [1x1 struct]
           PixelScale: [3x1 double]
           SpatialRef: [1x1 map.rasterref.MapCellsReference]
            RefMatrix: [3x2 double]
          BoundingBox: [2x2 double]
         CornerCoords: [1x1 struct]
         GeoTIFFCodes: [1x1 struct]
          GeoTIFFTags: [1x1 struct]

[A,R] = geotiffread(filename)中的R即为info.SpatiRef

（2）shp文件的读取

ENVI中的样本文件xml不能被matlab使用，那种转到ENVI Classic的方法也不需要再用了。直接转出为shapefile即可，不需要通过arcmap。

ENVI导出为shp文件：打开ROI工具，File>>Export>>Export to Shapefile

[S,A] = shaperead(shapefilename)；

http://ww2.mathworks.cn/help/map/ref/shaperead.html?s_tid=doc_ta

S为一个结构体，里面包含记录个数那么多数量的结构体，存储shp中记录的几何信息

shp(1)

ans = 

       Geometry: 'Polygon'    %结构体类型
    BoundingBox: [2x2 double]    %包围盒坐标
              X: [805095 805095 804945 804945 805095 NaN]    
           %顺序为：右下 > 左下 > 左上 > 右上 > 左下 > 标志结束 为包围盒中polygon的顶点连接方式
              Y: [3446475 3446325 3446325 3446475 3446475 NaN]

shp(1).BoundingBox %第一行为左上，第二行为右下

ans =

      804945     3446325
      805095     3446475

A为与S相同数量的结构体，包括shp中每个记录的其他信息

shp_info(1)

ans = 

    CLASS_NAME: 'barren'  %记录名  
      CLASS_ID: '1'    %记录的ID，跟ROI转换成shp文件时选择的记录类别有关
    CLASS_CLRS: '255,255,29'    %显示的颜色

（3）tiff的显示

可以使用imshow(img)直接显示，由于int16的数据类型，显示效果不佳，如果使用单纯的线性缩放，影像部分全是白，因为ENVI在对填充区域的像素值为-9999，而影像位深为2^14位，线性拉伸就白瞎了。
使用mapshow要注意的地方，①格式必须为uint8，②不能只用影像，还要包括ref

>>figure
>> mapshow(img_test,ref)
>> axis image off

2018.07.23

从师姐那里学到了另一种方法

①在ENVI中，ROI Tools》File》Export》Export to Classic

②在ENVI Classic中，先打开图像

③在图像窗口中，Overlay》Region of Interest 打开ROI Tool

④删除掉ROI Tool中default的第一类，然后Options》Create Class Image from ROIS

⑤得到一个*.hdr和一个同名文件，这个部分的读取方式就又要照着nfreadenvi来做了

⑥优点是相当于一张图上就有分类了，位置和分类信息同时获取，效率是否比自己摸索的方法快就先不管了

2018.08.06

对于读取的结果的一些注解，以下分析基于我在Landsat影像中勾画的5*5的矩形样本

①shpread之后的shape中的BoundingBox的两个点分别是左下、右上

② ref中的XWorldLimits、YWorldLimits都是按照大小顺序排列的，坐标系是大地坐标系（右手系），x是东西向，y是南北向

③第一行第一列对应的是（row，col）=（ ref.YWorldLimits(2)，ref.XWorldLimits(1) ）

prefix = 'E:\WeiLai\DATA\2015_Landsat7&Landsat8_SR\LC081230392015102501T1-SC20180705064849\LC08_L1TP_123039_20151025_20170402_01_T1_sr_band';
suffix = '.tif';
shpn = 'E:\WeiLai\DATA\sample\shp\summary.shp';
fn = strcat(prefix,int2str(1),suffix);
info = geotiffinfo(fn);
ref = info.SpatialRef;
row = ref.RasterSize(1);col = ref.RasterSize(2);
img = zeros(row,col,7);
% 直接读取原始影像
for i = 1:7
   fn = strcat(prefix,int2str(i),suffix);
   temp = geotiffread(fn);
   temp = double(temp);
   img(:,:,i) = temp;
end
clear temp

% read shapefile
[shp,shp_info] = shaperead(shpn);
% define the start position  第一行第一列对应于  ref.YWorldLimits(2)  -  ref.XWorldLimits(1)
%BoundingBox 中的两个点分别是 左下 - 右上
start(1) = ref.YWorldLimits(2);start(2) = ref.XWorldLimits(1); 

%采集样本
%在一块中随机选择像素作为样本
num_shp = size(shp,1);
sample = zeros(num_shp,7);
sample_label = zeros(num_shp,1);
pos = unidrnd(4,num_shp,2);
coordinate = zeros(num_shp,2); %左下角点的行列号
%选取样本 设置标签
for i = 1:num_shp
    coordinate(i,1) = ( start(1)-shp(i).BoundingBox(1,2) )/30; %row
    coordinate(i,2) =  ( shp(i).BoundingBox(1,1)-start(2) )/30; %col
    sample_label(i) = str2double(shp_info(i).CLASS_ID);
    sample(i,1:7)  =  img( coordinate(i,1)-pos(i,1) , coordinate(i,2)+pos(i,2) , : );
end