-
Notifications
You must be signed in to change notification settings - Fork 58
/
geotiffread_lcog_matlab2020a.m
179 lines (158 loc) · 6.63 KB
/
geotiffread_lcog_matlab2020a.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
function varargout = geotiffread_lcog_matlab2020a(filename, varargin)
%GEOTIFFREAD Read GeoTIFF file
%
% GEOTIFFREAD is not recommended. Use READGEORASTER instead.
%
% [A, R] = GEOTIFFREAD(FILENAME) reads a georeferenced grayscale, RGB, or
% multispectral image or data grid from the GeoTIFF file specified by
% FILENAME into A and constructs a spatial referencing object, R.
%
% [X, CMAP, R] = GEOTIFFREAD(FILENAME) reads an indexed image into X and
% the associated colormap into CMAP, and constructs a spatial referencing
% object, R. Colormap values in the image file are rescaled into the
% range [0,1].
%
% FILENAME is a string scalar or character vector that specifies the name
% of the GeoTIFF file. FILENAME can include the folder name. Otherwise,
% the file must be in the current folder or in a folder on the MATLAB
% path. If the named file includes the extension '.TIF' or '.TIFF'
% (either upper or lower case), you can omit the extension from FILENAME.
%
% A is a two-dimensional array if the file contains a grayscale image or
% data grid. A is an M-by-N-by-P array if the file contains a color
% image, multispectral image, hyperspectral image, or data grid. The
% class of A depends on the storage class of the pixel data in the file
% which is related to the BitsPerSample property as returned by the
% IMFINFO function.
%
% R is a geographic raster reference object if the image or data grid is
% referenced to a geographic coordinate system, or a map raster reference
% object if it is referenced to a projected coordinate system.
%
% [A, REFMAT, BBOX] = GEOTIFFREAD(FILENAME) reads a georeferenced
% grayscale, RGB, or multispectral image or data grid into A; the
% corresponding referencing matrix into REFMAT; and the bounding box into
% BBOX.
%
% [X, CMAP, REFMAT, BBOX] = GEOTIFFREAD(FILENAME) reads an indexed image
% into X, the associated colormap into CMAP, the referencing matrix into
% REFMAT, and the bounding box into BBOX. The referencing matrix must
% be unambiguously defined by the GeoTIFF file, otherwise it and the
% bounding box are returned empty.
%
% [___] = GEOTIFFREAD(URL) reads the GeoTIFF image from a URL. The
% URL must include the protocol type (e.g., "https://").
%
% [___] = GEOTIFFREAD(___, IDX) reads one image from a multi-image
% GeoTIFF file (or URL). IDX is an integer value that specifies the order
% that the image appears in the file. For example, if IDX is 3,
% GEOTIFFREAD reads the third image in the file. If you omit this
% argument, GEOTIFFREAD reads the first image.
%
% Note
% ----
% GEOTIFFREAD imports pixel data using the TIFF-reading capabilities of
% the MATLAB function IMREAD and likewise shares any limitations of
% IMREAD. Consult the IMREAD documentation for specific information on
% TIFF image support.
%
% Example
% -------
% % Read and display the Boston GeoTIFF image.
% % Includes material (c) GeoEye, all rights reserved.
% [boston, R] = geotiffread('boston.tif');
% figure
% mapshow(boston, R)
% axis image off
%
% See also GEOSHOW, GEOTIFFINFO, GEOTIFFWRITE, IMREAD, MAPSHOW, READGEORASTER
% Copyright 1996-2019 The MathWorks, Inc.
% Verify the input and output argument count.
narginchk(1,2);
nargoutchk(0,4);
% Parse the inputs.
filename = convertStringsToChars(filename);
[filename, url, idx] = parseInputs(filename, varargin);
% Read the info fields from the filename.
info = geotiffinfo_lcog_matlab2020a(filename);
% Read the image from the filename.
[A, cmap] = imread(filename,idx);
% If the RefMatrix is empty, try to obtain the spatial information from a
% corresponding worldfile.
if isempty(info.RefMatrix)
worldfilename = getworldfilename(filename);
info = getSpatialInfoFromWorldfile(worldfilename, info, size(A));
end
% Delete temporary file from Internet download.
if (url)
deleteDownload(filename);
end
% Assign output arguments.
varargout = assignOutputArguments(A, cmap, info, nargout);
%--------------------------------------------------------------------------
function [filename, url, idx] = parseInputs(filename, inputs)
% Parse the inputs from the cell array, INPUTS.
% Verify the filename and obtain the full pathname.
extensions = {'tif', 'tiff'};
[filename, url] = internal.map.checkfilename(filename, extensions, mfilename, 1, true);
% Check and set the image index number.
if ~isempty(inputs)
idx = inputs{1};
attributes = {'real' 'scalar' 'positive'};
validateattributes(idx, {'numeric'}, attributes, mfilename, 'IDX', 2);
else
idx = 1;
end
%--------------------------------------------------------------------------
function info = getSpatialInfoFromWorldfile(worldfilename, info, rasterSize)
% Obtain the referencing matrix from a world file, if it exits. If so,
% compute the bounding box and construct a spatial referencing object from
% the referencing matrix and update the fields of the INFO structure.
% WORLDFILENAME is a string denoting the name of the world file.
if exist(worldfilename,'file')
% Obtain the referencing matrix from the world file.
refmat = worldfileread(worldfilename);
% Calculate the spatial referencing object and bounding box from the
% referencing matrix if it is not empty.
if ~isempty(refmat)
if strcmp(info.ModelType, 'ModelTypeGeographic')
R = refmatToGeoRasterReference(refmat, rasterSize);
elseif strcmp(info.ModelType, 'ModelTypeProjected')
R = refmatToMapRasterReference(refmat, rasterSize);
else
R = [];
end
info.BoundingBox = mapbbox(refmat, rasterSize);
info.RefMatrix = refmat;
info.SpatialRef = R;
end
end
%--------------------------------------------------------------------------
function outputs = assignOutputArguments(A, cmap, info, numOutputs)
% Assign the output arguments based on the number of arguments requested.
outputs{1} = A;
switch numOutputs
case 2
if strcmp(info.ColorType, 'indexed')
% [X, CMAP] = GEOTIFFREAD(...)
outputs{2} = cmap;
else
% [A, R] = GEOTIFFREAD(...)
outputs{2} = info.SpatialRef;
end
case 3
if strcmp(info.ColorType, 'indexed')
% [X, CMAP, R] = GEOTIFFREAD(...)
outputs{2} = cmap;
outputs{3} = info.SpatialRef;
else
% [A, REFMAT, BBOX] = GEOTIFFREAD(...)
outputs{2} = info.RefMatrix;
outputs{3} = info.BoundingBox;
end
case 4
% [X, CMAP, REFMAT, BBOX] = GEOTIFFREAD(...)
outputs{2} = cmap;
outputs{3} = info.RefMatrix;
outputs{4} = info.BoundingBox;
end