Commit e742b18e authored by Daniel Bourgault's avatar Daniel Bourgault

Addition of the function g_geotiff.m to help make geotif figures

parent 4f1faf35
......@@ -25,7 +25,7 @@ hfov = 62.00; % Field of view of the camera
lambda = 53.0; % Dip angle below horizontal (e.g. straight down = 90, horizontal = 0)
phi = 1.0; % Tilt angle (generally close to 0).
H = 1.755; % Camera altitude (m)
theta = 180.0; % View angle anticlockwise from North (e.g. straight East = 270)
theta = 180.0; % View angle clockwise from North (e.g. straight East = 90)
% Uncertainty in parameters. Set the uncertainty to 0.0 for fixed parameters.
dhfov = 5.0;
......
function [] = g_geotiff(LON,LAT,dx,lon_min,lon_max,lat_min,lat_max)
%
% This function creates geotiff images given the longitudes (LON)
% and latitudes (LAT) of every pixels of the associated jpeg images found
% in the current folder. The matrices LON and LAT typically come from the
% g_rect package.
%
% The geotiff images are constructed by interpolating on a regular grid of
% size dx (in m) the irregularly-spaced pixel intensities defined by the
% LON-LAT pair.
%
% INPUT PARAMETERS:
%
% LON: A matrix of size identical to the associated images in the current
% folder that gives the longitude of every pixel of those images.
% This matrix would typically come from the g_rect package.
%
% LAT: Same as LON but for the latitude.
%
% dx: The grid size (in m) of the regularly-spaced interpolated
% geotif image
%
% lon_min: The longitude of the southwest corner of the desired
% geotif image.
%
% lon_max: The longitude of the northeastcorner of the desired
% geotif image.
%
% lat_min: The latitude of the southwest corner of the desired
% geotif image.
%
% lat_max: The latitude of the northeastcorner of the desired
% geotif image.
%
% OUTPUT:
%
% This function writes geotiff images that have the same name as the
% jpeg images but with the extension .tif
%
% LAST REVISION: 5 May 2020
%
% AUTHOR:
% - Daniel Bourgault (daniel_bourgault@uqar.ca)
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Construct the vectors 'lon' and 'lat' that define the interpolating
% grid of constant resolution dx (in m)
lat0 = (lat_min + lat_max)/2;
dlat = dx/(1852*60);
dlon = dlat/cosd(lat0);
lon = [lon_min:dlon:lon_max];
lat = [lat_min:dlat:lat_max];
% Find the indices of only the finite values in the matrix LON
% so that the interpolation is only done on finite values and does
% not spend time interpolating NaNs.
ifinite = find(isfinite(LON) == 1);
% Find all .jpg images in the current folder
img_fname = dir('*.jpg');
N_img = length(img_fname);
% Loop over all images
for i = 1:N_img
img = imread(img_fname(i).name);
% Convert to gray scale
img = rgb2gray(img);
img = double(img);
img = img';
% Interpolate on the regular grid defined above by the vector
% lon and lat
img_interp = griddata(LON(ifinite),LAT(ifinite),img(ifinite),lon,lat');
% Convert real numbers to unsigned integers
img_interp = uint8(img_interp);
% Write the geotiff file
geotif_img_fname = [img_fname(i).name(1:end-3),'tif']
bbox = [lon_min, lat_min; lon_max, lat_max];
geotiffwrite(geotif_img_fname, bbox, flipud(img_interp),8);
end
\ No newline at end of file
......@@ -116,7 +116,7 @@ end
fclose(fid);
%% Import the GCP data (4 column) at the end of the parameter file
gcp = importdata(inputFname,' ',nHeaderLine);
gcp = importdata(inputFname,' ',nHeaderLine)
i_gcp = gcp.data(:,1);
j_gcp = gcp.data(:,2);
lon_gcp = gcp.data(:,3);
......
%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]];