GIS算法——点是否在指定区域内

GIS算法——点是否在指定区域内

代码:
GPoint = record
X: Double; //经度，也使用于任何座标
Y: Double; //纬度
end;


GRange = array of GPoint;

函数返回true 点为在区域内


function GPointInPolygonX(P: GPoint; L: GRange): Boolean;
var
Num, i: Integer;
a, da, suma, pa, a0: double;
dp0, dp1: double;
twopi, pi: double;
begin
Num := Length(L);
pi := 3.1415926535897932384626433832795;
twopi := 2.0 * pi;

suma := 0.0;
dp0 := L[0].x - p.x;
dp1 := L[0].y - p.y;
a0 := ArcTan2(dp1, dp0);
if a0 < 0.0 then
a0 := a0 + twopi;

pa := a0;
for i := 1 to Num - 1 do
begin
dp0 := L[i].x - p.x;
dp1 := L[i].y - p.y;
a := ArcTan2(dp1, dp0);
if a < 0.0 then
a := a + twopi;
da := a - pa;
if da < -pi then
da := da + twopi; //insure the direction is croccet
if da > pi then
da := da - twopi;
suma := suma + da;
pa := a;
end;

da := a0 - pa;
if da < -pi then
da := da + twopi;
if da > pi then
da := da - twopi;
suma := suma + da;
//return (int)((fabs(suma)+1.0e-5)/twopi);
if floor((abs(suma) + 1.0E-5) / twopi) = 0 then
Result := False
else
Result := True;
end;