1 #include "point.h"
2
3 /* This file contains definitions for spherical point functions. */
4
5 PG_FUNCTION_INFO_V1(spherepoint_in);
6 PG_FUNCTION_INFO_V1(spherepoint_from_long_lat);
7 PG_FUNCTION_INFO_V1(spherepoint_distance);
8 PG_FUNCTION_INFO_V1(spherepoint_long);
9 PG_FUNCTION_INFO_V1(spherepoint_lat);
10 PG_FUNCTION_INFO_V1(spherepoint_x);
11 PG_FUNCTION_INFO_V1(spherepoint_y);
12 PG_FUNCTION_INFO_V1(spherepoint_z);
13 PG_FUNCTION_INFO_V1(spherepoint_xyz);
14 PG_FUNCTION_INFO_V1(spherepoint_equal);
15
16 bool
spoint_eq(const SPoint * p1,const SPoint * p2)17 spoint_eq(const SPoint *p1, const SPoint *p2)
18 {
19 Vector3D a,
20 b;
21
22 spoint_vector3d(&a, p1);
23 spoint_vector3d(&b, p2);
24 return (vector3d_eq(&a, &b));
25 }
26
27 void
spoint_check(SPoint * spoint)28 spoint_check(SPoint *spoint)
29 {
30 bool lat_is_neg = (spoint->lat < 0) ? true : false;
31
32 if (spoint->lng < 0 || spoint->lng > PID)
33 spoint->lng = spoint->lng - floor(spoint->lng / (PID)) * PID;
34 if (spoint->lat < -PIH || spoint->lat > PIH)
35 spoint->lat = spoint->lat - floor(spoint->lat / (PID)) * PID;
36 if (spoint->lng < 0.0)
37 {
38 spoint->lng += (PID);
39 }
40 if (spoint->lat > PI)
41 {
42 spoint->lat -= (2 * PI);
43 }
44 if (spoint->lat > PIH)
45 {
46 spoint->lat = (PI - spoint->lat);
47 spoint->lng += ((spoint->lng < PI) ? (PI) : (-PI));
48 }
49 if (spoint->lat < -PIH)
50 {
51 spoint->lat = (-PI - spoint->lat);
52 spoint->lng += ((spoint->lng < PI) ? (PI) : (-PI));
53 }
54 if (FPeq(spoint->lat, PIH) && lat_is_neg)
55 spoint->lat = -PIH;
56
57 if (FPeq(spoint->lng, PID))
58 {
59 spoint->lng = 0.0;
60 }
61 if (FPzero(spoint->lng))
62 {
63 spoint->lng = 0.0;
64 }
65 if (FPzero(spoint->lat))
66 {
67 spoint->lat = 0.0;
68 }
69 }
70
71 void
vector3d_spoint(SPoint * p,const Vector3D * v)72 vector3d_spoint(SPoint *p, const Vector3D *v)
73 {
74 double rho = sqrt((v->x) * (v->x) + (v->y) * (v->y));
75
76 if (0.0 == rho)
77 {
78 if (FPzero(v->z))
79 {
80 p->lat = 0.0;
81 }
82 else if (v->z > 0)
83 {
84 p->lat = PIH;
85 }
86 else if (v->z < 0)
87 {
88 p->lat = -PIH;
89 }
90 }
91 else
92 {
93 p->lat = atan(v->z / rho);
94 }
95
96 p->lng = atan2(v->y, v->x);
97 if (FPzero(p->lng))
98 {
99 p->lng = 0.0;
100 }
101 else if (p->lng < 0.0)
102 {
103 p->lng += PID;
104 }
105 }
106
107 void
spoint_vector3d(Vector3D * v,const SPoint * p)108 spoint_vector3d(Vector3D *v, const SPoint *p)
109 {
110 v->x = cos(p->lng) * cos(p->lat);
111 v->y = sin(p->lng) * cos(p->lat);
112 v->z = sin(p->lat);
113 }
114
115 Datum
spherepoint_in(PG_FUNCTION_ARGS)116 spherepoint_in(PG_FUNCTION_ARGS)
117 {
118 SPoint *sp = (SPoint *) palloc(sizeof(SPoint));
119 char *c = PG_GETARG_CSTRING(0);
120 double lng,
121 lat;
122
123 void sphere_yyparse(void);
124
125 init_buffer(c);
126 sphere_yyparse();
127 if (get_point(&lng, &lat))
128 {
129 sp->lng = lng;
130 sp->lat = lat;
131 spoint_check(sp);
132 }
133 else
134 {
135 reset_buffer();
136 pfree(sp);
137 sp = NULL;
138 elog(ERROR, "spherepoint_in: parse error");
139 }
140 reset_buffer();
141 PG_RETURN_POINTER(sp);
142 }
143
144
145 Datum
spherepoint_from_long_lat(PG_FUNCTION_ARGS)146 spherepoint_from_long_lat(PG_FUNCTION_ARGS)
147 {
148 SPoint *p = (SPoint *) palloc(sizeof(SPoint));
149
150 p->lng = PG_GETARG_FLOAT8(0);
151 p->lat = PG_GETARG_FLOAT8(1);
152 spoint_check(p);
153 PG_RETURN_POINTER(p);
154 }
155
156 static double
norm2(double a,double b)157 norm2(double a, double b)
158 {
159 return sqrt(a * a + b * b);
160 }
161
162 float8
spoint_dist(const SPoint * p1,const SPoint * p2)163 spoint_dist(const SPoint *p1, const SPoint *p2)
164 {
165 float8 dl = p1->lng - p2->lng;
166 /* use Vincenty's formula for the inverse geodesic problem on the sphere */
167 float8 f = atan2(norm2(cos(p2->lat) * sin(dl),
168 cos(p1->lat) * sin(p2->lat)
169 - sin(p1->lat) * cos(p2->lat) * cos(dl)),
170 sin(p1->lat) * sin(p2->lat)
171 + cos(p1->lat) * cos(p2->lat) * cos(dl));
172 if (FPzero(f))
173 {
174 return 0.0;
175 }
176 else
177 {
178 return f;
179 }
180 }
181
182 Datum
spherepoint_distance(PG_FUNCTION_ARGS)183 spherepoint_distance(PG_FUNCTION_ARGS)
184 {
185 SPoint *p1 = (SPoint *) PG_GETARG_POINTER(0);
186 SPoint *p2 = (SPoint *) PG_GETARG_POINTER(1);
187
188 PG_RETURN_FLOAT8(spoint_dist(p1, p2));
189
190 }
191
192 Datum
spherepoint_long(PG_FUNCTION_ARGS)193 spherepoint_long(PG_FUNCTION_ARGS)
194 {
195 SPoint *p = (SPoint *) PG_GETARG_POINTER(0);
196
197 PG_RETURN_FLOAT8(p->lng);
198 }
199
200 Datum
spherepoint_lat(PG_FUNCTION_ARGS)201 spherepoint_lat(PG_FUNCTION_ARGS)
202 {
203 SPoint *p = (SPoint *) PG_GETARG_POINTER(0);
204
205 PG_RETURN_FLOAT8(p->lat);
206 }
207
208 Datum
spherepoint_x(PG_FUNCTION_ARGS)209 spherepoint_x(PG_FUNCTION_ARGS)
210 {
211 SPoint *p = (SPoint *) PG_GETARG_POINTER(0);
212 Vector3D v;
213
214 spoint_vector3d(&v, p);
215 PG_RETURN_FLOAT8(v.x);
216 }
217
218 Datum
spherepoint_y(PG_FUNCTION_ARGS)219 spherepoint_y(PG_FUNCTION_ARGS)
220 {
221 SPoint *p = (SPoint *) PG_GETARG_POINTER(0);
222 Vector3D v;
223
224 spoint_vector3d(&v, p);
225 PG_RETURN_FLOAT8(v.y);
226 }
227
228 Datum
spherepoint_z(PG_FUNCTION_ARGS)229 spherepoint_z(PG_FUNCTION_ARGS)
230 {
231 SPoint *p = (SPoint *) PG_GETARG_POINTER(0);
232 Vector3D v;
233
234 spoint_vector3d(&v, p);
235 PG_RETURN_FLOAT8(v.z);
236 }
237
238 Datum
spherepoint_xyz(PG_FUNCTION_ARGS)239 spherepoint_xyz(PG_FUNCTION_ARGS)
240 {
241 SPoint *p = (SPoint *) PG_GETARG_POINTER(0);
242 Datum dret[3];
243 ArrayType *result;
244 Vector3D v;
245
246 spoint_vector3d(&v, p);
247 dret[0] = Float8GetDatumFast(v.x);
248 dret[1] = Float8GetDatumFast(v.y);
249 dret[2] = Float8GetDatumFast(v.z);
250 #ifdef USE_FLOAT8_BYVAL
251 result = construct_array(dret, 3, FLOAT8OID, sizeof(float8), true, 'd');
252 #else
253 result = construct_array(dret, 3, FLOAT8OID, sizeof(float8), false, 'd');
254 #endif
255 PG_RETURN_ARRAYTYPE_P(result);
256 }
257
258 Datum
spherepoint_equal(PG_FUNCTION_ARGS)259 spherepoint_equal(PG_FUNCTION_ARGS)
260 {
261 SPoint *p1 = (SPoint *) PG_GETARG_POINTER(0);
262 SPoint *p2 = (SPoint *) PG_GETARG_POINTER(1);
263
264 PG_RETURN_BOOL(spoint_eq(p1, p2));
265 }
266