Skip to content

Commit b02ac6e

Browse files
derickrderickr
committed
Added interpolation for linestrings and polygon.
1 parent 31fe8f4 commit b02ac6e

File tree

4 files changed

+207
-34
lines changed

4 files changed

+207
-34
lines changed

‎geo_array.c

Lines changed: 16 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -29,12 +29,28 @@ geo_array *geo_array_ctor(int element_count)
2929

3030
tmp = malloc(sizeof(geo_array));
3131
tmp->count = element_count;
32+
tmp->allocated = element_count;
3233
tmp->status = calloc(1, element_count);
3334
tmp->x = (double*) calloc(1, element_count * sizeof(double));
3435
tmp->y = (double*) calloc(1, element_count * sizeof(double));
3536

3637
return tmp;
3738
}
39+
40+
void geo_array_add(geo_array *points, double lat, double lon)
41+
{
42+
if (points->count >= points->allocated) {
43+
points->allocated = 1 + (points->allocated * 2);
44+
points->status = realloc(points->status, points->allocated);
45+
points->x = (double*) realloc(points->x, points->allocated * sizeof(double));
46+
points->y = (double*) realloc(points->y, points->allocated * sizeof(double));
47+
}
48+
points->x[points->count] = lat;
49+
points->y[points->count] = lon;
50+
points->status[points->count] = 1;
51+
52+
points->count++;
53+
}
3854

3955
void geo_array_dtor(geo_array *points)
4056
{

‎geo_array.h

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -24,7 +24,9 @@ typedef struct geo_array {
2424
double *y;
2525
char *status;
2626
int count;
27+
int allocated;
2728
} geo_array;
2829

2930
geo_array *geo_array_ctor(int element_count);
31+
void geo_array_add(geo_array *points, double lat, double lon);
3032
void geo_array_dtor(geo_array *points);

‎geospatial.c

Lines changed: 187 additions & 34 deletions
Original file line numberDiff line numberDiff line change
@@ -97,6 +97,16 @@ ZEND_BEGIN_ARG_INFO_EX(rdp_simplify_args, 0, 0, 2)
9797
ZEND_ARG_INFO(0, epsilon)
9898
ZEND_END_ARG_INFO()
9999

100+
ZEND_BEGIN_ARG_INFO_EX(interpolate_linestring_args, 0, 0, 2)
101+
ZEND_ARG_INFO(0, GeoJSONLineString)
102+
ZEND_ARG_INFO(0, epsilon)
103+
ZEND_END_ARG_INFO()
104+
105+
ZEND_BEGIN_ARG_INFO_EX(interpolate_polygon_args, 0, 0, 2)
106+
ZEND_ARG_INFO(0, GeoJSONPolygon)
107+
ZEND_ARG_INFO(0, epsilon)
108+
ZEND_END_ARG_INFO()
109+
100110
/* {{{ geospatial_functions[]
101111
*
102112
* Every user visible function must have an entry in geospatial_functions[].
@@ -113,6 +123,8 @@ const zend_function_entry geospatial_functions[] = {
113123
PHP_FE(decimal_to_dms, decimal_to_dms_args)
114124
PHP_FE(vincenty, vincenty_args)
115125
PHP_FE(rdp_simplify, rdp_simplify_args)
126+
PHP_FE(interpolate_linestring, interpolate_linestring_args)
127+
PHP_FE(interpolate_polygon, interpolate_polygon_args)
116128
/* End of functions */
117129
{ NULL, NULL, NULL }
118130
};
@@ -186,24 +198,12 @@ void retval_point_from_coordinates(zval *return_value, double lon, double lat)
186198
add_assoc_zval_ex(return_value, "coordinates", sizeof("coordinates"), coordinates);
187199
}
188200

189-
int geojson_point_to_lon_lat(zval *point, double *lon, double *lat)
201+
static int parse_point_pair(zval *coordinates, double *lon, double *lat)
190202
{
191-
zval **type, **coordinates, **z_lon, **z_lat;
192203
HashTable *coords;
204+
zval **z_lon, **z_lat;
193205

194-
if (zend_hash_find(HASH_OF(point), "type", sizeof("type"), (void**) &type) != SUCCESS) {
195-
return 0;
196-
}
197-
if (Z_TYPE_PP(type) != IS_STRING || strcmp(Z_STRVAL_PP(type), "Point") != 0) {
198-
return 0;
199-
}
200-
if (zend_hash_find(HASH_OF(point), "coordinates", sizeof("coordinates"), (void**) &coordinates) != SUCCESS) {
201-
return 0;
202-
}
203-
if (Z_TYPE_PP(coordinates) != IS_ARRAY) {
204-
return 0;
205-
}
206-
coords = HASH_OF(*coordinates);
206+
coords = HASH_OF(coordinates);
207207
if (coords->nNumOfElements != 2) {
208208
return 0;
209209
}
@@ -220,6 +220,26 @@ int geojson_point_to_lon_lat(zval *point, double *lon, double *lat)
220220
return 1;
221221
}
222222

223+
int geojson_point_to_lon_lat(zval *point, double *lon, double *lat)
224+
{
225+
zval **type, **coordinates;
226+
227+
if (zend_hash_find(HASH_OF(point), "type", sizeof("type"), (void**) &type) != SUCCESS) {
228+
return 0;
229+
}
230+
if (Z_TYPE_PP(type) != IS_STRING || strcmp(Z_STRVAL_PP(type), "Point") != 0) {
231+
return 0;
232+
}
233+
if (zend_hash_find(HASH_OF(point), "coordinates", sizeof("coordinates"), (void**) &coordinates) != SUCCESS) {
234+
return 0;
235+
}
236+
if (Z_TYPE_PP(coordinates) != IS_ARRAY) {
237+
return 0;
238+
}
239+
240+
return parse_point_pair(*coordinates, lon, lat);
241+
}
242+
223243
/* }}} */
224244

225245
double php_geo_haversine(double from_lat, double from_long, double to_lat, double to_long)
@@ -516,7 +536,9 @@ PHP_FUNCTION(transform_datum)
516536
return;
517537
}
518538

519-
geojson_point_to_lon_lat(geojson, &longitude, &latitude);
539+
if (!geojson_point_to_lon_lat(geojson, &longitude, &latitude)) {
540+
RETURN_FALSE;
541+
}
520542

521543
geo_ellipsoid eli_from = get_ellipsoid(from_reference_ellipsoid);
522544
geo_ellipsoid eli_to = get_ellipsoid(to_reference_ellipsoid);
@@ -630,8 +652,8 @@ var brng = Math.atan2(y, x).toDeg();
630652
*/
631653
double x, y;
632654

633-
y = sin(to_long - from_long) * cos(to_lat);
634-
x = cos(from_lat) * sin(to_lat) - sin(from_lat) * cos(to_lat) * cos(to_long - to_long);
655+
y = sin(fabs(to_long - from_long)) * cos(to_lat);
656+
x = (cos(from_lat) * sin(to_lat)) - (sin(from_lat) * cos(to_lat) * cos(fabs(to_long - from_long)));
635657

636658
return atan2(y, x);
637659
}
@@ -667,7 +689,7 @@ geo_array *geo_hashtable_to_array(zval *array)
667689
int element_count;
668690
HashPosition pos;
669691
zval **entry;
670-
zval **z_lon, **z_lat;
692+
double lon, lat;
671693
int i = 0;
672694

673695
element_count = zend_hash_num_elements(Z_ARRVAL_P(array));
@@ -676,24 +698,12 @@ geo_array *geo_hashtable_to_array(zval *array)
676698
zend_hash_internal_pointer_reset_ex(Z_ARRVAL_P(array), &pos);
677699
while (zend_hash_get_current_data_ex(Z_ARRVAL_P(array), (void **)&entry, &pos) == SUCCESS) {
678700

679-
if (Z_TYPE_PP(entry) != IS_ARRAY) {
680-
goto failure;
681-
}
682-
if (zend_hash_num_elements(Z_ARRVAL_PP(entry)) != 2)
683-
{
701+
if (!parse_point_pair(*entry, &lon, &lat)) {
684702
goto failure;
685703
}
686-
if (zend_hash_index_find(HASH_OF(*entry), 0, (void**) &z_lon) != SUCCESS) {
687-
return 0;
688-
}
689-
if (zend_hash_index_find(HASH_OF(*entry), 1, (void**) &z_lat) != SUCCESS) {
690-
return 0;
691-
}
692-
convert_to_double_ex(z_lon);
693-
convert_to_double_ex(z_lat);
694704

695-
tmp->x[i] = Z_DVAL_PP(z_lon);
696-
tmp->y[i] = Z_DVAL_PP(z_lat);
705+
tmp->x[i] = lon;
706+
tmp->y[i] = lat;
697707
tmp->status[i] = 1;
698708

699709
zend_hash_move_forward_ex(Z_ARRVAL_P(array), &pos);
@@ -707,6 +717,37 @@ geo_array *geo_hashtable_to_array(zval *array)
707717
return NULL;
708718
}
709719

720+
int geojson_linestring_to_array(zval *line, geo_array **array)
721+
{
722+
geo_array *tmp;
723+
zval **type, **coordinates;
724+
725+
if (Z_TYPE_P(line) != IS_ARRAY) {
726+
return 0;
727+
}
728+
729+
if (zend_hash_find(HASH_OF(line), "type", sizeof("type"), (void**) &type) != SUCCESS) {
730+
return 0;
731+
}
732+
if (Z_TYPE_PP(type) != IS_STRING || strcmp(Z_STRVAL_PP(type), "Linestring") != 0) {
733+
return 0;
734+
}
735+
if (zend_hash_find(HASH_OF(line), "coordinates", sizeof("coordinates"), (void**) &coordinates) != SUCCESS) {
736+
return 0;
737+
}
738+
if (Z_TYPE_PP(coordinates) != IS_ARRAY) {
739+
return 0;
740+
}
741+
742+
tmp = geo_hashtable_to_array(*coordinates);
743+
if (tmp && array) {
744+
*array = tmp;
745+
return 1;
746+
}
747+
748+
return 0;
749+
}
750+
710751
double rdp_find_perpendicular_distable(double pX, double pY, double p1X, double p1Y, double p2X, double p2Y)
711752
{
712753
double slope, intercept, result;
@@ -797,6 +838,118 @@ PHP_FUNCTION(rdp_simplify)
797838
}
798839
/* }}} */
799840

841+
static geo_array *interpolate_line(geo_array *points, double epsilon)
842+
{
843+
int i;
844+
geo_array *new_array;
845+
double dx, dy, distance, step_size, res_lat, res_long, fraction;
846+
847+
new_array = geo_array_ctor(0);
848+
849+
for (i = 0; i < points->count - 1; i++) {
850+
dx = fabs(points->x[i] - points->x[i + 1]);
851+
dy = fabs(points->y[i] - points->y[i + 1]);
852+
distance = sqrt((dx * dx) + (dy * dy));
853+
if (distance > epsilon) {
854+
step_size = epsilon/distance;
855+
for (fraction = 0; fraction < 1; fraction += step_size) {
856+
php_geo_fraction_along_gc_line(
857+
points->y[i] * GEO_DEG_TO_RAD,
858+
points->x[i] * GEO_DEG_TO_RAD,
859+
points->y[i + 1] * GEO_DEG_TO_RAD,
860+
points->x[i + 1] * GEO_DEG_TO_RAD,
861+
fraction, GEO_EARTH_RADIUS,
862+
&res_lat, &res_long
863+
);
864+
geo_array_add(new_array, res_long / GEO_DEG_TO_RAD, res_lat / GEO_DEG_TO_RAD);
865+
}
866+
} else {
867+
geo_array_add(new_array, points->x[i], points->y[i]);
868+
}
869+
}
870+
geo_array_add(new_array, points->x[points->count - 1], points->y[points->count - 1]);
871+
872+
return new_array;
873+
}
874+
875+
/* {{{ proto array interpolate_linestring(GeoJSONLineString line, float epsilon)
876+
Interpolates lines with intermediate points to show line segments as GC lines */
877+
PHP_FUNCTION(interpolate_linestring)
878+
{
879+
zval *line;
880+
double epsilon;
881+
geo_array *points;
882+
int i;
883+
zval *pair;
884+
geo_array *new_array;
885+
886+
if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, "zd", &line, &epsilon) == FAILURE) {
887+
return;
888+
}
889+
890+
if (!geojson_linestring_to_array(line, &points)) {
891+
RETURN_FALSE;
892+
}
893+
894+
array_init(return_value);
895+
896+
new_array = interpolate_line(points, epsilon);
897+
898+
for (i = 0; i < new_array->count; i++) {
899+
if (new_array->status[i]) {
900+
MAKE_STD_ZVAL(pair);
901+
array_init(pair);
902+
add_next_index_double(pair, new_array->x[i]);
903+
add_next_index_double(pair, new_array->y[i]);
904+
add_next_index_zval(return_value, pair);
905+
}
906+
}
907+
908+
geo_array_dtor(points);
909+
geo_array_dtor(new_array);
910+
}
911+
/* }}} */
912+
913+
/* {{{ proto array interpolate_polygon(GeoJSONPolygon polygon, float epsilon)
914+
Interpolates polygons with intermediate points to show line segments as GC lines */
915+
PHP_FUNCTION(interpolate_polygon)
916+
{
917+
zval *polygon;
918+
double epsilon;
919+
geo_array *points;
920+
int i;
921+
zval *pair;
922+
923+
if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, "zd", &polygon, &epsilon) == FAILURE) {
924+
return;
925+
}
926+
927+
if (!Z_TYPE_P(polygon) == IS_ARRAY) {
928+
return;
929+
}
930+
931+
if (!geojson_linestring_to_array(polygon, &points)) {
932+
RETURN_FALSE;
933+
}
934+
935+
array_init(return_value);
936+
937+
rdp_simplify(points, epsilon, 0, points->count - 1);
938+
939+
for (i = 0; i < points->count; i++) {
940+
if (points->status[i]) {
941+
MAKE_STD_ZVAL(pair);
942+
array_init(pair);
943+
add_next_index_double(pair, points->x[i]);
944+
add_next_index_double(pair, points->y[i]);
945+
add_next_index_zval(return_value, pair);
946+
}
947+
}
948+
949+
geo_array_dtor(points);
950+
}
951+
/* }}} */
952+
800953
/*
801954
* Local variables:
802955
* tab-width: 4

‎php_geospatial.h

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -139,6 +139,8 @@ PHP_FUNCTION(dms_to_decimal);
139139
PHP_FUNCTION(decimal_to_dms);
140140
PHP_FUNCTION(vincenty);
141141
PHP_FUNCTION(rdp_simplify);
142+
PHP_FUNCTION(interpolate_linestring);
143+
PHP_FUNCTION(interpolate_polygon);
142144

143145
#endif /* PHP_GEOSPATIAL_H */
144146

0 commit comments

Comments
 (0)