Skip to content

Commit 14d61d5

Browse files
author
Nathaniel McHugh
committed
* better set up of constants
* more tests
1 parent 51c5910 commit 14d61d5

File tree

5 files changed

+172
-50
lines changed

5 files changed

+172
-50
lines changed

‎geospatial.c

Lines changed: 114 additions & 37 deletions
Original file line numberDiff line numberDiff line change
@@ -54,6 +54,13 @@ ZEND_BEGIN_ARG_INFO_EX(cartesian_to_polar_args, 0, 0, 4)
5454
ZEND_ARG_INFO(0, reference_ellipsoid)
5555
ZEND_END_ARG_INFO()
5656

57+
ZEND_BEGIN_ARG_INFO_EX(change_datum_args, 0, 0, 4)
58+
ZEND_ARG_INFO(0, latitude)
59+
ZEND_ARG_INFO(0, longitude)
60+
ZEND_ARG_INFO(0, from_reference_ellipsoid)
61+
ZEND_ARG_INFO(0, to_reference_ellipsoid)
62+
ZEND_END_ARG_INFO()
63+
5764
/* {{{ geospatial_functions[]
5865
*
5966
* Every user visible function must have an entry in geospatial_functions[].
@@ -63,6 +70,7 @@ const zend_function_entry geospatial_functions[] = {
6370
PHP_FE(helmert, helmert_args)
6471
PHP_FE(polar_to_cartesian, polar_to_cartesian_args)
6572
PHP_FE(cartesian_to_polar, cartesian_to_polar_args)
73+
PHP_FE(change_datum, change_datum_args)
6674
/* End of functions */
6775
{ NULL, NULL, NULL }
6876
};
@@ -122,17 +130,25 @@ PHP_MINFO_FUNCTION(geospatial)
122130
}
123131
/* }}} */
124132

125-
PHP_FUNCTION(helmert)
133+
geo_ellipsoid get_ellipsoid(long ellipsoid_const)
126134
{
127-
double x, y, z;
128-
double xOut, yOut, zOut;
129-
double rX, rY, rZ;
130-
131-
if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, "ddd", &x, &y, &z) == FAILURE) {
132-
return;
135+
switch (ellipsoid_const) {
136+
case GEO_AIRY_1830:
137+
return airy_1830;
138+
break;
139+
case GEO_WGS84:
140+
default:
141+
return wgs84;
142+
break;
133143
}
144+
}
145+
146+
geo_cartesian php_helmert(double x, double y, double z)
147+
{
148+
double rX, rY, rZ;
149+
double xOut, yOut, zOut;
150+
geo_cartesian point;
134151

135-
array_init(return_value);
136152
rX = ROTATION_X / GEO_SEC_IN_DEG * GEO_DEG_TO_RAD;
137153
rY = ROTATION_Y / GEO_SEC_IN_DEG * GEO_DEG_TO_RAD;
138154
rZ = ROTATION_Z / GEO_SEC_IN_DEG * GEO_DEG_TO_RAD;
@@ -146,25 +162,33 @@ PHP_FUNCTION(helmert)
146162
zOut = WGS84_OSGB36_Z;
147163
zOut += ((-1 * rY * x) + (rX * y) + z) * SCALE_CHANGE;
148164

149-
add_assoc_double(return_value, "x", xOut);
150-
add_assoc_double(return_value, "y", yOut);
151-
add_assoc_double(return_value, "z", zOut);
165+
point.x = xOut;
166+
point.y = yOut;
167+
point.z = zOut;
168+
return point;
152169
}
153170

154-
PHP_FUNCTION(polar_to_cartesian)
171+
172+
PHP_FUNCTION(helmert)
155173
{
156-
double latitude, longitude;
157174
double x, y, z;
158-
long reference_ellipsoid;
159-
if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, "dd|l", &latitude, &longitude, &reference_ellipsoid) == FAILURE) {
175+
geo_cartesian point;
176+
if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, "ddd", &x, &y, &z) == FAILURE) {
160177
return;
161178
}
162-
geo_ellipsoid eli = {WGS84_A, WGS84_B};
163-
if (reference_ellipsoid == GEO_AIRY_1830) {
164-
eli.a = AIRY_1830_A;
165-
eli.b = AIRY_1830_B;
166-
}
179+
167180
array_init(return_value);
181+
point = php_helmert(x, y, z);
182+
add_assoc_double(return_value, "x", point.x);
183+
add_assoc_double(return_value, "y", point.y);
184+
add_assoc_double(return_value, "z", point.z);
185+
}
186+
187+
geo_cartesian php_polar_to_cartesian(double latitude, double longitude, geo_ellipsoid eli)
188+
{
189+
double x, y, z;
190+
191+
geo_cartesian point;
168192
double phi = latitude * GEO_DEG_TO_RAD;
169193
double lambda = longitude * GEO_DEG_TO_RAD;
170194
double eSq = ((eli.a * eli.a) - (eli.b * eli.b)) / (eli.a * eli.a);
@@ -175,43 +199,96 @@ PHP_FUNCTION(polar_to_cartesian)
175199
y *= cos(phi) * sin(lambda);
176200
z = ((1 - eSq) * nu) + HEIGHT;
177201
z*= sin(phi);
178-
add_assoc_double(return_value, "x", x);
179-
add_assoc_double(return_value, "y", y);
180-
add_assoc_double(return_value, "z", z);
202+
point.x = x;
203+
point.y = y;
204+
point.z = z;
205+
return point;
181206
}
182207

183-
PHP_FUNCTION(cartesian_to_polar)
208+
PHP_FUNCTION(polar_to_cartesian)
184209
{
185210
double latitude, longitude;
186-
double x, y, z;
187-
double nu, lambda, h;
188211
long reference_ellipsoid;
189-
if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, "ddd|l", &x, &y, &z, &reference_ellipsoid) == FAILURE) {
212+
geo_cartesian point;
213+
214+
if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, "dd|l", &latitude, &longitude, &reference_ellipsoid) == FAILURE) {
190215
return;
191216
}
192-
geo_ellipsoid eli = {WGS84_A, WGS84_B};
193-
if (reference_ellipsoid == GEO_AIRY_1830) {
194-
eli.a = AIRY_1830_A;
195-
eli.b = AIRY_1830_B;
196-
}
217+
geo_ellipsoid eli = get_ellipsoid(reference_ellipsoid);
218+
array_init(return_value);
219+
point = php_polar_to_cartesian(latitude, longitude, eli);
220+
add_assoc_double(return_value, "x", point.x);
221+
add_assoc_double(return_value, "y", point.y);
222+
add_assoc_double(return_value, "z", point.z);
223+
}
224+
225+
geo_lat_long php_cartesian_to_polar(double x, double y, double z, geo_ellipsoid eli)
226+
{
227+
228+
double latitude, longitude;
229+
double nu, lambda, h;
230+
geo_lat_long polar;
231+
197232
//aiming for 1m accuracy
198233
double precision = 0.1 / eli.a;
199-
array_init(return_value);
200234
double eSq = ((eli.a * eli.a) - (eli.b * eli.b)) / (eli.a * eli.a);
201235
double p = sqrt(x * x + y * y);
202236
double phi = atan2(z, p * (1 - eSq));
203237
double phiP = 2 * M_PI;
204238
while (abs(phi - phiP) > precision) {
205-
nu = AIRY_1830_A / sqrt(1 - eSq * sin(phi) * sin(phi));
239+
nu = eli.a / sqrt(1 - eSq * sin(phi) * sin(phi));
206240
phiP = phi;
207241
phi = atan2(z + eSq * nu * sin(phi), p);
208242
}
209243
lambda = atan2(y ,x);
210244
h = p / cos(phi) - nu;
245+
polar.latitude = phi / GEO_DEG_TO_RAD;
246+
polar.longitude = lambda / GEO_DEG_TO_RAD;
247+
polar.height = h;
248+
249+
return polar;
250+
}
211251

212-
add_assoc_double(return_value, "lat", phi / GEO_DEG_TO_RAD);
213-
add_assoc_double(return_value, "long", lambda / GEO_DEG_TO_RAD);
214-
add_assoc_double(return_value, "height", h);
252+
PHP_FUNCTION(cartesian_to_polar)
253+
{
254+
double x, y, z;
255+
long reference_ellipsoid;
256+
geo_lat_long polar;
257+
258+
if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, "ddd|l", &x, &y, &z, &reference_ellipsoid) == FAILURE) {
259+
return;
260+
}
261+
geo_ellipsoid eli = get_ellipsoid(reference_ellipsoid);
262+
array_init(return_value);
263+
polar = php_cartesian_to_polar(x, y, z, eli);
264+
add_assoc_double(return_value, "lat", polar.latitude);
265+
add_assoc_double(return_value, "long", polar.longitude);
266+
add_assoc_double(return_value, "height", polar.height);
267+
}
268+
269+
270+
271+
272+
273+
PHP_FUNCTION(change_datum)
274+
{
275+
double latitude, longitude;
276+
long from_reference_ellipsoid, to_reference_ellipsoid;
277+
geo_cartesian point, converted_point;
278+
geo_lat_long polar;
279+
if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, "ddll", &latitude, &longitude, &from_reference_ellipsoid, &to_reference_ellipsoid) == FAILURE) {
280+
return;
281+
}
282+
geo_ellipsoid eli_from = get_ellipsoid(from_reference_ellipsoid);
283+
geo_ellipsoid eli_to = get_ellipsoid(to_reference_ellipsoid);
284+
point = php_polar_to_cartesian(latitude, longitude, eli_from);
285+
converted_point = php_helmert(point.x, point.y, point.z);
286+
polar = php_cartesian_to_polar(converted_point.x, converted_point.y, converted_point.z, eli_to);
287+
288+
array_init(return_value);
289+
add_assoc_double(return_value, "lat", polar.latitude);
290+
add_assoc_double(return_value, "long", polar.longitude);
291+
add_assoc_double(return_value, "height", polar.height);
215292
}
216293

217294
/* {{{ proto haversine(double fromLat, double fromLong, double toLat, double toLong [, double radius ])

‎php_geospatial.h

Lines changed: 13 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -37,15 +37,24 @@ extern zend_module_entry geospatial_module_entry;
3737
#endif
3838

3939
typedef struct {
40-
double latitude;
41-
double longitude;
42-
} latLong;
40+
double latitude;
41+
double longitude;
42+
double height;
43+
} geo_lat_long;
4344

4445
typedef struct {
4546
double a;
4647
double b;
4748
} geo_ellipsoid;
4849

50+
typedef struct {
51+
double x;
52+
double y;
53+
double z;
54+
} geo_cartesian;
55+
56+
const geo_ellipsoid wgs84 = {6378137.000, 6356752.3142};
57+
const geo_ellipsoid airy_1830 = {6377563.396, 6356256.910};
4958

5059
#define GEO_DEG_TO_RAD 0.017453292519943295769236907684886
5160
/**
@@ -58,13 +67,6 @@ typedef struct {
5867
#define GEO_WGS84 0x0001
5968
#define GEO_AIRY_1830 0x0002
6069

61-
62-
#define AIRY_1830_A 6377563.396
63-
#define AIRY_1830_B 6356256.910
64-
65-
#define WGS84_A 6378137.000
66-
#define WGS84_B 6356752.3142
67-
6870
#define WGS84_OSGB36_X -446.448
6971
#define WGS84_OSGB36_Y 125.157
7072
#define WGS84_OSGB36_Z -542.060
@@ -85,6 +87,7 @@ PHP_FUNCTION(haversine);
8587
PHP_FUNCTION(helmert);
8688
PHP_FUNCTION(polar_to_cartesian);
8789
PHP_FUNCTION(cartesian_to_polar);
90+
PHP_FUNCTION(change_datum);
8891

8992
#endif /* PHP_GEOSPATIAL_H */
9093

‎tests/JodrellBank.phpt

Lines changed: 22 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,22 @@
1+
--TEST--
2+
WGS84 to OSGB36
3+
--FILE--
4+
<?php
5+
6+
function getDecimal($d, $m, $s) {
7+
$mul = $d < 0 ? -1 : 1;
8+
return $mul * (abs($d)+ $m / 60 + $s /3600);
9+
}
10+
11+
$lat = getDecimal(53, 14, 10.5);
12+
$long = getDecimal(-2, 18, 25.7);
13+
14+
$polar = change_datum($lat, $long, GEO_WGS84, GEO_AIRY_1830);
15+
16+
echo round($polar['lat'] ,6),PHP_EOL;
17+
echo round($polar['long'] ,6),PHP_EOL;
18+
echo round($polar['height'] ,3),PHP_EOL;
19+
--EXPECT--
20+
53.235974
21+
-2.305717
22+
-25.649

‎tests/OSGB36_to_WGS84.phpt

Lines changed: 22 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,22 @@
1+
--TEST--
2+
WGS84 to OSGB36
3+
--FILE--
4+
<?php
5+
6+
function getDecimal($d, $m, $s) {
7+
$mul = $d < 0 ? -1 : 1;
8+
return $mul * (abs($d)+ $m / 60 + $s /3600);
9+
}
10+
11+
$lat = getDecimal(53, 14, 10.5);
12+
$long = getDecimal(-2, 18, 25.7);
13+
14+
$polar = change_datum($lat, $long, GEO_AIRY_1830, GEO_WGS84);
15+
16+
echo round($polar['lat'] ,6),PHP_EOL;
17+
echo round($polar['long'] ,6),PHP_EOL;
18+
echo round($polar['height'] ,3),PHP_EOL;
19+
--EXPECT--
20+
53.236525
21+
-2.308560
22+
75.062

‎tests/WGS84_to_OSGB36.phpt

Lines changed: 1 addition & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -11,10 +11,8 @@ function getDecimal($d, $m, $s) {
1111
$lat = getDecimal(53, 23, 1);
1212
$long = getDecimal(-1, 28, 1);
1313

14-
extract(polar_to_cartesian($lat, $long, GEO_WGS84));
15-
extract(helmert($x, $y, $z));
14+
$polar = change_datum($lat, $long, GEO_WGS84, GEO_AIRY_1830);
1615

17-
$polar = cartesian_to_polar($x, $y, $z, GEO_AIRY_1830);
1816
var_dump($polar);
1917
--EXPECT--
2018
array(3) {

0 commit comments

Comments
 (0)