-
Notifications
You must be signed in to change notification settings - Fork 110
/
getoutline.m
135 lines (117 loc) · 2.86 KB
/
getoutline.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
function varargout = getoutline(DEM,nnan,wm)
%GETOUTLINE get or plot extent of GRIDobj
%
% Syntax
%
% getoutline(DEM)
% MS = getoutline(DEM)
% [x,y] = getoutline(DEM)
% ... = getoutline(DEM,removenans)
% getoutline(DEM,removenans,wm)
%
% Description
%
% getoutline plots the extent of a GRIDobj or returns the coordinate
% vectors that generate the plot. By default, getoutline returns the
% coordinate vectors of the DEM edges. By setting removenans = true,
% you can get the outline around the valid (non-nan) data in the DEM.
%
% Input arguments
%
% DEM GRIDobj
% removenans set to true, if the outline should be around the valid
% data in the DEM.
% wm plot in webmap (true or {false})
%
% Output arguments
%
% MS mapping structure that can be displayed using mapshow or
% exported with shapewrite.
% x,y coordinate vectors that can be used to plot the extent
% rectangle (plot(x,y))
%
% Example
%
% DEM = GRIDobj('srtm_bigtujunga30m_utm11.tif');
% getoutline(DEM,true,true)
%
% Author: Wolfgang Schwanghart (w.schwanghart[at]geo.uni-potsdam.de)
% Date: 17. October, 2018
nargoutchk(0,2);
if nargin >= 2
if nnan
I = isnan(DEM);
if ~any(I)
nnan = false;
end
end
else
nnan = false;
end
if ~nnan
[x,y] = getcoordinates(DEM);
maxx = max(x);
minx = min(x);
maxy = max(y);
miny = min(y);
x = [minx minx maxx maxx minx];
y = [miny maxy maxy miny miny];
else
I = ~I;
B = bwboundaries(I.Z);
x = [];
y = [];
for k = 1:length(B)
boundary = B{k};
[xb,yb] = sub2coord(I,boundary(:,1),boundary(:,2));
x = [x;xb;nan];
y = [y;yb;nan];
end
end
% simplify lines
xy = dpsimplify([x(:) y(:)],100*eps);
x = xy(:,1);
y = xy(:,2);
if nargout == 0
% No output. Plot outline
if nargin <= 2
wm = false;
end
if ~wm
plot(x,y,'LineWidth',3,'Color',[.8 .8 .8]);
hh = ishold;
hold on
plot(x,y,'k--','LineWidth',1);
if ~hh
hold off;
end
else
[lat,lon] = minvtran(DEM.georef.mstruct,x,y);
wmline(lat,lon)
end
elseif nargout == 1
% One output argument, create mapping structure
% split at nans
if ~isnan(x(end))
x = [x;nan];
y = [y;nan];
end
I = isnan(x);
nrlines = nnz(I);
ix = find(I)';
ixs = [1; ix(1:end-1)'+1];
ixe = ix-1;
for r = 1:nrlines
MS(r).Geometry = 'Polygon';
MS(r).X = x(ixs(r):ixe(r));
MS(r).Y = y(ixs(r):ixe(r));
MS(r).ID = r;
MS(r).name = DEM.name;
end
varargout{1} = MS;
elseif nargout == 2
% Two outputs, return x and y coordinates
varargout{1} = x;
varargout{2} = y;
end
end