pgLatLon
changeset 24:f765aed73421
Bugfix: consider meridian compression when projecting points onto edges of polygons/paths/outlines
author | jbe |
---|---|
date | Thu Sep 22 21:22:08 2016 +0200 (2016-09-22) |
parents | 6b4bdf99eb64 |
children | d8ac13006162 |
files | latlon-v0006.c |
line diff
1.1 --- a/latlon-v0006.c Wed Sep 21 19:32:39 2016 +0200 1.2 +++ b/latlon-v0006.c Thu Sep 22 21:22:08 2016 +0200 1.3 @@ -802,6 +802,7 @@ 1.4 1.5 /* calculate (approximate) distance between point and cluster */ 1.6 static double pgl_point_cluster_distance(pgl_point *point, pgl_cluster *cluster) { 1.7 + double comp; /* square of compression of meridians */ 1.8 int i, j, k; /* i: entry, j: point in entry, k: next point in entry */ 1.9 int entrytype; /* type of entry */ 1.10 int npoints; /* number of points in entry */ 1.11 @@ -819,6 +820,10 @@ 1.12 double min_dist = INFINITY; /* minimum distance */ 1.13 /* distance is zero if point is contained in cluster */ 1.14 if (pgl_point_in_cluster(point, cluster, false)) return 0; 1.15 + /* calculate (approximate) square compression of meridians */ 1.16 + /* TODO: use more exact formula based on WGS-84 */ 1.17 + comp = cos((lat0 / 180.0) * M_PI); 1.18 + comp *= comp; 1.19 /* iterate over all entries */ 1.20 for (i=0; i<cluster->nentries; i++) { 1.21 /* get properties of entry */ 1.22 @@ -887,8 +892,8 @@ 1.23 if (lat1 == lat2 && lon1 == lon2) continue; 1.24 /* otherwise test if point can be projected onto edge of polygon */ 1.25 s = ( 1.26 - ((lat0-lat1) * (lat2-lat1) + (lon0-lon1) * (lon2-lon1)) / 1.27 - ((lat2-lat1) * (lat2-lat1) + (lon2-lon1) * (lon2-lon1)) 1.28 + ((lat0-lat1) * (lat2-lat1) + comp * (lon0-lon1) * (lon2-lon1)) / 1.29 + ((lat2-lat1) * (lat2-lat1) + comp * (lon2-lon1) * (lon2-lon1)) 1.30 ); 1.31 /* go to next vertex and edge if point cannot be projected */ 1.32 if (!(s > 0 && s < 1)) continue;