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;