Setup: Define a Data Folder and File Name Utility Function

This example creates several temporary GeoTIFF files and uses the variabledatadirto denote their location. The value used here is determined by the output of thetempdircommand, but you could easily customize this. The contents ofdatadirare deleted at the end of the example.

datadir = fullfile(tempdir,'datadir');if~exist(datadir,'dir') mkdir(datadir)end

Define an anonymous function to prependdatadir输入文件名：

dataFile = @（filename）fullfile（datadir，filename）;

Example 1: Write an Image Referenced to Geographic Coordinates

编写引用WGS84地理坐标的图像到Geotiff文件。原始图像（boston_ovr.jpg）以JPEG格式存储，并在“ World File”（Boston_ovr.jgw）中使用图像文件外部的引用信息。该图像提供了波士顿，马萨诸塞州和周边地区的低分辨率“概述”。

从JPEG文件中读取图像。

Basename ='boston_ovr'; imagefile = [basename'.jpg']; A1 = imread(imagefile);

Obtain a referencing object from the world file.

世界file = getworldfilename(imagefile); R1 = worldfileread(worldfile,'geographic',size(A1));

将图像写入Geotiff文件。

filename1 = datafile([basename'.tif']); geotiffwrite(filename1,A1,R1)

再保险turn information about the file as arasterinfoobject. Note that the value of坐标系is a geographic coordinate reference system object.

info1 = georasterinfo(filename1); info1.CoordinateReferenceSystem

ans = geocrs with properties: Name: "WGS 84" Datum: "World Geodetic System 1984" Spheroid: [1×1 referenceEllipsoid] PrimeMeridian: 0 AngleUnit: "degree"

再保险-import the new GeoTIFF file and display the Boston overview image, correctly located, in a map axes.

figure usamap(R1.LatitudeLimits,R1.LongitudeLimits) setm(gca,“识别”,0.05,'PlabelRound',-2,'PlineLocation',0.05) geoshow(filename1) title('Boston Overview')

Example 2: Write a Data Grid Referenced to Geographic Coordinates

加载高程栅格数据和地理单元参考对象。将数据网格写入Geotiff文件。

loadtopo60cZ2 = topo60c; R2 = topo60cR; filename2 = datafile('topo60c.tif');geotiffwrite(filename2,Z2,R2)

The values in the data grid range from -7473 to 5731. Display the grid as a texture-mapped surface rather than as an intensity image.

figure worldmap世界gridm离开setm(gca,“ mlabelparallel”,-90,'MLabelLocation',90) tmap = geoshow(filename2,'显示类型','texturemap');demcmap（tmap.cdata）标题（'Elevation Data Grid')

Example 3: Change Data Organization of GeoTIFF Files

When you write data usinggeotiffwrite或者read data usingReadgeoraster, the columns of the data grid start from north and the rows start from west. For example, the input data fromtopo60c.mat从南方开始，但来自topo60c.tifstarts from north.

R2.ColumnsStartFrom [Z3,R3] = readgeoraster(filename2); R3.ColumnsStartFrom

ans = 'south' ans = 'north'

Therefore, the input data and data in the GeoTIFF file is flipped.

quequal（Z2，Flipud（Z3））

ans = logical 1

If you need the data in your workspace to match again, then flip the Z values and set the referencing object such that the columns start from the south:

R3.ColumnsStartFrom ='south'; Z3 = flipud(Z3); isequal(Z2,Z3)

ans = logical 1

The data in the GeoTIFF file is encoded with positive scale values. Therefore, when you view the file with ordinary TIFF-viewing software, the northern edge of the data set is at the top. To make the data layout in the output file match the data layout of the input, you can construct a Tiff object and use it to reset some of the tags and the image data.

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,'Compression', Tiff.Compression.None) write(t,Z2); rewriteDirectory(t) close(t)

Verify the referencing object and data grid from the input data match the output data file. To do this, read the Tiff image and create a reference object. Then, compare the grids.

t = tiff（filename2）;atiff = read（t）;关闭（t）rtiff = georefcells（r2.latitudelimits，r2.longitudelimits，size（atiff））;iSequal（Z2，ATIFF）iSequal（R2，RTIFF）

ans = logical 1 ans = logical 1

Example 4: Write an Image Referenced to a Projected Coordinate System

Write the Concord orthophotos to a single GeoTIFF file. The two adjacent (west-to-east) georeferenced grayscale (panchromatic) orthophotos cover part of Concord, Massachusetts, USA. The concord_ortho.txt file indicates that the data are referenced to the Massachusetts Mainland (NAD83) State Plane Projected Coordinate System. Units are meters. This corresponds to the GeoTIFF code number 26986 as noted in the GeoTIFF specification athttp://geotiff.maptools.org/spec/geotiff6.html#6.3.3.1under PCS_NAD83_Massachusetts.

再保险ad the two orthophotos.

[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文件写数据。使用num的代码ber, 26986, indicating the PCS_NAD83_Massachusetts Projected Coordinate System.

coordRefSysCode = 26986; filename4 = datafile('concord_ortho.tif');geotiffwrite(filename4,X4,R4,'CoordRefSysCode',coordRefSysCode)

再保险turn information about the file as arasterinfoobject. Note that the value of坐标系is a projected coordinate reference system object.

info4 = georasterinfo(filename4); info4.CoordinateReferenceSystem

ans = projcrs with properties: Name: "NAD83 / Massachusetts Mainland" GeographicCRS: [1×1 geocrs] ProjectionMethod: "Lambert Conic Conformal (2SP)" LengthUnit: "meter" ProjectionParameters: [1×1 map.crs.ProjectionParameters]

Display the combined Concord orthophotos.

图mapshow（filename4）标题（“矫形器的组合”）xlabel（“马兰大陆国家升级，仪表”）ylabel（'MA大陆国家飞机北方，米)

Example 5: Write a Cropped Image from a GeoTIFF File

Write a subset of a GeoTIFF file to a new GeoTIFF file.

再保险ad the RGB image and referencing information from theboston.tifGeoTIFF file.

[A5,R5] = readgeoraster('boston.tif');

裁剪图像。

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

Write the cropped image to a GeoTIFF file. Use the GeoKeyDirectoryTag from the original GeoTIFF file.

info5 = geotiffinfo('boston.tif');filename5 =丢失('boston_subimage.tif');geotiffwrite(filename5,A5crop,R5crop,。。。'GeoKeyDirectoryTag'，info5.geotifftags.geokeydirectorytag）

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

figure mapshow(filename5) title('Fenway Park - Cropped Image from GeoTIFF File'）xlabel（'MA Mainland State Plane easting, survey feet'）ylabel（'MA Mainland State Plane northing, survey feet')

示例6：将高程数据写入Geotiff文件

Write elevation data for an area around South Boulder Peak in Colorado to a GeoTIFF file.

elevFilename ='n39_w106_3arc_v2.dt1';

再保险ad the DEM from the file. To plot the data usinggeoshow，数据必须是类型single或者double。Specify the data type for the raster using the'OutputType'name-value pair.

[Z6,R6] = readgeoraster(elevFilename,'OutputType','double');

创建一个结构以保存GeokeyDirectoryTag信息。

key = struct(。。。'gtmodeltypegeokey'，[]，，。。。“GTRasterTypeGeoKey”，[]，，。。。'GeographicTypeGeoKey',[]);

通过指定数据，指示数据在地理坐标系中GTModelTypeGeoKeyfield as 2. Indicate that the reference object uses postings (rather than cells) by specifying theGTRasterTypeGeoKeyfield as 2. Indicate the data is referenced to a geographic coordinate reference system by specifying theGeographicTypeGeoKeyfield as 4326.

key.GTModelTypeGeoKey = 2; key.GTRasterTypeGeoKey = 2; key.GeographicTypeGeoKey = 4326;

Write the elevation data to a GeoTIFF file.

filename6 = datafile('southboulder.tif');geotiffwrite(filename6,Z6,R6,'GeoKeyDirectoryTag'，钥匙）

通过显示数据，验证数据已写入文件。首先，使用使用科罗拉多州状态边界的导入向量数据readgeotable。Then, display the boundary and GeoTIFF file.

GT = readgeotable('usastatelo.shp');行= gt.name =='Colorado'; colorado = GT(row,:); figure usamap'Colorado'hold上geoshow(colorado,'FaceColor','none') g = geoshow(filename6,'显示类型','mesh');demcmap（g.zdata）标题（'South Boulder Peak Elevation')

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

If you are working with a data grid that is class double with values that range outside the limits required of a floating point intensity image (0 <= data <= 1), and if you store the data in a TIFF file usingimwrite, then your data will be truncated to the interval [0,1], scaled, and converted to uint8. Obviously it is possible for some or even all of the information in the original data to be lost. To avoid these problems, and preserve the numeric class and range of your data grid, usegeotiffwriteto write the data.

创建示例Z数据。

n = 512; Z7 = peaks(n);

Create a referencing object to reference the rows and columns to X and Y.

R7 = maprasterref('RasterSize'，[n n]，'ColumnsStartFrom','north');R7.XWorldLimits = R7.XIntrinsicLimits; R7.YWorldLimits = R7.YIntrinsicLimits;

创建一个结构以保存GeokeyDirectoryTag信息。Set the model type to 1 indicating Projected Coordinate System (PCS).

key.GTModelTypeGeoKey = 1;

将光栅类型设置为1指示Pixelisarea（单元格）的1。

key.GTRasterTypeGeoKey = 1;

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

key.ProjectedCSTypeGeoKey = 32767;

Write the data to a GeoTIFF file withgeotiffwrite。For comparison, write a second file usingimwrite。

filename_geotiff = datafile('zdata_geotiff.tif');filename_tiff = datafile('zdata_tiff.tif');GeotiffWrite（filename_geotiff，z7，r7，'GeoKeyDirectoryTag'，钥匙）imwrite(Z7, filename_tiff);

When you read the file usingimreadthe data values and class type are preserved. You can see that the data values in the TIFF file are not preserved.

geoZ = imread(filename_geotiff); tiffZ = imread(filename_tiff); fprintf('类型z：％s \ n', class(Z7)) fprintf('Class type of data in GeoTIFF file: %s\n', class(geoZ)) fprintf('Class type of data in TIFF file: %s\n', class(tiffZ)) fprintf('Does data in GeoTIFF file equal Z: %d\n', isequal(geoZ, Z7)) fprintf('在tiff文件中做数据等于z：％d \ n', isequal(tiffZ, Z7))

Class type of Z: double Class type of data in GeoTIFF file: double Class type of data in TIFF file: uint8 Does data in GeoTIFF file equal Z: 1 Does data in TIFF file equal Z: 0

You can view the data grid usingmapshow。

figure mapshow(filename_geotiff,'显示类型','texturemap') title('Peaks - Stored in GeoTIFF File')

Example 8: Modify an Existing File While Preserving Meta Information

您可能需要修改现有文件，但保留TIFF标签中的元信息的大部分（如果不是全部）。此示例从boston.tiffile into an indexed image and writes the new data to an indexed GeoTIFF file. The TIFF meta-information, with the exception of the values of the BitDepth, BitsPerSample, and PhotometricInterpretation tags, is preserved.

再保险ad the image from theboston.tifGeoTIFF file.

[A8,R8] = readgeoraster('boston.tif');

Use the MATLAB function,rgb2ind，将RGB图像转换为索引图像Xusing minimum variance quantization.

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

使用imfinfo。

info8 = imfinfo('boston.tif');

创建一个TIFF标签结构，以保留从info结构。

标签= struct（。。。'Compression'，info8.compression，。。。'RowsPerStrip', info8.RowsPerStrip,。。。'XResolution', info8.XResolution,。。。'YResolution', info8.YResolution,。。。'ImageDescription', info8.ImageDescription,。。。'DateTime', info8.DateTime,。。。'Copyright', info8.Copyright,。。。'方向', info8.Orientation);

PlanarConfiguration和andolutionUnit标签的值必须是数字而不是符号的字符串，如由imfinfo。您可以使用TIFF类的常数属性来设置这些标签。例如，以下是PlanarConfiguration常数属性的可能值：

Tiff.PlanarConfiguration

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

使用从infostructure to obtain the desired value.

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

Set the ResolutionUnit value in the same manner.

标签。再保险solutionUnit = Tiff.ResolutionUnit.(info8.ResolutionUnit);

软件标签未在boston.tiffile. However,geotiffwritewill set theSoftwaretag by default. To preserve the information, set the value to the empty string which prevents the tag from being written to the file.

标签。Software ='';

Copy the GeoTIFF information fromboston.tif。

geoinfo = geotiffinfo('boston.tif');key = geoinfo.GeoTIFFTags.GeoKeyDirectoryTag;

Write to the GeoTIFF file.

filename8 = datafile('boston_indexed.tif');geotiffwrite(filename8,X8,cmap,R8,'GeoKeyDirectoryTag'，钥匙，'TiffTags',tags)

View the indexed image.

figure mapshow(filename8) title('Boston - Indexed Image'）xlabel（'MA Mainland State Plane easting, survey feet'）ylabel（'MA Mainland State Plane northing, survey feet')

Compare the information in the structures that should be equal by printing a table of the values.

info_rgb = imfinfo（'boston.tif');info_indexed = imfinfo(filename8); tagNames = fieldnames(tags); tagNames(strcmpi('Software', tagNames)) = []; names = [{'Height''Width'}，tagnames'];间距= 2;fieldNameLength = max（cellfun（@length，names）） +间距;格式pec = ['%-'int2str（fieldnameLength）'s']; fprintf([formatSpec formatSpec formatSpec'\ n'],。。。'Fieldname','RGB Information','Indexed Information') fprintf([formatSpec formatSpec formatSpec'\ n'],。。。'-------------','-----------------------,'-------------------')fork = 1：长度（名称）fprintf（[formatspec格式PEMATSPEC格式PEMATSPEC'\ n'],。。。names{k},。。。num2str(info_rgb.(names{k})),。。。num2str(info_indexed.(names{k})))end

Fieldname RGB Information Indexed Information --------- --------------- ------------------- Height 28812881宽度4481 4481压缩未压缩的未压缩牛皮条256 256 256 XResolution 300 300 yresolution 300 300 300 ImagedScription“ Geoeye”“ Geoeye”“ Geoeye” DateTime 2007：02：23 21:46:13 2007：02：23 21:46:46:13 21:46:13 pocyright“”（c）geoeye”方向1 1平面配置厚实的矮小的分辨率英寸英寸英寸英寸英寸

Compare the information that should be different, since you converted an RGB image to an indexed image, by printing a table of values.

names = {'FileSize','ColorType','BitDepth',。。。'BitsPerSample','PhotometricInterpretation'}; fieldnameLength = max(cellfun(@length, names)) + spacing; formatSpec = ['%-'int2str（fieldnameLength）'s']; formatSpec2 ='%-17s'; fprintf(['\ n'格式PEMATSPEC2 FORMATSPEC2'\ n'],。。。'Fieldname','RGB Information','Indexed Information'）fprintf（[格式PEMATSATSPEC2 FORMATSPEC2'\ n'],。。。'-------------','-----------------------,'-------------------')fork = 1:length(names) fprintf([formatSpec formatSpec2 formatSpec2'\ n'],。。。names{k},。。。num2str(info_rgb.(names{k})),。。。num2str(info_indexed.(names{k})))end

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

Cleanup: Remove Data Folder

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

rmdir(datadir,'s')

Data Set Information

The filesboston.tif和boston_ovr.jpg包括Geoeye版权所有的材料，保留所有权利。Geoeye于2013年1月29日合并为DigitalGlobe Corporation。有关数据集的更多信息，请使用命令type boston.txt和type boston_ovr.txt。

The filesconcord_orthow_w.tif和concord_ortho_e.tif使用来自地理信息局（Massgis），马萨诸塞州联邦，技术和安全服务执行办公室的正常光瓷砖得出。有关数据集的更多信息，请使用命令键入concord_ortho.txt。For an updated link to the data provided by MassGIS, seehttps://www.mass.gov/info-details/massgis-data-layers。

DTED文件n39_w106_3arc_v2.dt1is courtesy of the US Geological Survey.