Description of geo

0001 function [olon, olat] = geo_simplify(lon, lat, nmax, tolerance)
0002 % GEO_SIMPLIFY -- Simplify (reduce the points of) a lat/lon polyline
0003 %  [olon, olat] = geo_simplify(lon, lat, nmax, tolerance)
0004 %   simplifies the (lon, lat) contour, returning nmax points
0005 %   if (nmax > 0) or those points that are more significant
0006 %   than the given tolerance.  The gs_sort() routine is used,
0007 %   in which the deviation, not the curvature, is the key
0008 %   principle.  NaNs are ignored.
0009 %  geo_simplity(...) returns [olon olat] in one array.
0010 %  geo_simplify(nPoints) demonstrates itself with nPoints
0011 %   (random); default = 24.
0012  
0013 % Copyright (C) 1997 Dr. Charles R. Denham, ZYDECO.
0014 %  All Rights Reserved.
0015 %   Disclosure without explicit written consent from the
0016 %    copyright owner does not constitute publication.
0017  
0018 % Version of 06-Jul-1998 06:40:26.
0019 % Revised    04-Sep-1998 15:29:19.
0020 
0021 % Reference: Douglas, D.H. and T.K. Peucker, Algorithms for
0022 %  the reduction of the number of points required to represent
0023 %  a digitized line of (sic?) its caricature, Can. Cartogr.,
0024 %  10, 112-122, 1973.  We have adopted a similar scheme for
0025 %  selecting "significant" points.  Unlike the D-P algorithm,
0026 %  we sort all the points first, then isolate those that satisfy
0027 %  our criteria for tolerance or total number of points retained.
0028 
0029 if nargin < 1, help(mfilename), lon = 'demo'; end
0030 
0031 if isequal(lon, 'demo'), lon = 24; end
0032 
0033 if ischar(lon), lon = eval(lon); end
0034 
0035 if length(lon) == 1
0036         n = lon;
0037         lon = 1:n;
0038         lat = cumsum(rand(size(lon))-0.5) .* 10;
0039         lat(ceil(n/2)) = NaN;
0040         nmax = fix((n+2) ./ 3);
0041         tolerance = 0;
0042         tic
0043         [lon1, lat1] = geo_simplify(lon, lat, nmax, tolerance);
0044         toc
0045         plot(lon1, lat1, '-*', lon, lat, 'o')
0046         xlabel('lon (degrees)'), ylabel('lat (degrees)')
0047         title('geo\_simplify')
0048         axis equal
0049         set(gcf, 'Name', 'geo_simplify')
0050         figure(gcf)
0051         return
0052 end
0053 
0054 if nargin < 3, nmax = length(lon); end
0055 if nargin < 4
0056         tolerance = 0;
0057         if nmax == 0, nmax = length(lon); end
0058 end
0059 
0060 if nmax > 0 & nmax <= 1
0061         nmax = round(nmax*length(lon));
0062 end
0063 
0064 RCF = 180 ./ pi;
0065 x = cos(lat/RCF).*cos(lon/RCF);
0066 y = cos(lat/RCF).*sin(lon/RCF);
0067 z = sin(lat/RCF);
0068 
0069 % Save NaNs.
0070 
0071 xx = x(:);
0072 yy = y(:);
0073 zz = z(:);
0074 
0075 i = find(isfinite(xx) & isfinite(yy) & isfinite(zz));
0076 j = find(~isfinite(xx) | ~isfinite(yy) | ~isfinite(zz));
0077 if any(i), xx = xx(i); yy = yy(i); zz = zz(i); end
0078 
0079 % Sort.
0080 
0081 [ind, dev] = gs_sort(xx, yy, zz);
0082 
0083 % Restore NaNs.
0084 
0085 ind = i(ind);
0086 ind = [ind; j(:)];
0087 dev = [dev; NaN*j];
0088 
0089 % Deviation angle: always < 90 degrees.
0090 
0091 dev = asin(dev) .* 180 ./ pi;
0092 
0093 % Trim by count or tolerance.
0094 
0095 if nmax > 0
0096         if nmax <= 1, nmax = round(nmax .* length(lon)); end
0097         nmax = max(3, min(nmax, length(lon)));
0098         ind = ind(1:nmax);
0099         dev = dev(1:nmax);
0100 else
0101         f = find(dev >= tolerance);
0102         if any(f)
0103                 ind = ind(f);
0104                 dev = dev(f);
0105         end
0106 end
0107 
0108 % Unsort.
0109 
0110 [ind, s] = sort(ind);
0111 dev = dev(s);
0112 
0113 lon1 = lon(ind);
0114 lat1 = lat(ind);
0115 
0116 if nargout == 1
0117         olon = [lon1 lat1];
0118 elseif nargout > 1
0119         olon = lon1;
0120         olat = lat1;
0121 end
0122 
0123 function [theIndices, theDeviations] = gs_sort(x, y, z)
0124 
0125 % gs_sort -- Sort (x, y, z) data for line-simplification.
0126 %  [theIndices, theDeviations] = gs_sort(x, y, z) returns
0127 %   indices and deviations of the given (x, y, z) data,
0128 %   in order of decreasing importance.  The (x, y, z)
0129 %   data are geographic, expressed as direction-cosines.
0130 %   The deviations are the sines of the angles of
0131 %   deviation.
0132  
0133 % Copyright (C) 1997 Dr. Charles R. Denham, ZYDECO.
0134 %  All Rights Reserved.
0135 %   Disclosure without explicit written consent from the
0136 %    copyright owner does not constitute publication.
0137  
0138 % Version of 02-Jul-1998 09:36:15.
0139 
0140 if nargin < 1, help(mfilename), x = 'demo'; end
0141 
0142 if isequal(x, 'demo'), x = 24; end
0143 
0144 if ischar(x), x = eval(x); end
0145 
0146 if length(x) == 1
0147         n = x;
0148         a = cumsum(rand(n, 3)-0.5);
0149         for i = 1:n
0150                 a(i, :) = a(i, :) ./ norm(a(i, :));
0151         end
0152         x = a(:, 1); y = a(:, 2); z = a(:, 3);
0153         tic
0154         [i, d] = gs_sort(x, y, z);
0155         toc
0156         reduce = 3;   % Reduction-factor.
0157         k = fix((n+reduce-1)./reduce);
0158         ii = i(1:k);
0159         dd = d(1:k);
0160         xx = x(ii);
0161         yy = y(ii);
0162         zz = z(ii);
0163         [ii, s] = sort(ii);
0164         dd = dd(s);
0165         xx = xx(s);
0166         yy = yy(s);
0167         zz = zz(s);
0168         subplot(2, 1, 1)
0169         plot3(x, y, z, 'b+', xx, yy, zz, 'go-')
0170         title('Data and Subset')
0171         xlabel('x'), ylabel('y'), zlabel('z')
0172         view(2)
0173         subplot(2, 1, 2)
0174         plot(i, d, '.', xx, dd, 'o')
0175         title('Deviation')
0176         xlabel('x'), ylabel('y'), zlabel('z')
0177         view(2)
0178         set(gcf, 'Name', 'gs_sort demo')
0179         figure(gcf)
0180         return
0181 end
0182 
0183 ind = zeros(size(x));
0184 dev = zeros(size(x));
0185 
0186 nstack = ceil(1 + 2*log(length(x))/log(2));
0187 stack = zeros(nstack, 2);
0188 nstack = 1;
0189 stack(nstack, :) = [1 length(x)];
0190 
0191 ind(1) = 1;
0192 ind(2) = length(x);
0193 nout = 2;
0194 
0195 % Call gs_pivot() repeatedly and stack the
0196 %  longer interval before the shorter.
0197 
0198 while nstack > 0
0199         a = stack(nstack, 1);
0200         b = stack(nstack, 2);
0201         nstack = nstack - 1;
0202         if b-a > 1
0203                 [p, d] = gs_pivot(x([a:b]), y([a:b]), z([a:b]));
0204                 if any(p)
0205                         p = a+p-1;
0206                         nout = nout+1;
0207                         ind(nout) = p;
0208                         dev(nout) = d;
0209                         if b-p >= p-a
0210                                 if b-p > 1
0211                                         nstack = nstack + 1;
0212                                         stack(nstack, :) = [p b];
0213                                 end
0214                         end
0215                         if p-a > 1
0216                                 nstack = nstack + 1;
0217                                 stack(nstack, :) = [a p];
0218                         end
0219                         if b-p < p-a
0220                                 if b-p > 1
0221                                         nstack = nstack + 1;
0222                                         stack(nstack, :) = [p b];
0223                                 end
0224                         end
0225                 end
0226         end
0227 end
0228 
0229 % Sort in descending order of importance.
0230 
0231 dev(1:2) = max(dev);
0232 [ignore, s] = sort(-dev);
0233 ind = ind(s);
0234 dev = dev(s);
0235 
0236 % Output.
0237 
0238 if nargout > 0
0239         theIndices = ind;
0240         theDeviations = dev;
0241 else
0242         ind, dev
0243 end
0244 
0245 function [thePivot, theDeviation] = gs_pivot(x, y, z)
0246 
0247 % gs_pivot -- Geographic line-simplification pivot.
0248 %  [thePivot, theDeviation] = gs_pivot(x, y, z) returns
0249 %   the index of the (x, y, z) point that deviates most
0250 %   from the great-circle path between the first and last
0251 %   of the geographic points, given as direction-cosines.
0252 %   Only finite points are allowed.  The deviation is
0253 %   returned as the sine of the deviation-angle.
0254  
0255 % Copyright (C) 1997 Dr. Charles R. Denham, ZYDECO.
0256 %  All Rights Reserved.
0257 %   Disclosure without explicit written consent from the
0258 %    copyright owner does not constitute publication.
0259  
0260 % Version of 02-Jul-1998 08:32:40.
0261 
0262 if nargin < 1, help(mfilename), x = 'demo'; end
0263 
0264 if isequal(x, 'demo'), x = 10; end
0265 
0266 if ischar(x), x = eval(x); end
0267 
0268 if length(x) == 1
0269         n = x;
0270         a = cumsum(rand(n, 3)-0.5);
0271         for i = 1:n
0272                 a(i, :) = a(i, :) ./ norm(a(i, :));
0273         end
0274         x = a(:, 1); y = a(:, 2); z = a(:, 3);
0275         tic
0276         [piv, dev] = gs_pivot(x, y, z)
0277         toc
0278         plot3(x, y, z, 'b-', ...
0279                         x([1 n]), y([1 n]), z([1, n]), 'g-', ...
0280                         x(piv), y(piv), z(piv), 'r*')
0281         title(['gs\_pivot(' int2str(n) ')'])
0282         xlabel('x'), ylabel('y'), zlabel('z')
0283         set(gca, 'Box', 'on')
0284         figure(gcf)
0285         return
0286 end
0287 
0288 thePivot = 0;
0289 theDeviation = -inf;
0290 
0291 n = length(x);
0292 if n < 3, return, end
0293 
0294 % Pole of the spanning great-circle.
0295 %  WARNING: no such pole will exist if the two
0296 %  end-points are the same.  In that case, use
0297 %  point farthest from the ends.
0298 
0299 if any([x(1) y(1) z(1)] ~= [x(n) y(n) z(n)])
0300         c(1) = y(1)*z(n)-y(n)*z(1);
0301         c(2) = z(1)*x(n)-z(n)*x(1);
0302         c(3) = x(1)*y(n)-x(n)*y(1);
0303         c = c.' ./ norm(c);
0304         d = [x y z] * c;
0305 else
0306         c = [x(1) y(1) z(1)].';
0307         d = [x y z] * c;
0308         d = sqrt(1 - d.*d);
0309 end
0310 
0311 % Find the pivot and deviation.
0312 
0313 temp = abs(d(2:n-1));
0314 f = find(temp == max(temp));
0315 
0316 piv = f(1) + 1;
0317 dev = abs(d(piv));
0318 
0319 % Output.
0320 
0321 if nargout > 0
0322         thePivot = piv;
0323         theDeviation = dev;
0324 else
0325         piv, dev
0326 end
geo_simplify

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SUBFUNCTIONS

SOURCE CODE
