export images and raster grids to geotiff -凯发k8网页登录
this example shows how to write data referenced to standard geographic and projected coordinate systems to geotiff files, using geotiffwrite
. the tagged-image file format (tiff) has emerged as a popular format to store raster data. the geotiff specification defines a set of tiff tags that describe "cartographic" information associated with the tiff raster data. using these tags, geolocated imagery or raster grids with coordinates referenced to a geographic coordinate system (latitude and longitude) or a (planar) projected coordinate system can be stored in a geotiff file.
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 under pcs_nad83_massachusetts.
read the two orthophotos.
[x_west,r_west] = readgeoraster('concord_ortho_w.tif'); [x_east,r_east] = readgeoraster('concord_ortho_e.tif');
combine the images and reference objects.
x4 = [x_west x_east]; r4 = r_west; r4.xworldlimits = [r_west.xworldlimits(1) r_east.xworldlimits(2)]; r4.rastersize = size(x4);
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 .
the dted file n39_w106_3arc_v2.dt1
is courtesy of the us geological survey.
see also
| | |