Skip to content

Commit 790030f

Browse files
author
Anders Krogh Mortensen
committed
First public commit
1 parent 92dda21 commit 790030f

File tree

172 files changed

+28025
-0
lines changed

Some content is hidden

Large Commits have some content hidden by default. Use the searchbox below for content that may be hidden.

172 files changed

+28025
-0
lines changed

‎Utilities/main_cropMap.m‎

Lines changed: 162 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,162 @@
1+
clearvars;
2+
close all;
3+
4+
%%
5+
addpath(genpath('../common/'));
6+
7+
%% Select map
8+
9+
disp('Select geotiff map...');
10+
[fileName, pathName] = uigetfile({'*.tif','*.tiff'},'Select geotiff map');
11+
fileNamePath = fullfile(pathName, fileName);
12+
disp('Map selected:');
13+
disp([' File name: ' fileName]);
14+
disp([' Folder: ' pathName]);
15+
16+
%% Select crop polygon file
17+
18+
[polygonFileName, polygonPathName] = uigetfile(fullfile(pathName,'*.mat'),'Select mat-file with polygon');
19+
20+
if all(polygonFileName ~= 0) && all(polygonPathName ~= 0)
21+
polygon = load(fullfile(polygonPathName,polygonFileName));
22+
if isfield(polygon,'X') && isfield(polygon,'Y')
23+
X = polygon.X;
24+
Y = polygon.Y;
25+
disp(['Loaded polygon with ' num2str(length(X)) ' coordinates.']);
26+
else
27+
disp('Mat-file does not contain X- and Y-coordinates...');
28+
end
29+
end
30+
31+
%% Load map
32+
33+
fWaitbar = waitbar(0.25,'Loading map...');
34+
[mapRaster, mapRasterInfo] = geotiffread(fileNamePath);
35+
mapInfo = geotiffinfo(fileNamePath);
36+
mapMaxChanIdx = 1;
37+
if (size(mapRaster,3) >= 3)
38+
mapMaxChanIdx = 3;
39+
end
40+
41+
%%
42+
if (~exist('fWaitbar','var'))
43+
fWaitbar = waitbar(0,'','WindowStyle','modal');
44+
else
45+
if (~isvalid(fWaitbar))
46+
fWaitbar = waitbar(0,'','WindowStyle','modal');
47+
end
48+
end
49+
waitbar(0.5, fWaitbar, 'Copying map to outputmap...','WindowStyle','modal');
50+
mapRasterCropped = mapRaster;
51+
mapRasterCroppedInfo = mapRasterInfo;
52+
53+
waitbar(0.75, fWaitbar, 'Plot maps...','WindowStyle','modal');
54+
fMaps = figure('Name','Crop map','WindowState','maximized');
55+
ax1 = subplot(1,2,1);
56+
hMap1 = mapshow(im2double(mapRaster(:,:,1:mapMaxChanIdx)), mapRasterInfo);
57+
if (size(mapRaster,3) > mapMaxChanIdx)
58+
set(hMap1,'AlphaData',mapRaster(:,:,mapMaxChanIdx+1))
59+
end
60+
title('Input map');
61+
xlabel([upper(mapRasterInfo.RowsStartFrom(1)) lower(mapRasterInfo.RowsStartFrom(2:end))])
62+
ylabel([upper(mapRasterInfo.ColumnsStartFrom(1)) lower(mapRasterInfo.ColumnsStartFrom(2:end))])
63+
hold on;
64+
ax2 = subplot(1,2,2);
65+
hMap2 = mapshow(im2double(mapRasterCropped(:,:,1:mapMaxChanIdx)), mapRasterCroppedInfo);
66+
if (size(mapRasterCropped,3) > mapMaxChanIdx)
67+
set(hMap2,'AlphaData',mapRasterCropped(:,:,mapMaxChanIdx+1))
68+
end
69+
title('Output map');
70+
xlabel([upper(mapRasterCroppedInfo.RowsStartFrom(1)) lower(mapRasterCroppedInfo.RowsStartFrom(2:end))])
71+
ylabel([upper(mapRasterCroppedInfo.ColumnsStartFrom(1)) lower(mapRasterCroppedInfo.ColumnsStartFrom(2:end))])
72+
hold on;
73+
drawnow;
74+
75+
waitbar(1, fWaitbar, 'Ready for user interaction','WindowStyle','modal');
76+
pause(0.25);
77+
close(fWaitbar);
78+
79+
button = 0;
80+
updateLines = true;
81+
while(~isempty(button))
82+
if (updateLines)
83+
updateLines = false;
84+
ch = get(ax1,'Children');
85+
for c = 1:length(ch)
86+
if isa(ch(c),'matlab.graphics.chart.primitive.Line')
87+
delete(ch(c));
88+
end
89+
end
90+
if (exist('X','var'))
91+
if(length(X) > 0)
92+
plot(ax1,[X X(1)],[Y Y(1)],'-o','Color',[1 1 1],'LineWidth',2);
93+
plot(ax1,[X X(1)],[Y Y(1)],'-o','Color',[1 0 0],'LineWidth',1);
94+
end
95+
end
96+
end
97+
98+
[x,y,button] = ginput(1);
99+
thisAxes = gca;
100+
101+
if (button == 1) % Left mouse button
102+
% Add point to polygon
103+
if (thisAxes == ax1)
104+
if (~exist('X','var'))
105+
X = x;
106+
Y = y;
107+
else
108+
X(end+1) = x;
109+
Y(end+1) = y;
110+
end
111+
updateLines = true;
112+
end
113+
elseif (button == 2) % Shift+mouse button or middle mouse button
114+
% Move neasest point on polygon
115+
if (thisAxes == ax1)
116+
if (length(X) > 0)
117+
[minDist,minDistIdx] = min(pdist2([X' Y'],[x y]));
118+
X(minDistIdx) = x;
119+
Y(minDistIdx) = y;
120+
updateLines = true;
121+
end
122+
end
123+
elseif (button == 3) % Right mouse button
124+
% Remove nearest point in polygon
125+
if (thisAxes == ax1)
126+
if (length(X) > 0)
127+
[minDist,minDistIdx] = min(pdist2([X' Y'],[x y]));
128+
X(minDistIdx) = [];
129+
Y(minDistIdx) = [];
130+
updateLines = true;
131+
end
132+
end
133+
elseif (button == 32) % Space button
134+
% Update map crop
135+
if (length(X) > 2)
136+
fWaitbar = waitbar(0.5,'Cropping map. Please wait...');
137+
[ mapRasterCropped, mapRasterCroppedMask, mapRasterCroppedInfo, ~] = mapcrop( mapRaster, mapRasterInfo, X, Y, 'native' );
138+
cla(ax2);
139+
hMap2 = mapshow(ax2, im2double(mapRasterCropped(:,:,1:mapMaxChanIdx)), mapRasterCroppedInfo);
140+
if (size(mapRasterCropped,3) > mapMaxChanIdx)
141+
set(hMap2,'AlphaData',mapRasterCroppedMask)
142+
end
143+
close(fWaitbar);
144+
end
145+
end
146+
disp(num2str([x,y,button],' %i'));
147+
end
148+
close(fMaps)
149+
150+
fWaitbar = waitbar(0.5,'Saving cropped map...');
151+
[path,filename,ext] = fileparts(fileNamePath);
152+
outputFilePath = fullfile(path, [filename '_cropped' ext]);
153+
% geotiffwrite(outputFilePath,mapRasterCropped,mapRasterCroppedInfo);
154+
GeoKeyDirectoryTag = mapInfo.GeoTIFFTags.GeoKeyDirectoryTag;
155+
saveGeotiffWithMask( outputFilePath, mapRasterCropped, mapRasterCroppedMask, mapRasterCroppedInfo, GeoKeyDirectoryTag,'saveAsPseudoColor',false);
156+
157+
waitbar(0.75,fWaitbar,'Saving polygon...');
158+
save(fullfile(path, [filename '_polygon.mat']),'X','Y');
159+
160+
close(fWaitbar);
161+
162+
disp('Done');

‎Utilities/main_showROIsOnMap.m‎

Lines changed: 85 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,85 @@
1+
clearvars
2+
% close all
3+
4+
ROIlinewidth = 1; % Line width (in px) used when plotting ROIs. Increase if, the are hard to see.
5+
6+
%%
7+
addpath(genpath('../common'));
8+
9+
%% Select map
10+
11+
[mapFile, mapPath] = uigetfile('*.tif','Select map (geotif)');
12+
mapFilePath = fullfile(mapPath, mapFile);
13+
14+
%% Spreadsheet selection
15+
16+
[xlsROIfile, xlsROIfolder] = uigetfile({'*.xlsx';'*.xls'},'Select Excel spreadsheet with ROIs.');
17+
xlsFullfile = fullfile(xlsROIfolder, xlsROIfile);
18+
19+
% Sheet selection
20+
disp(['ROI spreadsheet folder : ' xlsROIfolder]);
21+
disp(['ROI spreadsheet : ' xlsROIfile]);
22+
23+
[ sheet, sheets ] = xlsSelectSheet( xlsFullfile );
24+
25+
disp(['ROI sheet : ' sheet]);
26+
27+
%% Read ROIs from XLS sheet
28+
29+
[ ROIs, ROIheader ] = readROIfromXLS( xlsFullfile, sheet );
30+
31+
disp('ROI header:');
32+
disp(ROIheader)
33+
disp('ROI areas:');
34+
disp(ROIheader.Area)
35+
36+
%% Load geotif
37+
38+
disp('Loading geotiff...');
39+
[map, mapRasterInfo] = geotiffread(mapFilePath);
40+
[mapInfo] = geotiffinfo(mapFilePath);
41+
if (isa(map,'single'))
42+
map = double(map);
43+
end
44+
if (size(map,3) == 1)
45+
chns = 1;
46+
elseif (size(map,3) == 2)
47+
chns = 1;
48+
elseif (size(map,3) == 3)
49+
chns = 1:3;
50+
elseif (size(map,3) == 4)
51+
chns = 1:3;
52+
else
53+
chns = 1;
54+
end
55+
disp('Loading geotiff completed!');
56+
57+
%% Convert ROIs to UTM coordinates
58+
59+
mstruct = geotiff2mstruct(mapInfo);
60+
for r = 1:length(ROIs)
61+
ROI = ROIs(r);
62+
[X, Y, Z] = mfwdtran(mstruct, ROI.latitude, ROI.longitude, ROI.altitude);
63+
ROIs(r).X = X;
64+
ROIs(r).Y = Y;
65+
ROIs(r).Z = Z;
66+
end
67+
68+
%% Plot
69+
disp('Displaying map...');
70+
figure;
71+
mapshow(map(:,:,chns), mapRasterInfo)
72+
xlabel('UTM East, m');
73+
ylabel('UTM North, m');
74+
title(['Map: ' mapFile],'Interpreter','none');
75+
disp('Plotting ROIs on map...');
76+
hold on;
77+
cmap = lines(length(ROIs));
78+
for r = 1:length(ROIs)
79+
X = ROIs(r).X;
80+
Y = ROIs(r).Y;
81+
plot([X X(1)],[Y Y(1)],'-','color',[1 1 1],'LineWidth',ROIlinewidth+1);
82+
plot([X X(1)],[Y Y(1)],'-','color',cmap(r,:),'LineWidth',ROIlinewidth);
83+
text(mean(X),mean(Y),ROIs(r).name,'VerticalAlignment','middle','HorizontalAlignment','center','Background',cmap(r,:));
84+
end
85+
disp('Done!');

0 commit comments

Comments
 (0)