%GEOTIFFWRITE Write a 2D or 3D array to a single or multi-band GeoTIFF file
%
% MATLAB's Mapping Toolbox only provides a "geotiffread" function, but
% it does not have a "geotiffwrite" function (Note). This is the MATLAB
% program to write a 2D or 3D array to a single or multi-band GeoTIFF
% file, where data can be either 1-bit monochrome data (i.e. binary or
% bilevel), 8-bit unsigned integer, 16-bit signed / unsigned integer,
% 32-bit signed integer, or 32-bit floating point.
%
% Note: Starting from version R2011a, MATLAB's Mapping Toolbox also
% provides its own "geotiffwrite" function
% (http://www.mathworks.com/help/toolbox/map/rn/bsq4us7-1.html#bsu1ro5-1).
%
% This program is based on GeoTIFF Format Specification under:
% http://www.remotesensing.org/geotiff/spec/geotiffhome.html
% or http://download.osgeo.org/geotiff/spec
%
% It does not need MATLAB's Mapping Toolbox, or any other library.
%
% The motivation to create such a program is because I want to convert
% the free satelite data from NASA & USGS to the GARMIN topo map for my
% Garmin GPS. DEM2TOPO (http://people.uleth.ca/~brad.gom/dem2topo) is a
% wonderful program that can nicely extract contours from the elevation
% data. However, only 1201x1201 pixels data are supported. For USGS24K
% DEM data, 1-arc-second SRTM1 data, etc, DEM2TOPO cannot process them
% directly, which must be converted into GeoTIFF format first (see
% "Known Issues" in the above site). Although there is a free Windows
% program "3DEM" that can convert USGS24K DEM & SRTM1 data into GeoTIFF
% format, sometimes you will still need to interpolate those DEM data
% for void or missing elevation points.
%
% If you just want to write DEM data to a GeoTIFF file, all you have to
% provide is a bounding box (2 latitude & 2 longitude), and a 2D or 3D
% array. Otherwise, you can provide proper GeoTIFF fields to the option
% argument. To make your life easier, a GUI program "make_option.m"
% is also provided. So you can select the proper GeoTIFF fields to
% generate the option argument.
%
% I recommend that you use "make_option" program to generate option
% argument instead of constructing it manually.
%
% Usage: [option bbox] = make_option([bbox]);
% geotiffwrite(filename, bbox, image, [bit_depth, option]);
%
% filename - GeoTIFF file name to be written.
%
% bbox - Bounding box with West/East of longitude and South/North of
% latitude. The latitude & longitude must be in Decimal Degree
% (DD) format, and they must be in a 2x2 array with the order
% of:
%
% [ longitude_West, latitude_South;
% longitude_East, latitude_North ]
%
% Note: If GTModelTypeGeoKey is not ModelTypeGeographic, bbox
% will be reset to empty.
%
% image - 2D or 3D array. Although any orientation can be specified
% through option.Orientation, it is simpler to just use the
% following default orientation:
%
% First row at North edge; last row at South edge;
% First column at West edge; last column at East edge;
%
% Note: By default, "geotiffwrite" assumes intensity image and
% option.FullColor is 0. So the third dimension of a 3D
% image represents different bands. If option.FullColor
% is 1, the third dimension of 3D image must be 3, which
% indicates R,G,B. i.e. when the third dimension of a 3D
% array equals to 3, it represents a 3-band GeoTIFF
% image, unless the option.FullColor is set to 1. When
% option.ColorMap is set, the specified palette will be
% used, and you cannot set option.FullColor to 1 at the
% same time.
%
% bit_depth - (optional) Number of bits to represent a data point (i.e.
% bits per sample):
%
% 1 for 1-bit monochrome data (i.e. binary or bilevel)
% 8 for 8-bit unsigned integer
% 16 for 16-bit signed integer
% 32 for 32-bit floating point (default)
% -16 for 16-bit unsigned integer
% -32 for 32-bit signed integer
%
% Note: If you set values for option.ColorMap, 16 will be for
% 16-bit unsigned integer, and you cannot write in 32-bit
% depth or 1-bit depth. If you set option.FullColor to 1,
% 8-bit unsigned integer will be used, and you cannot
% change bit_depth for Full Color image.
%
% option - (optional) A structure of GeoTIFF fields:
%
% First, please run: [option bbox] = make_option( [bbox] ); This
% will include most fields, except the following:
%
% option.NaN
% By default, NaN value (void data point) in the image
% variable will be kept untouched. However, if you decide
% to replace it with other value (e.g. -32768), you can
% define it here.
%
% option.FullColor = [0] | 1
% 0: Intensity or ColorMap image; (default)
% 1: FullColor (24-bit RGB) image;
%
% option.ColorMap [Nx3 array]
% Where N = 2^bit_depth, and bit_depth can only be 8 or 16.
% Values in option.ColorMap range from 0 to 65535, and value
% 0 in image points to the first row in the ColorMap. Columns
% in option.ColorMap are [R G B]. Black is represented by
% [0 0 0], and white is represented by [65535 65535 65535].
%
% option.Orientation = [1] | 2 | 3 | 4 | 5 | 6 | 7 | 8
% 1: Row from Top, Col from Left; (default)
% 2: Row from Top, Col from Right;
% 3: Row from Bottom, Col from Right;
% 4: Row from Bottom, Col from Left;
% 5: Row from Left, Col from Top;
% 6: Row from Right, Col from Top;
% 7: Row from Right, Col from Bottom;
% 8: Row from Left, Col from Bottom;
%
% option.ModelTransformationTag [4x4 array]
% In most case, ModelPixelScaleTag and ModelTiepointTag from
% "make_option.m" program is sufficient to transform from raster
% to model space. If raster image requires rotation or shearing,
% you can use this ModelTransformationTag to define a 4x4 affine
% transformation matrix.
%
% Examples:
% - Although not illustrated, it would be easier to use "make_option.m"
% for option argument
%
% Data from UTM projected coordinate systems:
% ftp://ftp.remotesensing.org/pub/geotiff/samples/spot/chicago/UTM2GTIF.TIF
% or http://download.osgeo.org/geotiff/samples/spot/chicago/UTM2GTIF.TIF
% img = imread('UTM2GTIF.TIF');
% clear option
% option.GTModelTypeGeoKey = 1;
% option.ModelPixelScaleTag = [10;10;0];
% option.ModelTiepointTag = [0;0;0;444650;4640510;0];
% option.GeogEllipsoidGeoKey = 7008;
% option.ProjectedCSTypeGeoKey = 26716;
% option.GTCitationGeoKey = 'UTM Zone 16N NAD27"';
% option.GeogCitationGeoKey = 'Clarke, 1866 by Default';
% geotiffwrite('UTM2GTIFX.TIF',[],img,8,option);
%
% Data from colormap indexed scan image:
% ftp://ftp.remotesensing.org/pub/geotiff/samples/usgs/o41078a5.tif
% or http://download.osgeo.org/geotiff/samples/usgs/o41078a5.tif
% [img cmap] = imread('o41078a5.tif');
% clear option
% option.GTModelTypeGeoKey = 1;
% option.ModelPixelScaleTag = [2.4384;2.4384;0];
% option.ModelTiepointTag = [0;0;0;698753.304798;4556059.506392;0];
% option.ProjectedCSTypeGeoKey = 32617;
% option.PCSCitationGeoKey = 'UTM Zone 17 N with WGS84';
% option.ColorMap = 65535*cmap;
% geotiffwrite('o41078a5X.tif',[],img,8,option);
%
% DEM data:
% ftp://ftp.remotesensing.org/pub/geotiff/samples/usgs/mlatlon.tif
% or http://download.osgeo.org/geotiff/samples/usgs/mlatlon.tif
% img = imread('mlatlon.tif');
% clear option
% option.GeographicTypeGeoKey = 4267;
% option.GTCitationGeoKey = 'Geographic Model/GeoTIFF 1.0';
% option.GeogCitationGeoKey = 'NAD 27 Datum';
% geotiffwrite('mlatlonX.tif',[-122.6261,37.4531;-122.0777,38],img,16,option);
%
% Even simpler:
% http://dds.cr.usgs.gov/srtm/version2_1/SRTM3/North_America/N46W082.hgt.zip
% srtm = load_srtm('N46W082.hgt');
% geotiffwrite('N46W082.tif',[-82,46;-81,47],flipud(srtm.dat));
%
% - Jimmy Shen (jimmy@rotman-baycrest.on.ca)
%
function geotiffwrite(filename, bbox, image, bit_depth, option)
% Set default inputs
%
if nargin < 5
option = [];
if isempty(bbox)
help geotiffwrite;
error('Because you did not specify bounding box, please assign proper values to option argument');
end
end
if nargin < 4
bit_depth = [];
end
if nargin < 3
error('Usage: geotiffwrite(filename, bbox, image, [bit_depth, option])');
end
if isempty(bit_depth)
bit_depth = 32;
end
if isfield(option, 'FullColor') & option.FullColor == 1
bit_depth = 8;
if isfield(option, 'ColorMap')
error('If you set option.FullColor to 1, you cannot set option.ColorMap.');
end
else
option.FullColor = 0;
end
if isfield(option, 'ColorMap') & ( abs(bit_depth)==32 | abs(bit_depth)==1 )
error('If you set values to option.ColorMap, you cannot write in 32-bit depth or 1-bit depth.');
end
% Set void data: GeoTIFF set void data to -32768, while data from
% other software can be -32767 or NaN.
%
% if strcmpi( class(image), 'double' )
% image(find(isnan(image))) = -32768;
% end
% Only for int16 DEM data
%
% if bit_depth == 16
% image(find(image==-32767)) = -32768;
% end
if isfield(option, 'NaN')
image(find(isnan(image))) = option.NaN;
end
% Set Orientation, 274
%
if isfield(option, 'Orientation')
ifd.Orientation = option.Orientation;
for i = 1:size(image,3)
switch ifd.Orientation
case 2
image(:,:,i) = fliplr(image(:,:,i));
case 3
image(:,:,i) = fliplr(flipud(image(:,:,i)));
case 4
image(:,:,i) = flipud(image(:,:,i));
case 5
image(:,:,i) = flipud(rot90(image(:,:,i)));
case 6
image(:,:,i) = rot90(image(:,:,i));
case 7
image(:,:,i) = fliplr(rot90(image(:,:,i)));
case 8
image(:,:,i) = fliplr(flipud(rot90(image(:,:,i))));
end
end
else
% First row at North edge; last row at South edge;
% First column at West edge; last column at East edge;
%
ifd.Orientation = 1;
end
if abs(bit_depth) == 1 & mod(size(image,2),8) ~= 0
error('Please use integral multiple of 8 for 1-bit image size.');
end
GeoAsciiParamsTag = '';
GeoDoubleParamsTag = [];
GeoAsciiOffset = 0;
GeoDoubleOffset = 0;
NumberOfKeys = 0;
GeoKeyDirectoryTag = [1 1 0 NumberOfKeys];
% Set GTModelTypeGeoKey, 1024
%
NumberOfKeys = NumberOfKeys + 1;
GeoKeyDirectoryTag(1, 4) = NumberOfKeys;
if isfield(option, 'GTModelTypeGeoKey') & option.GTModelTypeGeoKey ~= 2
code = option.GTModelTypeGeoKey;
bbox = [];
disp('Warning: Since your Model Type is not Geographic, bbox is reset to empty.');
else
code = 2; % ModelTypeGeographic
end
GeoKeyDirectoryTag = [GeoKeyDirectoryTag; [1024 0 1 code]];
% Set GTRasterTypeGeoKey, 1025
%
NumberOfKeys = NumberOfKeys + 1;
GeoKeyDirectoryTag(1, 4) = NumberOfKeys;
if isfield(option, 'GTRasterTypeGeoKey')
code = option.GTRasterTypeGeoKey;
else
code = 1; % RasterPixelIsArea
end
gifd.GTRasterTypeGeoKey = code;
GeoKeyDirectoryTag = [GeoKeyDirectoryTag; [1025 0 1 code]];
% Set GTCitationGeoKey, 1026
%
if isfield(option, 'GTCitationGeoKey')
NumberOfKeys = NumberOfKeys + 1;
GeoKeyDirectoryTag(1, 4) = NumberOfKeys;
code = [option.GTCitationGeoKey(:); '|']; cnt = length(code);
GeoKeyDirectoryTag = [GeoKeyDirectoryTag; [1026 34737 cnt GeoAsciiOffset]];
GeoAsciiParamsTag = [GeoAsciiParamsTag; code];
GeoAsciiOffset = GeoAsciiOffset + cnt;
end
% Set GeographicTypeGeoKey, 2048
%
if isfield(option, 'GeographicTypeGeoKey')
NumberOfKeys = NumberOfKeys + 1;
GeoKeyDirectoryTag(1, 4) = NumberOfKeys;
code = option.GeographicTypeGeoKey;
GeoKeyDirectoryTag = [GeoKeyDirectoryTag; [2048 0 1 code]];
elseif ~isempty(bbox)
NumberOfKeys = NumberOfKeys + 1;
GeoKeyDirectoryTag(1, 4) = NumberOfKeys;
% code = 4267; % GCS_NAD27
% code = 4322; % GCS_WGS_72
code = 4326; % GCS_WGS_84
% code = 4269; % GCS_NAD83
GeoKeyDirectoryTag = [GeoKeyDirectoryTag; [2048 0 1 code]];
end
% Set GeogCitationGeoKey, 2049
%
if isfield(option, 'GeogCitationGeoKey')
NumberOfKeys = NumberOfKeys + 1;
GeoKeyDirectoryTag(1, 4) = NumberOfKeys;
code = [option.GeogCitationGeoKey(:); '|']; cnt = length(code);
GeoKeyDirectoryTag = [GeoKeyDirectoryTag; [2049 34737 cnt GeoAsciiOffset]];
GeoAsciiParamsTag = [GeoAsciiParamsTag; code];
GeoAsciiOffset = GeoAsciiOffset + cnt;
end
% Set GeogGeodeticDatumGeoKey, 2050
%
if isfield(option, 'GeogGeodeticDatumGeoKey')
NumberOfKeys = NumberOfKeys + 1;
GeoKeyDirectoryTag(1, 4) = NumberOfKeys;
code = option.GeogGeodeticDatumGeoKey;
GeoKeyDirectoryTag = [GeoKeyDirectoryTag; [2050 0 1 code]];
end
% Set GeogPrimeMeridianGeoKey, 2051
%
if isfield(option, 'GeogPrimeMeridianGeoKey')
NumberOfKeys = NumberOfKeys + 1;
GeoKeyDirectoryTag(1, 4) = NumberOfKeys;
code = option.GeogPrimeMeridianGeoKey;
GeoKeyDirectoryTag = [GeoKeyDirectoryTag; [2051 0 1 code]];
end
% Set GeogLinearUnitsGeoKey, 2052
%
if isfield(option, 'GeogLinearUnitsGeoKey')
NumberOfKeys = NumberOfKeys + 1;
GeoKeyDirectoryTag(1, 4) = NumberOfKeys;
code = option.GeogLinearUnitsGeoKey;
GeoKeyDirectoryTag = [GeoKeyDirectoryTag; [2052 0 1 code]];
end
% Set GeogLinearUnitSizeGeoKey, 2053
%
if isfield(option, 'GeogLinearUnitSizeGeoKey')
NumberOfKeys = NumberOfKeys + 1;
GeoKeyDirectoryTag(1, 4) = NumberOfKeys;
code = option.GeogLinearUnitSizeGeoKey;
GeoKeyDirectoryTag = [GeoKeyDirectoryTag; [2053 34736 1 GeoDoubleOffset]];
GeoDoubleParamsTag = [GeoDoubleParamsTag; code];
GeoDoubleOffset = GeoDoubleOffset + 1;
end
% Set GeogAngularUnitsGeoKey, 2054
%
if isfield(option, 'GeogAngularUnitsGeoKey')
NumberOfKeys = NumberOfKeys + 1;
GeoKeyDirectoryTag(1, 4) = NumberOfKeys;
code = option.GeogAngularUnitsGeoKey;
GeoKeyDirectoryTag = [GeoKeyDirectoryTag; [2054 0 1 code]];
end
% Set GeogAngularUnitSizeGeoKey, 2055
%
if isfield(option, 'GeogAngularUnitSizeGeoKey')
NumberOfKeys = NumberOfKeys + 1;
GeoKeyDirectoryTag(1, 4) = NumberOfKeys;
code = option.GeogAngularUnitSizeGeoKey;
GeoKeyDirectoryTag = [GeoKeyDirectoryTag; [2055 34736 1 GeoDoubleOffset]];
GeoDoubleParamsTag = [GeoDoubleParamsTag; code];
GeoDoubleOffset = GeoDoubleOffset + 1;
end
% Set GeogEllipsoidGeoKey, 2056
%
if isfield(option, 'GeogEllipsoidGeoKey')
NumberOfKeys = NumberOfKeys + 1;
GeoKeyDirectoryTag(1, 4) = NumberOfKeys;
code = option.GeogEllipsoidGeoKey;
GeoKeyDirectoryTag = [GeoKeyDirectoryTag; [2056 0 1 code]];
end
% Set GeogSemiMajorAxisGeoKey, 2057
%
if isfield(option, 'GeogSemiMajorAxisGeoKey')
NumberOfKeys = NumberOfKeys + 1;
GeoKeyDirectoryTag(1, 4) = NumberOfKeys;
code = option.GeogSemiMajorAxisGeoKey;
GeoKeyDirectoryTag = [GeoKeyDirectoryTag; [2057 34736 1 GeoDoubleOffset]];
GeoDoubleParamsTag = [GeoDoubleParamsTag; code];
GeoDoubleOffset = GeoDoubleOffset + 1;
end
% Set GeogSemiMinorAxisGeoKey, 2058
%
if isfield(option, 'GeogSemiMinorAxisGeoKey')
NumberOfKeys = NumberOfKeys + 1;
GeoKeyDirectoryTag(1, 4) = NumberOfKeys;
code = option.GeogSemiMinorAxisGeoKey;
GeoKeyDirectoryTag = [GeoKeyDirectoryTag; [2058 34736 1 GeoDoubleOffset]];
GeoDoubleParamsTag = [GeoDoubleParamsTag; code];
GeoDoubleOffset = GeoDoubleOffset + 1;
end
% Set GeogInvFlatteningGeoKey, 2059
%
if isfield(option, 'GeogInvFlatteningGeoKey')
NumberOfKeys = NumberOfKeys + 1;
GeoKeyDirectoryTag(1, 4) = NumberOfKeys;
code = option.GeogInvFlatteningGeoKey;
GeoKeyDirectoryTag = [GeoKeyDirectoryTag; [2059 34736 1 GeoDoubleOffset]];
GeoDoubleParamsTag = [GeoDoubleParamsTag; code];
GeoDoubleOffset = GeoDoubleOffset + 1;
end
% Set GeogAzimuthUnitsGeoKey, 2060
%
if isfield(option, 'GeogAzimuthUnitsGeoKey')
NumberOfKeys = NumberOfKeys + 1;
GeoKeyDirectoryTag(1, 4) = NumberOfKeys;
code = option.GeogAzimuthUnitsGeoKey;
GeoKeyDirectoryTag = [GeoKeyDirectoryTag; [2060 0 1 code]];
end
% Set GeogPrimeMeridianLongGeoKey, 2061
%
if isfield(option, 'GeogPrimeMeridianLongGeoKey')
NumberOfKeys = NumberOfKeys + 1;
GeoKeyDirectoryTag(1, 4) = NumberOfKeys;
code = option.GeogPrimeMeridianLongGeoKey;
GeoKeyDirectoryTag = [GeoKeyDirectoryTag; [2061 34736 1 GeoDoubleOffset]];
GeoDoubleParamsTag = [GeoDoubleParamsTag; code];
GeoDoubleOffset = GeoDoubleOffset + 1;
end
% Set ProjectedCSTypeGeoKey, 3072
%
if isfield(option, 'ProjectedCSTypeGeoKey')
NumberOfKeys = NumberOfKeys + 1;
GeoKeyDirectoryTag(1, 4) = NumberOfKeys;
code = option.ProjectedCSTypeGeoKey;
GeoKeyDirectoryTag = [GeoKeyDirectoryTag; [3072 0 1 code]];
end
% Set PCSCitationGeoKey, 3073
%
if isfield(option, 'PCSCitationGeoKey')
NumberOfKeys = NumberOfKeys + 1;
GeoKeyDirectoryTag(1, 4) = NumberOfKeys;
code = [option.PCSCitationGeoKey(:); '|']; cnt = length(code);
GeoKeyDirectoryTag = [GeoKeyDirectoryTag; [3073 34737 cnt GeoAsciiOffset]];
GeoAsciiParamsTag = [GeoAsciiParamsTag; code];
GeoAsciiOffset = GeoAsciiOffset + cnt;
end
% Set ProjectionGeoKey, 3074
%
if isfield(option, 'ProjectionGeoKey')
NumberOfKeys = NumberOfKeys + 1;
GeoKeyDirectoryTag(1, 4) = NumberOfKeys;
code = option.ProjectionGeoKey;
GeoKeyDirectoryTag = [GeoKeyDirectoryTag; [3074 0 1 code]];
end
% Set ProjCoordTransGeoKey, 3075
%
if isfield(option, 'ProjCoordTransGeoKey')
NumberOfKeys = NumberOfKeys + 1;
GeoKeyDirectoryTag(1, 4) = NumberOfKeys;
code = option.ProjCoordTransGeoKey;
GeoKeyDirectoryTag = [GeoKeyDirectoryTag; [3075 0 1 code]];
end
% Set ProjLinearUnitsGeoKey, 3076
%
if isfield(option, 'ProjLinearUnitsGeoKey')
NumberOfKeys = NumberOfKeys + 1;
GeoKeyDirectoryTag(1, 4) = NumberOfKeys;
code = option.ProjLinearUnitsGeoKey;
GeoKeyDirectoryTag = [GeoKeyDirectoryTag; [3076 0 1 code]];
end
% Set ProjLinearUnitSizeGeoKey, 3077
%
if isfield(option, 'ProjLinearUnitSizeGeoKey')
NumberOfKeys = NumberOfKeys + 1;
GeoKeyDirectoryTag(1, 4) = NumberOfKeys;
code = option.ProjLinearUnitSizeGeoKey;
GeoKeyDirectoryTag = [GeoKeyDirectoryTag; [3077 34736 1 GeoDoubleOffset]];
GeoDoubleParamsTag = [GeoDoubleParamsTag; code];
GeoDoubleOffset = GeoDoubleOffset + 1;
end
% Set ProjStdParallel1GeoKey, 3078
%
if isfield(option, 'ProjStdParallel1GeoKey')
NumberOfKeys = NumberOfKeys + 1;
GeoKeyDirectoryTag(1, 4) = NumberOfKeys;
code = option.ProjStdParallel1GeoKey;
GeoKeyDirectoryTag = [GeoKeyDirectoryTag; [3078 34736 1 GeoDoubleOffset]];
GeoDoubleParamsTag = [GeoDoubleParamsTag; code];
GeoDoubleOffset = GeoDoubleOffset + 1;
end
% Set ProjStdParallel2GeoKey, 3079
%
if isfield(option, 'ProjStdParallel2GeoKey')
NumberOfKeys = NumberOfKeys + 1;
GeoKeyDirectoryTag(1, 4) = NumberOfKeys;
code = option.ProjStdParallel2GeoKey;
GeoKeyDirectoryTag = [GeoKeyDirectoryTag; [3079 34736 1 GeoDoubleOffset]];
GeoDoubleParamsTag = [GeoDoubleParamsTag; code];
GeoDoubleOffset = GeoDoubleOffset + 1;
end
% Set ProjNatOriginLongGeoKey, 3080
%
if isfield(option, 'ProjNatOriginLongGeoKey')
NumberOfKeys = NumberOfKeys + 1;
GeoKeyDirectoryTag(1, 4) = NumberOfKeys;
code = option.ProjNatOriginLongGeoKey;
GeoKeyDirectoryTag = [GeoKeyDirectoryTag; [3080 34736 1 GeoDoubleOffset]];
GeoDoubleParamsTag = [GeoDoubleParamsTag; code];
GeoDoubleOffset = GeoDoubleOffset + 1;
end
% Set ProjNatOriginLatGeoKey, 3081
%
if isfield(option, 'ProjNatOriginLatGeoKey')
NumberOfKeys = NumberOfKeys + 1;
GeoKeyDirectoryTag(1, 4) = NumberOfKeys;
code = option.ProjNatOriginLatGeoKey;
GeoKeyDirectoryTag = [GeoKeyDirectoryTag; [3081 34736 1 GeoDoubleOffset]];
GeoDoubleParamsTag = [GeoDoubleParamsTag; code];
GeoDoubleOffset = GeoDoubleOffset + 1;
end
% Set ProjFalseEastingGeoKey, 3082
%
if isfield(option, 'ProjFalseEastingGeoKey')
NumberOfKeys = NumberOfKeys + 1;
GeoKeyDirectoryTag(1, 4) = NumberOfKeys;
code = option.ProjFalseEastingGeoKey;
GeoKeyDirectoryTag = [GeoKeyDirectoryTag; [3082 34736 1 GeoDoubleOffset]];
GeoDoubleParamsTag = [GeoDoubleParamsTag; code];
GeoDoubleOffset = GeoDoubleOffset + 1;
end
% Set ProjFalseNorthingGeoKey, 3083
%
if isfield(option, 'ProjFalseNorthingGeoKey')
NumberOfKeys = NumberOfKeys + 1;
GeoKeyDirectoryTag(1, 4) = NumberOfKeys;
code = option.ProjFalseNorthingGeoKey;
GeoKeyDirectoryTag = [GeoKeyDirectoryTag; [3083 34736 1 GeoDoubleOffset]];
GeoDoubleParamsTag = [GeoDoubleParamsTag; code];
GeoDoubleOffset = GeoDoubleOffset + 1;
end
% Set ProjFalseOriginLongGeoKey, 3084
%
if isfield(option, 'ProjFalseOriginLongGeoKey')
NumberOfKeys = NumberOfKeys + 1;
GeoKeyDirectoryTag(1, 4) = NumberOfKeys;
code = option.ProjFalseOriginLongGeoKey;
GeoKeyDirectoryTag = [GeoKeyDirectoryTag; [3084 34736 1 GeoDoubleOffset]];
GeoDoubleParamsTag = [GeoDoubleParamsTag; code];
GeoDoubleOffset = GeoDoubleOffset + 1;
end
% Set ProjFalseOriginLatGeoKey, 3085
%
if isfield(option, 'ProjFalseOriginLatGeoKey')
NumberOfKeys = NumberOfKeys + 1;
GeoKeyDirectoryTag(1, 4) = NumberOfKeys;
code = option.ProjFalseOriginLatGeoKey;
GeoKeyDirectoryTag = [GeoKeyDirectoryTag; [3085 34736 1 GeoDoubleOffset]];
GeoDoubleParamsTag = [GeoDoubleParamsTag; code];
GeoDoubleOffset = GeoDoubleOffset + 1;
end
% Set ProjFalseOriginEastingGeoKey, 3086
%
if isfield(option, 'ProjFalseOriginEastingGeoKey')
NumberOfKeys = NumberOfKeys + 1;
GeoKeyDirectoryTag(1, 4) = NumberOfKeys;
code = option.ProjFalseOriginEastingGeoKey;
GeoKeyDirectoryTag = [GeoKeyDirectoryTag; [3086 34736 1 GeoDoubleOffset]];
GeoDoubleParamsTag = [GeoDoubleParamsTag; code];
GeoDoubleOffset = GeoDoubleOffset + 1;
end
% Set ProjFalseOriginNorthingGeoKey, 3087
%
if isfield(option, 'ProjFalseOriginNorthingGeoKey')
NumberOfKeys = NumberOfKeys + 1;
GeoKeyDirectoryTag(1, 4) = NumberOfKeys;
code = option.ProjFalseOriginNorthingGeoKey;
GeoKeyDirectoryTag = [GeoKeyDirectoryTag; [3087 34736 1 GeoDoubleOffset]];
GeoDoubleParamsTag = [GeoDoubleParamsTag; code];
GeoDoubleOffset = GeoDoubleOffset + 1;
end
% Set ProjCenterLongGeoKey, 3088
%
if isfield(option, 'ProjCenterLongGeoKey')
NumberOfKeys = NumberOfKeys + 1;
GeoKeyDirectoryTag(1, 4) = NumberOfKeys;
code = option.ProjCenterLongGeoKey;
GeoKeyDirectoryTag = [GeoKeyDirectoryTag; [3088 34736 1 GeoDoubleOffset]];
GeoDoubleParamsTag = [GeoDoubleParamsTag; code];
GeoDoubleOffset = GeoDoubleOffset + 1;
end
% Set ProjCenterLatGeoKey, 3089
%
if isfield(option, 'ProjCenterLatGeoKey')
NumberOfKeys = NumberOfKeys + 1;
GeoKeyDirectoryTag(1, 4) = NumberOfKeys;
code = option.ProjCenterLatGeoKey;
GeoKeyDirectoryTag = [GeoKeyDirectoryTag; [3089 34736 1 GeoDoubleOffset]];
GeoDoubleParamsTag = [GeoDoubleParamsTag; code];
GeoDoubleOffset = GeoDoubleOffset + 1;
end
% Set ProjCenterEastingGeoKey, 3090
%
if isfield(option, 'ProjCenterEastingGeoKey')
NumberOfKeys = NumberOfKeys + 1;
GeoKeyDirectoryTag(1, 4) = NumberOfKeys;
code = option.ProjCenterEastingGeoKey;
GeoKeyDirectoryTag = [GeoKeyDirectoryTag; [3090 34736 1 GeoDoubleOffset]];
GeoDoubleParamsTag = [GeoDoubleParamsTag; code];
GeoDoubleOffset = GeoDoubleOffset + 1;
end
% Set ProjCenterOriginNorthingGeoKey, 3091
%
if isfield(option, 'ProjCenterOriginNorthingGeoKey')
NumberOfKeys = NumberOfKeys + 1;
GeoKeyDirectoryTag(1, 4) = NumberOfKeys;
code = option.ProjCenterOriginNorthingGeoKey;
GeoKeyDirectoryTag = [GeoKeyDirectoryTag; [3091 34736 1 GeoDoubleOffset]];
GeoDoubleParamsTag = [GeoDoubleParamsTag; code];
GeoDoubleOffset = GeoDoubleOffset + 1;
end
% Set ProjScaleAtNatOriginGeoKey, 3092
%
if isfield(option, 'ProjScaleAtNatOriginGeoKey')
NumberOfKeys = NumberOfKeys + 1;
GeoKeyDirectoryTag(1, 4) = NumberOfKeys;
code = option.ProjScaleAtNatOriginGeoKey;
GeoKeyDirectoryTag = [GeoKeyDirectoryTag; [3092 34736 1 GeoDoubleOffset]];
GeoDoubleParamsTag = [GeoDoubleParamsTag; code];
GeoDoubleOffset = GeoDoubleOffset + 1;
end
% Set ProjScaleAtCenterGeoKey, 3093
%
if isfield(option, 'ProjScaleAtCenterGeoKey')
NumberOfKeys = NumberOfKeys + 1;
GeoKeyDirectoryTag(1, 4) = NumberOfKeys;
code = option.ProjScaleAtCenterGeoKey;
GeoKeyDirectoryTag = [GeoKeyDirectoryTag; [3093 34736 1 GeoDoubleOffset]];
GeoDoubleParamsTag = [GeoDoubleParamsTag; code];
GeoDoubleOffset = GeoDoubleOffset + 1;
end
% Set ProjAzimuthAngleGeoKey, 3094
%
if isfield(option, 'ProjAzimuthAngleGeoKey')
NumberOfKeys = NumberOfKeys + 1;
GeoKeyDirectoryTag(1, 4) = NumberOfKeys;
code = option.ProjAzimuthAngleGeoKey;
GeoKeyDirectoryTag = [GeoKeyDirectoryTag; [3094 34736 1 GeoDoubleOffset]];
GeoDoubleParamsTag = [GeoDoubleParamsTag; code];
GeoDoubleOffset = GeoDoubleOffset + 1;
end
% Set ProjStraightVertPoleLongGeoKey, 3095
%
if isfield(option, 'ProjStraightVertPoleLongGeoKey')
NumberOfKeys = NumberOfKeys + 1;
GeoKeyDirectoryTag(1, 4) = NumberOfKeys;
code = option.ProjStraightVertPoleLongGeoKey;
GeoKeyDirectoryTag = [GeoKeyDirectoryTag; [3095 34736 1 GeoDoubleOffset]];
GeoDoubleParamsTag = [GeoDoubleParamsTag; code];
GeoDoubleOffset = GeoDoubleOffset + 1;
end
% Set VerticalCSTypeGeoKey, 4096
%
if isfield(option, 'VerticalCSTypeGeoKey')
NumberOfKeys = NumberOfKeys + 1;
GeoKeyDirectoryTag(1, 4) = NumberOfKeys;
code = option.VerticalCSTypeGeoKey;
GeoKeyDirectoryTag = [GeoKeyDirectoryTag; [4096 0 1 code]];
end
% Set VerticalCitationGeoKey, 4097
%
if isfield(option, 'VerticalCitationGeoKey')
NumberOfKeys = NumberOfKeys + 1;
GeoKeyDirectoryTag(1, 4) = NumberOfKeys;
code = [option.VerticalCitationGeoKey(:); '|']; cnt = length(code);
GeoKeyDirectoryTag = [GeoKeyDirectoryTag; [4097 34737 cnt GeoAsciiOffset]];
GeoAsciiParamsTag = [GeoAsciiParamsTag; code];
GeoAsciiOffset = GeoAsciiOffset + cnt;
end
% Set VerticalUnitsGeoKey, 4099
%
if isfield(option, 'VerticalUnitsGeoKey')
NumberOfKeys = NumberOfKeys + 1;
GeoKeyDirectoryTag(1, 4) = NumberOfKeys;
code = option.VerticalUnitsGeoKey;
GeoKeyDirectoryTag = [GeoKeyDirectoryTag; [4099 0 1 code]];
end
% Set ifd.ImageWidth, 256
%
ifd.ImageWidth = size(image, 2);
% Set ifd.ImageLength, 257
%
ifd.ImageLength = size(image, 1);
% Set ifd.BitsPerSample, 258
%
ifd.BitsPerSample = abs(bit_depth)*ones(size(image,3),1);
% Set ifd.Compression, 259
%
ifd.Compression = 1; % Uncompressed
% Set ifd.PhotometricInterpretation, 262
%
if isfield(option, 'ColorMap')
ifd.PhotometricInterpretation = 3; % Palette color
elseif isfield(option, 'FullColor') & option.FullColor == 1
ifd.PhotometricInterpretation = 2; % RGB
else
ifd.PhotometricInterpretation = 1; % BlackIsZero
end
% Set ifd.StripOffsets, 273
%
ifd.StripOffsets = ones(ifd.ImageLength,1) * ifd.ImageWidth ...
* abs(bit_depth)/8 * size(image,3); % Num of bytes in a strip
ifd.StripOffsets = cumsum(ifd.StripOffsets);
ifd.StripOffsets(end) = [];
ifd.StripOffsets = [0; ifd.StripOffsets];
ifd.StripOffsets = ifd.StripOffsets + 8; % 8 bytes before 1st strip
% Set ifd.SamplesPerPixel, 277
%
ifd.SamplesPerPixel = size(image, 3); % Gray level intensity
% Set ifd.RowsPerStrip, 278
%
ifd.RowsPerStrip = 1; % rows per strip
% Set ifd.StripByteCounts, 279
%
ifd.StripByteCounts = ones(ifd.ImageLength,1) * ifd.ImageWidth ...
* abs(bit_depth)/8 * size(image,3); % Num of bytes in a strip
% Set ifd.XResolution, 282
%
ifd.XResolution = 96; % 96 dpi
% Set ifd.YResolution, 283
%
ifd.YResolution = 96; % 96 dpi
% Set ifd.PlanarConfiguration, 284
%
ifd.PlanarConfiguration = 1; % Chunky
% Set ifd.ResolutionUnit, 296
%
ifd.ResolutionUnit = 2; % Inch
% Set ifd.ColorMap, 320
%
if isfield(option, 'ColorMap')
ifd.ColorMap = option.ColorMap(:);
end
% Set ifd.ExtraSamples, 338
%
if ndims(image) > 2
ifd.ExtraSamples = zeros((size(image,3)-1),1);
end
% Set ifd.SampleFormat, 339
%
switch bit_depth
case {1, -1}
ifd.SampleFormat = 4*ones(size(image,3),1); % Undefined
case {8, -16}
ifd.SampleFormat = ones(size(image,3),1); % Unsigned integer
case {16, -32}
ifd.SampleFormat = 2*ones(size(image,3),1); % Signed integer
case 32
ifd.SampleFormat = 3*ones(size(image,3),1); % Floating point
end
if isfield(ifd, 'ColorMap') & bit_depth == 16
ifd.SampleFormat = ones(size(image,3),1); % Unsigned integer
end
if ~isempty(bbox)
% Set gifd.ModelPixelScaleTag, 33550
%
if gifd.GTRasterTypeGeoKey == 1
Sx = (bbox(2) - bbox(1)) / ifd.ImageWidth;
Sy = (bbox(4) - bbox(3)) / ifd.ImageLength;
elseif gifd.GTRasterTypeGeoKey == 2
Sx = (bbox(2) - bbox(1)) / (ifd.ImageWidth-1);
Sy = (bbox(4) - bbox(3)) / (ifd.ImageLength-1);
else
error('Incorrect option.GTRasterTypeGeoKey value');
end
gifd.ModelPixelScaleTag = [Sx Sy 0];
% Set gifd.ModelTiepointTag, 33922
%
RasterTiepoint = [0 0 0];
ModelTiepoint = [bbox(1) bbox(4) 0];
gifd.ModelTiepointTag = [RasterTiepoint ModelTiepoint];
else
if isfield(option, 'ModelPixelScaleTag')
gifd.ModelPixelScaleTag = option.ModelPixelScaleTag(:);
else
% error('ModelPixelScale is required when your Model Type is not Geographic.');
end
if isfield(option, 'ModelTiepointTag')
gifd.ModelTiepointTag = option.ModelTiepointTag(:);
end
if isfield(option, 'ModelTransformationTag')
gifd.ModelTransformationTag = option.ModelTransformationTag(:);
end
if ~isfield(option, 'ModelTiepointTag') & ~isfield(option, 'ModelTransformationTag')
error('Either ModelTiepoint or ModelTransformation is required when your Model Type is not Geographic.');
end
end
% Write TIFF info
%
num_entry = 16; % 16 basic entries
if isfield(ifd, 'ColorMap')
num_entry = num_entry + 1;
end
if isfield(ifd, 'ExtraSamples')
num_entry = num_entry + 1;
end
if isfield(gifd, 'ModelPixelScaleTag')
num_entry = num_entry + 1;
end
if isfield(gifd, 'ModelTiepointTag')
num_entry = num_entry + 1;
end
if isfield(gifd, 'ModelTransformationTag')
num_entry = num_entry + 1;
end
if GeoAsciiOffset > 0
num_entry = num_entry + 1;
end
if GeoDoubleOffset > 0
num_entry = num_entry + 1;
end
ifd_end = 8 + sum(ifd.StripByteCounts) + 2 + num_entry*12 + 4;
fid = fopen(filename, 'wb', 'ieee-le');
if fid < 0,
msg = sprintf('Cannot create file "%s".', filename);
error(msg);
end
fwrite(fid, 'II', 'uint8'); % Little-endian
fwrite(fid, 42, 'uint16'); % TIFF signature
fwrite(fid, (8 + sum(ifd.StripByteCounts)), 'uint32'); % IFD offset
% Write data strip
%
image = squeeze(permute(image, [3 2 1]));
switch abs(bit_depth)
case 1
image = bit2byte(image);
fwrite(fid, image, 'uint8');
case 8
fwrite(fid, image, 'uint8');
case 16
if isfield(ifd, 'ColorMap') | bit_depth == -16
fwrite(fid, image, 'uint16');
else
fwrite(fid, image, 'int16');
end
case 32
if bit_depth == -32
fwrite(fid, image, 'int32');
else
fwrite(fid, image, 'single');
end
end
fwrite(fid, num_entry, 'uint16');
% Entry 1
%
fwrite(fid, 256, 'uint16');
fwrite(fid, 3, 'uint16'); % uint16
fwrite(fid, 1, 'uint32'); % count
fwrite(fid, ifd.ImageWidth, 'uint16');
fwrite(fid, 0, 'uint16');
% Entry 2
%
fwrite(fid, 257, 'uint16');
fwrite(fid, 3, 'uint16'); % uint16
fwrite(fid, 1, 'uint32'); % count
fwrite(fid, ifd.ImageLength, 'uint16');
fwrite(fid, 0, 'uint16');
% Entry 3
%
fwrite(fid, 258, 'uint16');
fwrite(fid, 3, 'uint16'); % uint16
fwrite(fid, length(ifd.BitsPerSample), 'uint32'); % count
if length(ifd.BitsPerSample) > 2
fwrite(fid, ifd_end, 'uint32');
ifd_end = ifd_end + length(ifd.BitsPerSample)*2;
else
fwrite(fid, ifd.BitsPerSample(1), 'uint16');
if length(ifd.BitsPerSample) > 1
fwrite(fid, ifd.BitsPerSample(2), 'uint16');
else
fwrite(fid, 0, 'uint16');
end
end
% Entry 4
%
fwrite(fid, 259, 'uint16');
fwrite(fid, 3, 'uint16'); % uint16
fwrite(fid, 1, 'uint32'); % count
fwrite(fid, ifd.Compression, 'uint16');
fwrite(fid, 0, 'uint16');
% Entry 5
%
fwrite(fid, 262, 'uint16');
fwrite(fid, 3, 'uint16'); % uint16
fwrite(fid, 1, 'uint32'); % count
fwrite(fid, ifd.PhotometricInterpretation, 'uint16');
fwrite(fid, 0, 'uint16');
% Entry 6
%
fwrite(fid, 273, 'uint16'); % StripOffsets
fwrite(fid, 4, 'uint16'); % uint32
fwrite(fid, ifd.ImageLength, 'uint32'); % count
fwrite(fid, ifd_end, 'uint32');
ifd_end = ifd_end + ifd.ImageLength*4;
% Entry 7
%
fwrite(fid, 274, 'uint16');
fwrite(fid, 3, 'uint16'); % uint16
fwrite(fid, 1, 'uint32'); % count
fwrite(fid, ifd.Orientation, 'uint16');
fwrite(fid, 0, 'uint16');
% Entry 8
%
fwrite(fid, 277, 'uint16');
fwrite(fid, 3, 'uint16'); % uint16
fwrite(fid, 1, 'uint32'); % count
fwrite(fid, ifd.SamplesPerPixel, 'uint16');
fwrite(fid, 0, 'uint16');
% Entry 9
%
fwrite(fid, 278, 'uint16');
fwrite(fid, 3, 'uint16'); % uint16
fwrite(fid, 1, 'uint32'); % count
fwrite(fid, ifd.RowsPerStrip, 'uint16');
fwrite(fid, 0, 'uint16');
% Entry 10
%
fwrite(fid, 279, 'uint16'); % StripByteCounts
fwrite(fid, 4, 'uint16'); % uint32
fwrite(fid, ifd.ImageLength, 'uint32'); % count
fwrite(fid, ifd_end, 'uint32');
ifd_end = ifd_end + ifd.ImageLength*4;
% Entry 11
%
fwrite(fid, 282, 'uint16'); % XResolution
fwrite(fid, 5, 'uint16'); % uint16
fwrite(fid, 1, 'uint32'); % count
fwrite(fid, ifd_end, 'uint32');
ifd_end = ifd_end + 2*4;
% Entry 12
%
fwrite(fid, 283, 'uint16'); % YResolution
fwrite(fid, 5, 'uint16'); % uint16
fwrite(fid, 1, 'uint32'); % count
fwrite(fid, ifd_end, 'uint32');
ifd_end = ifd_end + 2*4;
% Entry 13
%
fwrite(fid, 284, 'uint16');
fwrite(fid, 3, 'uint16'); % uint16
fwrite(fid, 1, 'uint32'); % count
fwrite(fid, ifd.PlanarConfiguration, 'uint16');
fwrite(fid, 0, 'uint16');
% Entry 14
%
fwrite(fid, 296, 'uint16');
fwrite(fid, 3, 'uint16'); % uint16
fwrite(fid, 1, 'uint32'); % count
fwrite(fid, ifd.ResolutionUnit, 'uint16');
fwrite(fid, 0, 'uint16');
if isfield(option, 'ColorMap')
fwrite(fid, 320, 'uint16');
fwrite(fid, 3, 'uint16'); % uint16
fwrite(fid, length(ifd.ColorMap), 'uint32'); % count
fwrite(fid, ifd_end, 'uint32');
ifd_end = ifd_end + length(ifd.ColorMap)*2;
end
if isfield(ifd, 'ExtraSamples')
fwrite(fid, 338, 'uint16');
fwrite(fid, 3, 'uint16'); % uint16
fwrite(fid, length(ifd.ExtraSamples), 'uint32'); % count
if length(ifd.ExtraSamples) > 2
fwrite(fid, ifd_end, 'uint32');
ifd_end = ifd_end + length(ifd.ExtraSamples)*2;
else
fwrite(fid, ifd.ExtraSamples(1), 'uint16');
if length(ifd.ExtraSamples) > 1
fwrite(fid, ifd.ExtraSamples(2), 'uint16');
else
fwrite(fid, 0, 'uint16');
end
end
end
% Entry 15
%
fwrite(fid, 339, 'uint16');
fwrite(fid, 3, 'uint16'); % uint16
fwrite(fid, length(ifd.SampleFormat), 'uint32'); % count
if length(ifd.SampleFormat) > 2
fwrite(fid, ifd_end, 'uint32');
ifd_end = ifd_end + length(ifd.SampleFormat)*2;
else
fwrite(fid, ifd.SampleFormat(1), 'uint16');
if length(ifd.SampleFormat) > 1
fwrite(fid, ifd.SampleFormat(2), 'uint16');
else
fwrite(fid, 0, 'uint16');
end
end
if isfield(gifd, 'ModelPixelScaleTag')
fwrite(fid, 33550, 'uint16'); % ModelPixelScaleTag
fwrite(fid, 12, 'uint16'); % double
fwrite(fid, 3, 'uint32'); % count
fwrite(fid, ifd_end, 'uint32');
ifd_end = ifd_end + 3*8;
end
if isfield(gifd, 'ModelTiepointTag')
fwrite(fid, 33922, 'uint16'); % ModelTiepointTag
fwrite(fid, 12, 'uint16'); % double
fwrite(fid, length(gifd.ModelTiepointTag), 'uint32'); % count
fwrite(fid, ifd_end, 'uint32');
ifd_end = ifd_end + length(gifd.ModelTiepointTag)*8;
end
if isfield(gifd, 'ModelTransformationTag')
fwrite(fid, 34264, 'uint16'); % ModelTransformationTag
fwrite(fid, 12, 'uint16'); % double
fwrite(fid, 16, 'uint32'); % count
fwrite(fid, ifd_end, 'uint32');
ifd_end = ifd_end + 16*8;
end
% Entry 16
%
fwrite(fid, 34735, 'uint16'); % GeoKeyDirectoryTag
fwrite(fid, 3, 'uint16'); % uint16
fwrite(fid, prod(size(GeoKeyDirectoryTag)), 'uint32'); % count
fwrite(fid, ifd_end, 'uint32');
ifd_end = ifd_end + prod(size(GeoKeyDirectoryTag))*2;
if GeoDoubleOffset > 0
fwrite(fid, 34736, 'uint16'); % GeoDoubleParamsTag
fwrite(fid, 12, 'uint16'); % double
fwrite(fid, length(GeoDoubleParamsTag), 'uint32'); % count
fwrite(fid, ifd_end, 'uint32');
ifd_end = ifd_end + length(GeoDoubleParamsTag)*8;
end
if GeoAsciiOffset > 0
GeoAsciiParamsTag = [GeoAsciiParamsTag; 0]; % TIFF6 requres NUL ending
fwrite(fid, 34737, 'uint16'); % GeoAsciiParamsTag
fwrite(fid, 2, 'uint16'); % ascii
fwrite(fid, length(GeoAsciiParamsTag), 'uint32'); % count
fwrite(fid, ifd_end, 'uint32');
ifd_end = ifd_end + length(GeoAsciiParamsTag);
end
% IFD is terminated with 4-byte offset to the next IFD,
% or 0 if there are none.
%
fwrite(fid, 0, 'uint32');
if length(ifd.BitsPerSample) > 2
fwrite(fid, ifd.BitsPerSample, 'uint16'); % 258
end
fwrite(fid, ifd.StripOffsets, 'uint32'); % 273
fwrite(fid, ifd.StripByteCounts, 'uint32'); % 279
fwrite(fid, ifd.XResolution, 'uint32'); % 282
fwrite(fid, 1, 'uint32'); % 282
fwrite(fid, ifd.YResolution, 'uint32'); % 283
fwrite(fid, 1, 'uint32'); % 283
if isfield(option, 'ColorMap')
fwrite(fid, ifd.ColorMap, 'uint16'); % 320
end
if isfield(ifd, 'ExtraSamples') & length(ifd.ExtraSamples) > 2
fwrite(fid, ifd.ExtraSamples, 'uint16'); % 338
end
if length(ifd.SampleFormat) > 2
fwrite(fid, ifd.SampleFormat, 'uint16'); % 339
end
if isfield(gifd,'ModelPixelScaleTag')
fwrite(fid, gifd.ModelPixelScaleTag, 'double'); % 33550
end
if isfield(gifd,'ModelTiepointTag')
fwrite(fid, gifd.ModelTiepointTag, 'double'); % 33922
end
if isfield(gifd,'ModelTransformationTag')
fwrite(fid, gifd.ModelTransformationTag, 'double'); % 34264
end
GeoKeyDirectoryTag = GeoKeyDirectoryTag';
fwrite(fid, GeoKeyDirectoryTag, 'uint16'); % 34735
if GeoDoubleOffset > 0
fwrite(fid, GeoDoubleParamsTag, 'double'); % 34736
end
if GeoAsciiOffset > 0
fwrite(fid, GeoAsciiParamsTag, 'uint8'); % 34737
end
fclose(fid);
return; % geotiffwrite
%------------------------------------------------------------------------
% Convert bits to bytes to store
%
function byte = bit2byte(bit)
% byte = bit(:);
% Pad with 0
%
% if mod(length(byte),8) ~= 0
% byte = [byte; zeros(8-mod(length(byte),8),1)];
% end
byte = reshape(bit, [8, length(bit(:))/8])';
byte(find(byte>0)) = 1;
byte(find(byte<1)) = 0;
% byte = byte*[1 2 4 8 16 32 64 128]';
byte = double(byte)*[128 64 32 16 8 4 2 1]';
return; % bit2byte
