Export Images and Raster Grids to GeoTIFF

Setup: Define a Data Folder and File Name Utility Function

This example creates several temporary GeoTIFF files and uses the variable datadir to denote their location. The value used here is determined by the output of the tempdir command, but you could easily customize this. The contents of datadir are deleted at the end of the example.

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

Define an anonymous function to prepend datadir to the input file name:

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

Example 1: Write an Image Referenced to Geographic Coordinates

Write an image referenced to WGS84 geographic coordinates to a GeoTIFF file. The original image (boston_ovr.jpg) is stored in JPEG format, with referencing information external to the image file, in the "world file" (boston_ovr.jgw). The image provides a low resolution "overview" of Boston, Massachusetts, and the surrounding area.

Read the image from the JPEG file.

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

Obtain a referencing object from the world file.

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

Write the image to a GeoTIFF file.

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

Return information about the file as a RasterInfo object. Note that the value of CoordinateReferenceSystem 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"

Re-import the new GeoTIFF file and display the Boston overview image, correctly located, on a map.

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

Example 2: Write a Data Grid Referenced to Geographic Coordinates

Load elevation raster data and a geographic cells reference object. Write the data grid to a GeoTIFF file.

load topo60c
Z2 = 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 world
gridm off
setm(gca,'MLabelParallel',-90,'MLabelLocation',90)
tmap = geoshow(filename2,'DisplayType','texturemap');
demcmap(tmap.CData)
title('Elevation Data Grid')

Example 3: Change Data Organization of GeoTIFF Files

When you write data using geotiffwrite or read data using readgeoraster, the columns of the data grid start from north and the rows start from west. For example, the input data from topo60c.mat starts from south, but the output data from topo60c.tif starts 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.

isequal(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);
close(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 at http://geotiff.maptools.org/spec/geotiff6.html#6.3.3.1 under PCS_NAD83_Massachusetts.

Read the two orthophotos. Create reference objects for the orthophotos by reading their world files.

X_west = imread('concord_ortho_w.tif');
X_east = imread('concord_ortho_e.tif');

R_west = worldfileread("concord_ortho_w.tfw","planar",size(X_west));
R_east = worldfileread("concord_ortho_e.tfw","planar",size(X_east));

Merge the orthophotos.

[X4,R4] = mergetiles(X_west,R_west,X_east,R_east);

Write the data to a GeoTIFF file. Use the code number, 26986, indicating the PCS_NAD83_Massachusetts Projected Coordinate System.

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

Return information about the file as a RasterInfo object. Note that the value of CoordinateReferenceSystem 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.

figure
mapshow(filename4)
title('Combined Orthophotos')
xlabel('MA Mainland State Plane easting, meters')
ylabel('MA Mainland State Plane northing, meters')

Example 5: Write a Cropped Image from a GeoTIFF File

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

Read the RGB image and referencing information from the boston.tif GeoTIFF file.

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

Crop the image.

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 = datafile('boston_subimage.tif');
geotiffwrite(filename5,A5crop,R5crop, ...
   'GeoKeyDirectoryTag',info5.GeoTIFFTags.GeoKeyDirectoryTag)

Display the GeoTIFF file containing the cropped image.

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')

Example 6: Write Elevation Data to GeoTIFF File

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

elevFilename = 'n39_w106_3arc_v2.dt1';

Read the DEM from the file. To plot the data using geoshow, the data must be of type single or double. Specify the data type for the raster using the 'OutputType' name-value pair.

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

Create a structure to hold the GeoKeyDirectoryTag information.

key = struct( ...
    'GTModelTypeGeoKey',[], ...
    'GTRasterTypeGeoKey',[], ...
    'GeographicTypeGeoKey',[]);

Indicate the data is in a geographic coordinate system by specifying the GTModelTypeGeoKey field as 2. Indicate that the reference object uses postings (rather than cells) by specifying the GTRasterTypeGeoKey field as 2. Indicate the data is referenced to a geographic coordinate reference system by specifying the GeographicTypeGeoKey field 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',key)

Verify the data has been written to a file by displaying it. First, import vector data that represents the state boundary of Colorado using readgeotable. Then, display the boundary and GeoTIFF file.

GT = readgeotable('usastatelo.shp');
row = GT.Name == 'Colorado';
colorado = GT(row,:);

figure
usamap 'Colorado'
hold on
geoshow(colorado,'FaceColor','none')
g = geoshow(filename6,'DisplayType','mesh');
demcmap(g.ZData)
title('South Boulder Peak Elevation')

Example 7: Write Non-Image Data to a TIFF File

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 using imwrite, 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, use geotiffwrite to write the data.

Create sample Z data.

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;

Create a structure to hold the GeoKeyDirectoryTag information. Set the model type to 1 indicating Projected Coordinate System (PCS).

key.GTModelTypeGeoKey = 1;

Set the raster type to 1 indicating PixelIsArea (cells).

key.GTRasterTypeGeoKey = 1;

Indicate a user-defined Projected Coordinate System by using a value of 32767.

key.ProjectedCSTypeGeoKey = 32767;

Write the data to a GeoTIFF file with geotiffwrite. For comparison, write a second file using imwrite.

filename_geotiff = datafile('zdata_geotiff.tif');
filename_tiff = datafile('zdata_tiff.tif');
geotiffwrite(filename_geotiff,Z7,R7,'GeoKeyDirectoryTag',key)
imwrite(Z7, filename_tiff);

When you read the file using imread the 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('Class type of 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('Does data in    TIFF file equal 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 using mapshow.

figure
mapshow(filename_geotiff,'DisplayType','texturemap')
title('Peaks - Stored in GeoTIFF File')

Example 8: Modify an Existing File While Preserving Meta Information

You may want to modify an existing file, but preserve most, if not all, of the meta information in the TIFF tags. This example converts the RGB image from the boston.tif file 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.

Read the image from the boston.tif GeoTIFF file.

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

Use the MATLAB function, rgb2ind, to convert the RGB image to an indexed image X using minimum variance quantization.

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

Obtain the TIFF tag information using imfinfo.

info8 = imfinfo('boston.tif');

Create a TIFF tags structure to preserve selected information from the info structure.

tags = struct( ...
    'Compression',  info8.Compression, ...
    'RowsPerStrip', info8.RowsPerStrip, ...
    'XResolution',  info8.XResolution, ...
    'YResolution',  info8.YResolution, ...
    'ImageDescription', info8.ImageDescription, ...
    'DateTime',    info8.DateTime, ...
    'Copyright',   info8.Copyright, ...
    'Orientation', info8.Orientation);

The values for the PlanarConfiguration and ResolutionUnit tags must be numeric rather than string valued, as returned by imfinfo. You can set these tags by using the constant properties from the Tiff class. For example, here are the possible values for the PlanarConfiguration constant property:

Tiff.PlanarConfiguration

ans = 

  struct with fields:

      Chunky: 1
    Separate: 2

Use the string value from the info structure to obtain the desired value.

tags.PlanarConfiguration = ...
    Tiff.PlanarConfiguration.(info8.PlanarConfiguration);

Set the ResolutionUnit value in the same manner.

tags.ResolutionUnit = Tiff.ResolutionUnit.(info8.ResolutionUnit);

The Software tag is not set in the boston.tif file. However, geotiffwrite will set the Software tag by default. To preserve the information, set the value to the empty string which prevents the tag from being written to the file.

tags.Software = '';

Copy the GeoTIFF information from boston.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',key,'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'];

spacing = 2;
fieldnameLength = max(cellfun(@length, names)) + spacing;
formatSpec = ['%-' int2str(fieldnameLength) 's'];

fprintf([formatSpec formatSpec formatSpec '\n'], ...
    'Fieldname', 'RGB Information', 'Indexed Information')
fprintf([formatSpec formatSpec formatSpec '\n'], ...
    '---------', '---------------', '-------------------')

for k = 1:length(names)
    fprintf([formatSpec formatSpec formatSpec '\n'], ...
        names{k}, ...
        num2str(info_rgb.(names{k})), ...
        num2str(info_indexed.(names{k})))
end

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               
ResolutionUnit       Inch                 Inch                 

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' formatSpec formatSpec2 formatSpec2 '\n'], ...
    'Fieldname', 'RGB Information', 'Indexed Information')
fprintf([formatSpec formatSpec2 formatSpec2 '\n'], ...
    '---------', '---------------', '-------------------')
for k = 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

Remove the temporary folder and data files.

rmdir(datadir, 's')

Data Set Information

The files boston.tif and boston_ovr.jpg include materials copyrighted by GeoEye, all rights reserved. GeoEye was merged into the DigitalGlobe corporation on January 29th, 2013. For more information about the data sets, use the commands type boston.txt and type boston_ovr.txt.

The files concord_orthow_w.tif and concord_ortho_e.tif are derived using orthophoto tiles from the Bureau of Geographic Information (MassGIS), Commonwealth of Massachusetts, Executive Office of Technology and Security Services. For more information about the data sets, use the command type concord_ortho.txt. For an updated link to the data provided by MassGIS, see https://www.mass.gov/info-details/massgis-data-layers.

The DTED file n39_w106_3arc_v2.dt1 is courtesy of the US Geological Survey.