导出图像和栅格到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_over .jpg)以JPEG格式存储，在“世界文件”(boston_over .jgw)中包含图像文件外部的引用信息。该图像提供了马萨诸塞州波士顿及其周边地区的低分辨率“概览”。

从JPEG文件中读取图像。

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

从世界文件中获取一个引用对象。

worldfile = getworldfilename (imagefile);R1 = worldfileread (worldfile,“地理”、尺寸(A1));

将图像写入GeoTIFF文件。

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

返回关于文件的信息RasterInfo对象。的值CoordinateReferenceSystem是一个地理坐标参考系统对象。

info1 = georasterinfo (filename1);info1。CoordinateReferenceSystem

名称:"WGS 84"基准:"World Geodetic System 1984"球体:[1×1 referenceEllipsoid] PrimeMeridian: 0 AngleUnit: "degree"

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

图usamap (R1.LatitudeLimits R1.LongitudeLimits) setm (gca),“PLabelLocation”, 0.05,“PlabelRound”2,“PlineLocation”0.05) geoshow (filename1)标题(“波士顿概述”）

示例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(filename);R3。ColumnsStartFrom

Ans = 'south' Ans = 'north'

因此，输入数据和GeoTIFF文件中的数据被翻转。

isequal (Z2 flipud (Z3))

符合逻辑的1

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

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

符合逻辑的1

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

t = Tiff (filename2' r + '）;pixelScale = getTag (t)“ModelPixelScaleTag”）;pixelScale (2) = -pixelScale (2);setTag (t)“ModelPixelScaleTag”, pixelScale);tiepoint = getTag (t)“ModelTiepointTag”）;tiepoint (5) = intrinsicToGeographic (R2, 0.5, 0.5);setTag (t)“ModelTiepointTag”, tiepoint);setTag (t)“压缩”Tiff.Compression.None)写(t, Z2);rewriteDirectory (t)关闭(t)

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

t = Tiff (filename2);Atiff =阅读(t);close(t) Rtiff = georefcells(R2.LatitudeLimits, r2 . longitude elimits,size(Atiff));Atiff isequal (Z2) isequal (R2, Rtiff)

Ans =合乎逻辑的1

例4:写入引用到投影坐标系的图像

将Concord正射影像写入单个GeoTIFF文件。这两张相邻的(从西向东)地理参考灰度(全色)正射影片覆盖了美国马萨诸塞州康科德的部分地区。txt文件表明数据引用了马萨诸塞州大陆(NAD83)州平面投影坐标系。单位是米。这对应于GeoTIFF代码编号26986，如GeoTIFF规范中所述http://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 (filename4);info4。CoordinateReferenceSystem

名称:"NAD83 / Massachusetts Mainland"地理crs: [1×1 geocrs] ProjectionMethod: "Lambert Conic Conformal (2SP)"LengthUnit: "meter" ProjectionParameters: [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”）;filename5 =丢失(“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信息。

关键=结构(.．.“GTModelTypeGeoKey”[],.．.“GTRasterTypeGeoKey”[],.．.“GeographicTypeGeoKey”[]);

通过指定地理坐标系统来指示数据GTModelTypeGeoKey字段2。属性来指示引用对象使用发帖(而不是单元格)GTRasterTypeGeoKey字段2。通过指定地理坐标参考系统来指示数据被引用GeographicTypeGeoKey在4326年。

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

将海拔数据写入GeoTIFF文件。

filename6 =丢失(“southboulder.tif”）;geotiffwrite (filename6 Z6 R6,“GeoKeyDirectoryTag”键),

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

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

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

如果您使用的数据网格是类双精度的，其值的范围超出浮点强度图像(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来指示用户定义的投影坐标系。

关键。ProjectedCSTypeGeoKey = 32767;

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

filename_geotiff =丢失(“zdata_geotiff.tif”）;filename_tiff =丢失(“zdata_tiff.tif”）;R7 geotiffwrite (filename_geotiff Z7,“GeoKeyDirectoryTag”键)imwrite (Z7 filename_tiff);

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

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

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

您可以使用mapshow．

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

例8:在保留元信息的情况下修改现有文件

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

阅读图片从boston.tifGeoTIFF文件。

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

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

[出数,提出]= rgb2ind (A8, 65536);

获取TIFF标签信息imfinfo．

info8 = imfinfo (“boston.tif”）;

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

标签=结构(.．.“压缩”, 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.ResolutionUnit);

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

标签。软件=”；

复制GeoTIFF信息boston.tif．

geoinfo = geotiffinfo (“boston.tif”）;关键= 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 (filename8);tagname =字段名(标签);tagname (strcmpi (“软件”， tagNames)) = [];名称= [{“高度”“宽度”}, tagname];间隔= 2;fieldnameLength = max(cellfun(@length, names)) + spacing;formatSpec = [% - - - - - -”int2str (fieldnameLength)“年代”];流([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})))结束

域名RGB信息索引信息--------- --------------- -------------------高度2881 2881宽度4481 4481压缩未压缩未压缩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版权"(c) GeoEye" "(c) GeoEye"定向1 1 PlanarConfiguration Chunky Chunky ResolutionUnit Inch Inch

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

名称= {“文件大小”，“ColorType”，“BitDepth”，.．.“BitsPerSample”，“PhotometricInterpretation”};fieldnameLength = max(cellfun(@length, names)) + spacing;formatSpec = [% - - - - - -”int2str (fieldnameLength)“年代”];formatSpec2 =“% -17年代”；流([' \ n 'formatSpec 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 truecolindexed BitDepth 24 16 BitsPerSample 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