Skip to content

Commit 02bbf4e

Browse files
author
Anders Krogh Mortensen
committed
Added script for creating reference areas from ROI
1 parent f968d76 commit 02bbf4e

File tree

1 file changed

+245
-0
lines changed

1 file changed

+245
-0
lines changed

‎Utilities/main_projecting_ROIs.m‎

Lines changed: 245 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,245 @@
1+
clearvars;
2+
close all;
3+
4+
margin = 1; % m. Distance to skip outside ROI
5+
width = 1; % m. Width of new ROI to either of the short ends
6+
7+
%
8+
% width margin margin width
9+
% <-----> <----------> <----------> <----->
10+
%
11+
% +-------+ +-------------------------------------------+ +-------+
12+
% | | | | | |
13+
% | New | | | | New |
14+
% | ROI | | ROI | | ROI |
15+
% | "LEFT"| | (parcel / plot) | |"RIGHT"|
16+
% | | | | | |
17+
% +-------+ +-------------------------------------------+ +-------+
18+
%
19+
%
20+
21+
22+
%%
23+
addpath(genpath('../common'));
24+
25+
%% Spreadsheet selection
26+
27+
[xlsROIfile, xlsROIfolder] = uigetfile({'*.xlsx';'*.xls'},'Select Excel spreadsheet with ROIs.');
28+
xlsFullfile = fullfile(xlsROIfolder, xlsROIfile);
29+
30+
% Sheet selection
31+
disp(['ROI spreadsheet folder : ' xlsROIfolder]);
32+
disp(['ROI spreadsheet : ' xlsROIfile]);
33+
34+
[ sheet, sheets ] = xlsSelectSheet( xlsFullfile );
35+
36+
disp(['ROI sheet : ' sheet]);
37+
38+
%% Read ROIs from XLS sheet
39+
40+
[ ROIs, ROIheader ] = readROIfromXLS( xlsFullfile, sheet );
41+
42+
disp('ROI header:');
43+
disp(ROIheader)
44+
disp('ROI areas:');
45+
disp(ROIheader.Area)
46+
47+
%% Calculate UTM coordinates of ROIs
48+
49+
% Assume all points are in the same UTM zone (if not, they should be close
50+
% enough to the border that the error is sufficiently small)
51+
zone = utmzone(ROIs(1).latitude, ROIs(1).longitude);
52+
[ellipsoid,estr] = utmgeoid(zone);
53+
mstruct = defaultm('utm');
54+
mstruct.zone = zone;
55+
mstruct.geoid = ellipsoid;
56+
mstruct = defaultm(mstruct);
57+
for r = 1:length(ROIs)
58+
ROI = ROIs(r);
59+
[X, Y, Z] = mfwdtran(mstruct, ROI.latitude, ROI.longitude, ROI.altitude);
60+
ROIs(r).E = X;
61+
ROIs(r).N = Y;
62+
ROIs(r).Z = Z;
63+
end
64+
65+
%%
66+
67+
for i = 1:length(ROIs)
68+
ROI = ROIs(i);
69+
70+
X = ROI.E;
71+
Y = ROI.N;
72+
73+
theta = atan2d(Y-mean(Y), X-mean(X));
74+
75+
[~, sort_idx] = sort(theta);
76+
77+
X = X(sort_idx);
78+
Y = Y(sort_idx);
79+
80+
d12 = sqrt((X(1)-X(2))^2 + (Y(1)-Y(2))^2);
81+
d23 = sqrt((X(2)-X(3))^2 + (Y(2)-Y(3))^2);
82+
d34 = sqrt((X(3)-X(4))^2 + (Y(3)-Y(4))^2);
83+
d40 = sqrt((X(4)-X(1))^2 + (Y(4)-Y(1))^2);
84+
85+
if (d12 + d34 > d23 + d40)
86+
longside_idx = [1 2 3 4];
87+
shortside_idx = [2 3 4 1];
88+
else
89+
longside_idx = [2 3 4 1];
90+
shortside_idx = [1 2 3 4];
91+
end
92+
93+
ROIs(i).E = X;
94+
ROIs(i).N = Y;
95+
ROIs(i).Z = ROI.Z(sort_idx);
96+
ROIs(i).longitude = ROI.longitude(sort_idx);
97+
ROIs(i).latitude = ROI.latitude(sort_idx);
98+
ROIs(i).altitude = ROI.altitude(sort_idx);
99+
ROIs(i).longside_idx = longside_idx;
100+
ROIs(i).shortside_idx = shortside_idx;
101+
102+
% ROIs(i) = ROI;
103+
104+
end
105+
106+
%%
107+
108+
ROIs_right_corners = nan(length(ROIs), 4, 2);
109+
ROIs_left_corners = nan(length(ROIs), 4, 2);
110+
111+
for i = 1:length(ROIs)
112+
ROI = ROIs(i);
113+
longside_idx = ROI.longside_idx;
114+
shortside_idx = ROI.shortside_idx;
115+
116+
X = ROI.E;
117+
Y = ROI.N;
118+
119+
X_offset = mean(X);
120+
Y_offset = mean(Y);
121+
122+
X = X-X_offset;
123+
Y = Y-Y_offset;
124+
125+
v1 = [X(longside_idx(2))-X(longside_idx(1)), Y(longside_idx(2))-Y(longside_idx(1))];
126+
u1 = v1/norm(v1);
127+
v2 = [X(longside_idx(3))-X(longside_idx(4)), Y(longside_idx(3))-Y(longside_idx(4))]; % Swap large and small index to ensure, that the two lines point in the same general direction
128+
u2 = v2/norm(v2);
129+
130+
roi_right_corners = [[X(longside_idx(2)), Y(longside_idx(2))]+[u1*margin; u1*(margin+width)]; [X(longside_idx(3)), Y(longside_idx(3))]+[u2*(margin+width); u2*margin]];
131+
roi_left_corners = [[X(longside_idx(1)), Y(longside_idx(1))]-[u1*margin; u1*(margin+width)]; [X(longside_idx(4)), Y(longside_idx(4))]-[u2*(margin+width); u2*margin]];
132+
133+
ROIs_right_corners(i,:,:) = roi_right_corners + [X_offset, Y_offset];
134+
ROIs_left_corners(i,:,:) = roi_left_corners + [X_offset, Y_offset];
135+
136+
end
137+
138+
%% Export to xlsx (RIGHT)
139+
name_parts = strsplit(xlsROIfile,'.');
140+
output_filename = fullfile(xlsROIfolder, [name_parts{1} '__RIGHT.' name_parts{2}]);
141+
142+
A = {'Version',ROIheader.Version;
143+
'Name', [ROIheader.Name '__RIGHT'];
144+
'Description','ROI_RIGHT';
145+
'Location',ROIheader.Location;
146+
'Date',ROIheader.Date;
147+
'Coordinate system','LL';
148+
'Unit','DEG';
149+
'UTM Zone',zone};
150+
xlswrite(output_filename, A, 'A1:B8');
151+
152+
xlswrite(output_filename, {'UTM East, m','UTM North, m','Height, m', 'ROI id'}, 'A10:D10');
153+
xlswrite(output_filename, {'Longitude','Latitude','Elevation, m', 'ROI id'}, 'A10:D10');
154+
155+
% Export data
156+
A = cell(length(ROIs)*4,4);
157+
for i = 1:length(ROIs)
158+
159+
[lat,lon] = minvtran(mstruct, ROIs_right_corners(i,1,1),ROIs_right_corners(i,1,2));
160+
A{(i-1)*4+1,1} = lon;
161+
A{(i-1)*4+1,2} = lat;
162+
A{(i-1)*4+1,3} = NaN;
163+
A{(i-1)*4+1,4} = ROIs(i).name{1};
164+
165+
[lat,lon] = minvtran(mstruct, ROIs_right_corners(i,2,1),ROIs_right_corners(i,2,2));
166+
A{(i-1)*4+2,1} = lon;
167+
A{(i-1)*4+2,2} = lat;
168+
A{(i-1)*4+2,3} = NaN;
169+
A{(i-1)*4+2,4} = ROIs(i).name{1};
170+
171+
[lat,lon] = minvtran(mstruct, ROIs_right_corners(i,3,1),ROIs_right_corners(i,3,2));
172+
A{(i-1)*4+3,1} = lon;
173+
A{(i-1)*4+3,2} = lat;
174+
A{(i-1)*4+3,3} = NaN;
175+
A{(i-1)*4+3,4} = ROIs(i).name{1};
176+
177+
[lat,lon] = minvtran(mstruct, ROIs_right_corners(i,4,1),ROIs_right_corners(i,4,2));
178+
A{(i-1)*4+4,1} = lon;
179+
A{(i-1)*4+4,2} = lat;
180+
A{(i-1)*4+4,3} = NaN;
181+
A{(i-1)*4+4,4} = ROIs(i).name{1};
182+
end
183+
184+
xlswrite(output_filename, A, ['A11:D' num2str(size(A,1)+10)]);
185+
186+
%% Export to xlsx (LEFT)
187+
name_parts = strsplit(xlsROIfile,'.');
188+
output_filename = fullfile(xlsROIfolder, [name_parts{1} '__LEFT.' name_parts{2}]);
189+
190+
A = {'Version',ROIheader.Version;
191+
'Name', [ROIheader.Name '__LEFT'];
192+
'Description','ROI_LEFT';
193+
'Location',ROIheader.Location;
194+
'Date',ROIheader.Date;
195+
'Coordinate system','LL';
196+
'Unit','DEG';
197+
'UTM Zone',zone};
198+
xlswrite(output_filename, A, 'A1:B8');
199+
200+
xlswrite(output_filename, {'UTM East, m','UTM North, m','Height, m', 'ROI id'}, 'A10:D10');
201+
xlswrite(output_filename, {'Longitude','Latitude','Elevation, m', 'ROI id'}, 'A10:D10');
202+
203+
% Export data
204+
A = cell(length(ROIs)*4,4);
205+
for i = 1:length(ROIs)
206+
207+
[lat,lon] = minvtran(mstruct, ROIs_left_corners(i,1,1),ROIs_left_corners(i,1,2));
208+
A{(i-1)*4+1,1} = lon;
209+
A{(i-1)*4+1,2} = lat;
210+
A{(i-1)*4+1,3} = NaN;
211+
A{(i-1)*4+1,4} = ROIs(i).name{1};
212+
213+
[lat,lon] = minvtran(mstruct, ROIs_left_corners(i,2,1),ROIs_left_corners(i,2,2));
214+
A{(i-1)*4+2,1} = lon;
215+
A{(i-1)*4+2,2} = lat;
216+
A{(i-1)*4+2,3} = NaN;
217+
A{(i-1)*4+2,4} = ROIs(i).name{1};
218+
219+
[lat,lon] = minvtran(mstruct, ROIs_left_corners(i,3,1),ROIs_left_corners(i,3,2));
220+
A{(i-1)*4+3,1} = lon;
221+
A{(i-1)*4+3,2} = lat;
222+
A{(i-1)*4+3,3} = NaN;
223+
A{(i-1)*4+3,4} = ROIs(i).name{1};
224+
225+
[lat,lon] = minvtran(mstruct, ROIs_left_corners(i,4,1),ROIs_left_corners(i,4,2));
226+
A{(i-1)*4+4,1} = lon;
227+
A{(i-1)*4+4,2} = lat;
228+
A{(i-1)*4+4,3} = NaN;
229+
A{(i-1)*4+4,4} = ROIs(i).name{1};
230+
end
231+
232+
xlswrite(output_filename, A, ['A11:D' num2str(size(A,1)+10)]);
233+
234+
%%
235+
236+
figure;
237+
hold on;
238+
for i = 1:length(ROIs)
239+
h1 = plot(ROIs(i).E, ROIs(i).N,'color',[0 1 0]);
240+
end
241+
h2 = plot(ROIs_right_corners(:,:,1)', ROIs_right_corners(:,:,2)','color',[1 0 0]);
242+
h3 = plot(ROIs_left_corners(:,:,1)', ROIs_left_corners(:,:,2)','color',[0 0 1]);
243+
legend([h1(1), h2(1), h3(1)],{'Plots','"Right"','"Left"'},'Location','EastOutside')
244+
245+
axis equal;

0 commit comments

Comments
 (0)