-
Notifications
You must be signed in to change notification settings - Fork 110
/
coord2sub.m
executable file
·84 lines (70 loc) · 1.6 KB
/
coord2sub.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
function [r,c,res] = coord2sub(DEM,x,y)
%COORD2SUB convert x and y coordinates to subscripts into a GRIDobj
%
% Syntax
%
% [r,c] = coord2sub(DEM,x,y)
% [r,c,res] = ...
%
% Description
%
% coord2sub converts vectors of x and y coordinates to subscripts
% into an instance of GRIDobj.
%
% Input arguments
%
% DEM instance of GRIDobj
% x,y x- and y-coordinates
%
% Output arguments
%
% r,c row and column subscripts
% res residual distance between coordinates and nearest grid cell
% centers for coordinate pair x, y
%
% See also: GRIDobj/coord2ind, GRIDobj/sub2coord, GRIDobj/getcoordinates
%
% Author: Wolfgang Schwanghart (w.schwanghart[at]geo.uni-potsdam.de)
% Date: 17. August, 2017
narginchk(3,3);
% get coordinate vectors
[X,Y] = refmat2XY(DEM.refmat,DEM.size);
X = X(:);
Y = Y(:);
% check input
np = numel(x);
if np ~= numel(y);
error('TopoToolbox:wronginput',...
'x and y must have the same number of elements.');
end
% force column vectors
x = x(:);
y = y(:);
dx = X(2)-X(1);
dy = Y(2)-Y(1);
IX1 = (x-X(1))./dx + 1;
IX2 = (y-Y(1))./dy + 1;
IX1 = round(IX1);
IX2 = round(IX2);
I = IX1>DEM.size(2) | IX1<1 | IX2>DEM.size(1) | IX2<1;
if any(I(:));
warning('TopoToolbox:outsidegrid',...
'There are some points outside the grid''s borders');
end
x(I) = [];
y(I) = [];
I = ~I;
if any(I)
r = nan(np,1);
r(I) = IX2(I);
c = nan(np,1);
c(I) = IX1(I);
if nargout == 3
res = nan(np,1);
res(I) = hypot(X(IX1(I))-x,Y(IX2(I))-y);
end
else
r = nan(np,1);
c = nan(np,1);
res = nan(np,1);
end