导出图像和栅格到GeoTIFF

打开脚本

此示例演示如何将引用到标准地理和投影坐标系的数据写入GeoTIFF文件，使用geotiffwrite．标记图像文件格式(TIFF)已成为存储栅格数据的流行格式。GeoTIFF规范定义了一组TIFF标记，用于描述与TIFF光栅数据相关的“制图”信息。使用这些标记，可以将地理定位图像或栅格网格与引用地理坐标系(纬度和经度)或(平面)投影坐标系的坐标存储在GeoTIFF文件中。

设置:定义一个数据文件夹和文件名实用函数

这个示例创建几个临时GeoTIFF文件并使用该变量datadir表示它们的位置。类的输出决定此处使用的值tempdir命令，但是您可以很容易地自定义它。的内容datadir在示例的末尾删除。

Datadir = fullfile(tempdir，“datadir”）;如果~存在(datadir“dir”mkdir (datadir)结束

定义一个匿名函数来前置datadir到输入文件名:

Datafile = @(filename)fullfile(datadir, filename);

例1:写一个参考地理坐标的图像

将引用WGS84地理坐标的图像写入GeoTIFF文件。原始图像(boston_ovr.jpg)以JPEG格式存储，在“世界文件”(boston_ovr.jgw)中包含图像文件外部的引用信息。该图像提供了马萨诸塞州波士顿及其周边地区的低分辨率“概览”。

从JPEG文件中读取图像。

: =“boston_ovr”；Imagefile = [basename .“jpg”];A1 = imread(imagefile);

从world文件中获取引用对象。

Worldfile = getworldfilename(图像文件);R1 = worldfileread(worldfile，“地理”、尺寸(A1));

将图像写入GeoTIFF文件。

Filename1 = datafile([basename .“.tif”]);geotiffwrite (filename1 A1, R1)

返回有关文件的信息RasterInfo对象。注意的值CoordinateReferenceSystem是地理坐标引用系统对象。

Info1 = georasterinfo(fileame1);info1。CoordinateReferenceSystem

ans = geoocrs与属性:名称:“WGS 84”数据:“世界大地测量系统1984”球体:[1×1 referenceEllipsoid]质点:0角度单位:“度”

重新导入新的GeoTIFF文件，并在地图上显示位置正确的波士顿概况图像。

figure usamap(R1.LatitudeLimits, r1 . longitude elimits) setm(gca，“PLabelLocation”, 0.05,“PlabelRound”2,“PlineLocation”，0.05) geoshow(filename1) title(“波士顿概述”）

例2:编写一个引用地理坐标的数据网格

加载高程栅格数据和地理单元格引用对象。将数据网格写入GeoTIFF文件。

负载topo60cZ2 = topo60c;R2 = topo60cR;Filename2 =数据文件(“topo60c.tif”）;geotiffwrite (filename2 Z2, R2)

数据网格中的值范围为-7473 ~ 5731。将网格显示为纹理映射的表面，而不是强度图像。

图worldmap世界gridm从setm (gca),“MLabelParallel”, -90,“MLabelLocation”，90) tmap = geoshow(文件名2，“DisplayType”，“texturemap”）;demcmap (tmap.CData)标题(“高程数据网格”）

例3:更改GeoTIFF文件的数据组织

写入数据时使用geotiffwrite或使用readgeoraster，数据网格的列从北开始，行从西开始。例如，来自的输入数据topo60c.mat起点从南，但输出数据从南topo60c.tif从北方出发。

R2。ColumnsStartFrom [Z3,R3] = readgeoraster(文件名2);R3。ColumnsStartFrom

Ans = '南' Ans = '北'

因此，输入数据和GeoTIFF文件中的数据是颠倒的。

isequal (Z2 flipud (Z3))

Ans =逻辑1

如果你需要工作空间中的数据再次匹配，那么翻转Z值并设置引用对象，以便列从南开始:

R3。ColumnsStartFrom =“南”；Z3 = flipd (Z3);Z3 isequal (Z2)

Ans =逻辑1

GeoTIFF文件中的数据是用正比例值编码的。因此，当您使用普通的tiff查看软件查看文件时，数据集的北部边缘位于顶部。要使输出文件中的数据布局与输入的数据布局相匹配，可以构造一个Tiff对象并使用它来重置一些标记和图像数据。

t = Tiff(文件名2，' r + '）;pixelScale = getTag(t，“ModelPixelScaleTag”）;pixelScale(2) = -pixelScale(2);setTag (t)“ModelPixelScaleTag”, pixelScale);tiepoint = getTag(t，“ModelTiepointTag”）;tiepoint(5) = intrinsic geography (R2,0.5,0.5);setTag (t)“ModelTiepointTag”, tiepoint);setTag (t)“压缩”， tiff . compress . none) write(t,Z2);rewriteDirectory (t)关闭(t)

验证来自输入数据的引用对象和数据网格与输出数据文件是否匹配。为此，读取Tiff图像并创建一个引用对象。然后，比较网格。

t = Tiff(文件名2);Atiff = read(t);close(t) Rtiff = georefcells(R2.LatitudeLimits, r2 . longitude elimits,size(Atiff));Atiff isequal (Z2) isequal (R2, Rtiff)

Ans = logical 1 Ans = logical 1

例4:写一个参考投影坐标系的图像

将协和正射影像写入一个GeoTIFF文件。这两张相邻的(从西到东)地理参考灰度(全色)正射影像覆盖了美国马萨诸塞州康科德的部分地区。concord_or .txt文件表示数据引用了马萨诸塞州大陆(NAD83)州平面投影坐标系统。单位是米。这对应于GeoTIFF规范中的GeoTIFF代码编号26986http://geotiff.maptools.org/spec/geotiff6.html#6.3.3.1PCS_NAD83_Massachusetts之下。

读这两个正射影像。

[X_west,R_west] = readgeoraster(“concord_ortho_w.tif”）;[X_east,R_east] = readgeoraster(“concord_ortho_e.tif”）;

组合图像和引用对象。

X4 = [X_west X_east];R4 = R_west;R4。XWorldLimits = [R_west.XWorldLimits(1) R_east.XWorldLimits(2)]; R4.RasterSize = size(X4);

将数据写入GeoTIFF文件。使用代号26986，表示PCS_NAD83_Massachusetts投影坐标系统。

coordRefSysCode = 26986;Filename4 =数据文件(“concord_ortho.tif”）;geotiffwrite (filename4 X4、R4、“CoordRefSysCode”coordRefSysCode)

返回有关文件的信息RasterInfo对象。注意的值CoordinateReferenceSystem投影坐标引用系统对象。

Info4 = georasterinfo(文件名4);info4。CoordinateReferenceSystem

ans =项目与属性:名称:“NAD83 /马萨诸塞州大陆”地理crs: [1×1 geocrs]投影方法:“兰伯特圆锥共形(2SP)”长度单元:"meter"投影参数:[1×1 map.crs.ProjectionParameters]

显示合并的协和正射影像。

图mapshow(filename4)“联合正色摄影”)包含(“MA大陆国家飞机东，米”) ylabel (“MA大陆国家飞机北，米”）

例5:从GeoTIFF文件中写入裁剪图像

将一个GeoTIFF文件的子集写入一个新的GeoTIFF文件。

中的RGB图像和引用信息boston.tifGeoTIFF文件。

[A5,R5] = readgeoraster(“boston.tif”）;

裁剪图像。

Xlimits = [764318 767677];Ylimits = [2951122 2954482];[A5crop,R5crop] = mapcrop(A5,R5,xlimits,ylimits);

将裁剪后的图像写入GeoTIFF文件。使用原始GeoTIFF文件中的GeoKeyDirectoryTag。

Info5 = geotiffinfo(“boston.tif”）;fileame5 =数据文件(“boston_subimage.tif”）;geotiffwrite (filename5 A5crop R5crop,.．.“GeoKeyDirectoryTag”info5.GeoTIFFTags.GeoKeyDirectoryTag)

显示包含裁剪图像的GeoTIFF文件。

图mapshow(filename5)芬威公园-来自GeoTIFF文件的裁剪图像)包含(“MA大陆国家飞机向东，测量英尺”) ylabel (“MA大陆国家飞机向北，测量英尺”）

例6:将高程数据写入GeoTIFF文件

将科罗拉多州南博尔德峰附近地区的海拔数据写入GeoTIFF文件。

elevFilename =“n39_w106_3arc_v2.dt1”；

从文件中读取DEM。绘制数据使用geoshow，则数据必须为类型单或双．属性指定栅格的数据类型“OutputType”名称-值对。

[Z6,R6] = readgeoraster(elevFilename，“OutputType”，“双”）;

创建一个结构来保存GeoKeyDirectoryTag信息。

键= struct(.．.“GTModelTypeGeoKey”[],.．.“GTRasterTypeGeoKey”[],.．.“GeographicTypeGeoKey”[]);

属性指示数据位于地理坐标系统中GTModelTypeGeoKey字段为2。属性指示引用对象使用贴子(而不是单元格)GTRasterTypeGeoKey字段为2。属性指示数据引用到地理坐标参考系统GeographicTypeGeoKey字段为4326。

关键。GTModelTypeGeoKey = 2;关键。GTRasterTypeGeoKey = 2;关键。geoictypegeokey = 4326;

将标高数据写入GeoTIFF文件。

Filename6 =数据文件(“southboulder.tif”）;geotiffwrite (filename6 Z6 R6,“GeoKeyDirectoryTag”键),

通过显示数据来验证数据是否已写入文件。首先，导入表示科罗拉多州边界的矢量数据readgeotable．然后，显示边界和GeoTIFF文件。

GT =重新定位表(“usastatelo.shp”）;row = GT.Name ==“科罗拉多”；科罗拉多= GT(行，:);图usamap“科罗拉多”持有在geoshow(科罗拉多州,“FaceColor”，“没有”) g = geoshow(文件名6，“DisplayType”，“网”）;demcmap (g.ZData)标题(“南博尔德峰高地”）

例7:将非图像数据写入TIFF文件

如果您使用的数据网格是类double，其值范围超出浮点强度图像所要求的限制(0 <= data <= 1)，并且如果您使用imwrite，那么你的数据将被截断到间隔[0,1]，缩放，并转换为uint8。显然，原始数据中的部分甚至全部信息都有可能丢失。要避免这些问题，并保留数据网格的数值类和范围，请使用geotiffwrite写入数据。

创建示例Z数据。

N = 512;Z7 =峰数(n);

创建一个引用对象，将行和列引用到X和Y。

R7 = maprasterref(“RasterSize”(n (n),“ColumnsStartFrom”，“北”）;R7。XWorldLimits = R7.XIntrinsicLimits; R7.YWorldLimits = R7.YIntrinsicLimits;

创建一个结构来保存GeoKeyDirectoryTag信息。将模型类型设置为1，表示投影坐标系(PCS)。

关键。GTModelTypeGeoKey = 1;

设置栅格类型为1，表示PixelIsArea(单元格)。

关键。GTRasterTypeGeoKey = 1;

使用值32767指示用户定义的投影坐标系。

关键。projtedcstypegeokey = 32767;

将数据写入GeoTIFF文件geotiffwrite．为了进行比较，使用imwrite．

fileame_geotiff =数据文件(“zdata_geotiff.tif”）;Filename_tiff =数据文件(“zdata_tiff.tif”）;R7 geotiffwrite (filename_geotiff Z7,“GeoKeyDirectoryTag”(Z7, fileame_tiff);

读取文件时使用imread数据值和类类型被保留。您可以看到TIFF文件中的数据值没有被保留。

geoZ = imread(fileame_geotiff);tiffZ = imread(filename_tiff);流(Z的类类型:%s\n，类(Z7)) fprintf(GeoTIFF文件中数据的类类型:%s\n，类(geoZ))TIFF文件中数据的类类型:%s\n，类(tiffZ))“GeoTIFF文件中的数据是否等于Z: %d\n”， isequal(geoZ, Z7))TIFF文件中的数据是否等于Z: %d\n， isequal(tiffZ, Z7))

Z类类型:double GeoTIFF文件中的数据类类型:double TIFF文件中的数据类类型:uint8 GeoTIFF文件中的数据是否等于Z: 1 TIFF文件中的数据是否等于Z: 0

您可以使用查看数据网格mapshow．

图mapshow (filename_geotiff,“DisplayType”，“texturemap”)标题(“峰值-存储在GeoTIFF文件”）

例8:修改现有文件，同时保留元信息

您可能希望修改现有文件，但保留TIFF标记中的大部分元信息(如果不是全部的话)。此示例将RGB图像从boston.tif文件写入索引图像，并将新数据写入索引GeoTIFF文件。TIFF元信息，除了BitDepth、BitsPerSample和PhotometricInterpretation标签的值，将被保留。

中读取图像boston.tifGeoTIFF文件。

[A8,R8] = readgeoraster(“boston.tif”）;

使用MATLAB函数，rgb2ind，将RGB图像转换为索引图像X使用最小方差量化。

[X8,cmap] = rgb2ind(A8,65536);

获取TIFF标签信息imfinfo．

Info8 = imfinfo(“boston.tif”）;

创建TIFF标记结构以保存从信息结构。

标签= struct(.．.“压缩”, info8。压缩,.．.“RowsPerStrip”, info8。RowsPerStrip,.．.“XResolution”, info8。XResolution,.．.“YResolution”, info8。YResolution,.．.“ImageDescription”, info8。ImageDescription,.．.“DateTime”, info8。DateTime,.．.“版权”, info8。版权,.．.“定位”, info8.Orientation);

PlanarConfiguration和ResolutionUnit标记的值必须是数值，而不是字符串值，由返回imfinfo．您可以使用Tiff类中的常量属性来设置这些标记。例如，以下是PlanarConfiguration常量属性的可能值:

Tiff。PlanarConfiguration

ans = struct with fields: Chunky: 1 Separate: 2

方法中的字符串值信息结构来获取所需的值。

标签。PlanarConfiguration =.．.Tiff.PlanarConfiguration。(info8.PlanarConfiguration);

以同样的方式设置ResolutionUnit的值。

标签。ResolutionUnit = Tiff.ResolutionUnit.(info8. resoltionunit);

的“软件”标签未设置boston.tif文件。然而,geotiffwrite将设置软件标签。为了保存信息，将值设置为空字符串，以防止标记被写入文件。

标签。软件=＇＇；

将GeoTIFF信息从boston.tif．

Geoinfo = geotiffinfo(“boston.tif”）;key = geoinfo.GeoTIFFTags.GeoKeyDirectoryTag;

写入GeoTIFF文件。

Filename8 =数据文件(“boston_indexed.tif”）;R8 geotiffwrite (filename8出数,提出,“GeoKeyDirectoryTag”的关键,“TiffTags”、标签)

查看索引图像。

图mapshow(filename8)“波士顿-索引图像”)包含(“MA大陆国家飞机向东，测量英尺”) ylabel (“MA大陆国家飞机向北，测量英尺”）

通过打印一个值表来比较结构中应该相等的信息。

Info_rgb = imfinfo(“boston.tif”）;Info_indexed = imfinfo(fileame8);tagNames =字段名(标签);tagname (strcmpi (“软件”， tagNames)) = [];名称= [{“高度”“宽度”}, tagname];间距= 2;fieldnameLength = max(cellfun(@length, names)) + spacing;formatSpec = [% - - - - - -”int2str (fieldnameLength)“年代”];fprintf([formatSpec formatSpec formatSpec' \ n '),.．.的字段名，“RGB信息”，“索引信息”) fprintf([formatSpec formatSpec formatSpec' \ n '),.．.'---------'，'---------------'，'-------------------'）为k = 1:length(names) fprintf([formatSpec formatSpec formatSpec' \ n '),.．.名字{k},.．.num2str (info_rgb。(名字{k})),.．.num2str (info_indexed。(名字{k})))结束

Fieldname RGB Information Indexed Information --------- --------------- ------------------- Height 2881 2881 Width 4481 4481 Compression Uncompressed Uncompressed RowsPerStrip 256 256 XResolution 300 300 YResolution 300 300 ImageDescription "GeoEye" "GeoEye" DateTime 2007:02:23 21:46:13 2007:02:23 21:46:13 Copyright "(c) GeoEye" "(c) GeoEye" Orientation 1 1 PlanarConfiguration Chunky Chunky resoltionunit Inch Inch

比较应该不同的信息，因为您是通过打印一个值表将RGB图像转换为索引图像。

名称= {“文件大小”，“ColorType”，“BitDepth”，.．.“BitsPerSample”，“PhotometricInterpretation”};fieldnameLength = max(cellfun(@length, names)) + spacing;formatSpec = [% - - - - - -”int2str (fieldnameLength)“年代”];formatSpec2 =“% -17年代”；流([' \ n 'formatSpec2 formatSpec2' \ n '),.．.的字段名，“RGB信息”，“索引信息”) fprintf([formatSpec formatSpec2 formatSpec2' \ n '),.．.'---------'，'---------------'，'-------------------'）为k = 1:length(names) fprintf([formatSpec formatSpec2 formatSpec2' \ n '),.．.名字{k},.．.num2str (info_rgb。(名字{k})),.．.num2str (info_indexed。(名字{k})))结束

Fieldname RGB Information Indexed Information --------- --------------- ------------------- FileSize 38729900 27925078 ColorType truecololindexed BitDepth 24 16 BitsPerSample 8 8 8 8 16 PhotometricInterpretation RGB RGB Palette

清除:删除数据文件夹

删除临时文件夹和数据文件。

删除目录(datadir“年代”）

数据集信息

的文件concord_orthow_w.tif而且concord_ortho_e.tif使用马萨诸塞州联邦地理信息局(MassGIS)的正射电像片，技术和安全服务执行办公室。有关数据集的更多信息，请使用该命令类型concord_ortho.txt．有关MassGIS提供的数据的更新链接，请参见https://www.mass.gov/info-details/massgis-data-layers．

DTED文件n39_w106_3arc_v2.dt1图片由美国地质调查局提供。

另请参阅

geotiffwrite|worldfileread|getworldfilename|geotiffinfo