-
Notifications
You must be signed in to change notification settings - Fork 3
/
tinitalh.m
124 lines (116 loc) · 3.66 KB
/
tinitalh.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
function varargout=tinitalh(dirp,diro,xver)
% [hdr,TV,TN,TA,bx,by,tt,tl]=TINITALH(dirp,diro,xver)
%
% Gets and displays all headers inside a TINITALY directory ready to feed
% into TINITALG which will get and display data and grids.
%
% INPUT:
%
% dirp Subdirectory [e.g. 'DATA'] of:
% diro Main directory [e.g. '/u/fjsimonsIFILES/TOPOGRAPHY/ITALY/TINITALY']
% xver 2 Provides a graphical test [default]
% 0 Does not provide a graphical test
%
% OUTPUT:
%
% hdr All the header name strings in a cell array
% TV All the header variables, in a cell
% TN The header variable names, as a matrix
% TA All the the header variables, in a matrix
% bx,by All the box corners, if you like
% tt,tl Object handles to the box titles and the overall title
%
% EXAMPLE:
%
% tinitalh([],[],2) % will bring up the map with available tiles
%
% Last modified by fjsimons-at-alum.mit.edu, 11/18/2020
% Bottom-level directory name, taken from the Tinitaly download
defval('dirp','DATA')
% Top-level directory name, where you keep the Tinitaly directory
defval('diro','/u/fjsimons/IFILES/TOPOGRAPHY/ITALY/TINITALY')
% Graphical checking of grid parameters
defval('xver',2)
% Get all the headers listed
try
% Find all the hdr files inside the directory
hdr=ls2cell(fullfile(fullfile(diro,dirp),'*.hdr'));
catch
% Some checks and balances
disp(sprintf('Looking inside %s I am finding\n',fullfile(diro,dirp)))
ls(fullfile(diro,dirp))
disp('which I expect to contain at least one hdr file')
end
% We know how many header lines there are in each of the hdr files
nhdr=6;
% Now read all the header variables
for index=1:length(hdr)
% The HDR filename
fhdr=fullfile(diro,dirp,hdr{index});
% Read it in
fid=fopen(fhdr);
H=textscan(fid,'%s %d',nhdr);
% Shove the values inside a growing cell array
TV{index}=H{2};
% Be sure to close the files
fclose(fid);
end
% Collate all the header information in TA and keep the names in TN
TN=H{1};
TA=[TV{:}];
% If you want the box corners
if nargout>4
for index=1:length(hdr)
nc=TV{index}(1);
nr=TV{index}(2);
xl=TV{index}(3);
yl=TV{index}(4);
sp=TV{index}(5);
% Plot the outer extent of the boxes, as I interpret it now
bx(index,:)=double([xl xl xl xl xl]+[0 0 nc*sp nc*sp 0]);
by(index,:)=double([yl yl yl yl yl]+[0 nr*sp nr*sp 0 0]);
end
else
[bx,by]=deal(NaN);
end
%%%%%%%%%% VISUAL CHECK %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Make a plot of all the metadata in your directory
if xver==2
% Plot ALL the boxes of the header, they are supposedly all in zone 32
% Compare to http://tinitaly.pi.ingv.it/immagini/Imm_TINITALY_DOWNLOAD_03.jpg
clf
ah=gca;
[BX,BY]=deal(nan(length(hdr),2));
for index=1:length(hdr)
nc=TV{index}(1);
nr=TV{index}(2);
xl=TV{index}(3);
yl=TV{index}(4);
sp=TV{index}(5);
% Plot the outer extent of the boxes, as I interpret it now
bbx=double([xl xl xl xl xl]+[0 0 nc*sp nc*sp 0]);
bby=double([yl yl yl yl yl]+[0 nr*sp nr*sp 0 0]);
plot(bbx,bby); hold on
tt(index)=text(bbx(1)+[bbx(3)-bbx(1)]/2,...
bby(1)+[bby(2)-bby(1)]/2,...
sprintf('%i %s',index,pref(pref(hdr{index}),'_')));
BX(index,:)=minmax(bbx);
BY(index,:)=minmax(bby);
end
hold off
axis image
xel=[min(BX(:,1)) max(BX(:,2))];
yel=[min(BY(:,1)) max(BY(:,2))];
xlim(xel+[-1 1]*range(xel)/20)
ylim(yel+[-1 1]*range(yel)/20)
% Annotate
shrink(ah,1.5,1.5)
tl=title(sprintf('From the headers inside\n %s',...
fullfile(diro,dirp)));
movev(tl,range(ylim)/10)
else
[tt,tl]=deal(NaN);
end
% All the outputs fit to print
varns={hdr,TV,TN,TA,bx,by,tt,tl};
varargout=varns(1:nargout);