geoencode.cc

Go to the documentation of this file.

/* Copyright (C) 2011 Richard Boulton
 * Based closely on a python version, copyright (C) 2010 Olly Betts
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to
 * deal in the Software without restriction, including without limitation the
 * rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
 * sell copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in
 * all copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
 * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS
 * IN THE SOFTWARE.
 */

#include <config.h>
#include "geoencode.h"

#include <cmath>

using namespace std;

struct DegreesMinutesSeconds {
    int degrees;

    int minutes;

    int seconds;

    int sec16ths;

    explicit DegreesMinutesSeconds(int angle_16th_secs) {
        degrees = angle_16th_secs / (3600 * 16);
        angle_16th_secs = angle_16th_secs % (3600 * 16);
        minutes = angle_16th_secs / (60 * 16);
        angle_16th_secs = angle_16th_secs % (60 * 16);
        seconds = angle_16th_secs / 16;
        sec16ths = angle_16th_secs % 16;
    }
};

bool
GeoEncode::encode(double lat, double lon, string & result)
{
    // Check range of latitude.
    if (rare(lat < -90.0 || lat > 90.0)) {
        return false;
    }

    // Wrap longitude to range [0,360).
    lon = fmod(lon, 360.0);
    if (lon < 0) {
        lon += 360;
    }

    int lat_16ths, lon_16ths;
    lat_16ths = lround((lat + 90.0) * 57600.0);
    if (lat_16ths == 0 || lat_16ths == 57600 * 180) {
        lon_16ths = 0;
    } else {
        lon_16ths = lround(lon * 57600.0);
        if (lon_16ths == 57600 * 360) {
            lon_16ths = 0;
        }
    }

    DegreesMinutesSeconds lat_dms(lat_16ths);
    DegreesMinutesSeconds lon_dms(lon_16ths);

    size_t old_len = result.size();
    result.resize(old_len + 6);

    // Add degrees parts as first two bytes.
    unsigned dd = lat_dms.degrees + lon_dms.degrees * 181;
    // dd is in range 0..180*360+359 = 0..65159
    result[old_len] = char(dd >> 8);
    result[old_len + 1] = char(dd & 0xff);

    // Add minutes next; 4 bits from each in the first byte.
    result[old_len + 2] = char(((lat_dms.minutes / 4) << 4) |
                               (lon_dms.minutes / 4)
                              );

    result[old_len + 3] = char(
                               ((lat_dms.minutes % 4) << 6) |
                               ((lon_dms.minutes % 4) << 4) |
                               ((lat_dms.seconds / 15) << 2) |
                               (lon_dms.seconds / 15)
                              );

    result[old_len + 4] = char(
                               ((lat_dms.seconds % 15) << 4) |
                               (lon_dms.seconds % 15)
                              );

    result[old_len + 5] = char(
                               (lat_dms.sec16ths << 4) |
                               lon_dms.sec16ths
                              );

    return true;
}

void
GeoEncode::decode(const char * value, size_t len,
                  double & lat_ref, double & lon_ref)
{
    const unsigned char * ptr
            = reinterpret_cast<const unsigned char *>(value);
    unsigned tmp = (ptr[0] & 0xff) << 8 | (ptr[1] & 0xff);
    lat_ref = tmp % 181;
    lon_ref = tmp / 181;
    if (len > 2) {
        tmp = ptr[2];
        double lat_m = (tmp >> 4) * 4;
        double lon_m = (tmp & 0xf) * 4;

        if (len > 3) {
            tmp = ptr[3];
            lat_m += (tmp >> 6) & 3;
            lon_m += (tmp >> 4) & 3;
            double lat_s = ((tmp >> 2) & 3) * 15;
            double lon_s = (tmp & 3) * 15;

            if (len > 4) {
                tmp = ptr[4];
                lat_s += (tmp >> 4) & 0xf;
                lon_s += tmp & 0xf;

                if (len > 5) {
                    tmp = ptr[5];
                    lat_s += ((tmp >> 4) / 16.0);
                    lon_s += ((tmp & 0xf) / 16.0);
                }
            }

            lat_m += lat_s / 60.0;
            lon_m += lon_s / 60.0;
        }

        lat_ref += lat_m / 60.0;
        lon_ref += lon_m / 60.0;
    }

    lat_ref -= 90.0;
}

static void
calc_latlon_16ths(double lat, double lon, int & lat_16ths, int & lon_16ths)
{
    lat_16ths = lround((lat + 90.0) * 57600.0);
    lon_16ths = lround(lon * 57600.0);
    if (lon_16ths == 57600 * 360) {
        lon_16ths = 0;
    }
}

GeoEncode::DecoderWithBoundingBox::DecoderWithBoundingBox
(double lat1, double lon1_, double lat2, double lon2_)
        : lon1(lon1_), lon2(lon2_),
          min_lat(lat1), max_lat(lat2),
          include_poles(false)
{
    // Wrap longitudes to range [0,360).
    lon1 = fmod(lon1, 360.0);
    if (lon1 < 0) {
        lon1 += 360;
    }

    lon2 = fmod(lon2, 360.0);
    if (lon2 < 0) {
        lon2 += 360;
    }

    // Calculate start1
    int lat_16ths, lon_16ths;
    calc_latlon_16ths(lat1, lon1, lat_16ths, lon_16ths);
    if (lat_16ths == 0 || lat_16ths == 57600 * 180) {
        include_poles = true;
    }
    unsigned dd = lat_16ths / (3600 * 16) + (lon_16ths / (3600 * 16)) * 181;
    start1 = char(dd >> 8);

    calc_latlon_16ths(lat2, lon2, lat_16ths, lon_16ths);
    if (lat_16ths == 0 || lat_16ths == 57600 * 180) {
        include_poles = true;
    }
    dd = lat_16ths / (3600 * 16) + (lon_16ths / (3600 * 16)) * 181;
    start2 = char(dd >> 8);

    discontinuous_longitude_range = (lon1 > lon2);
}

bool
GeoEncode::DecoderWithBoundingBox::decode(const std::string & value,
                                          double & lat_ref,
                                          double & lon_ref) const
{
    unsigned char start = value[0];
    if (discontinuous_longitude_range) {
        // start must be outside range of (start2..start1)
        // (start2 will be > start1)
        if (start2 < start && start < start1) {
            if (!(include_poles && start == 0))
                return false;
        }
    } else {
        // start must be inside range of [start1..start2] (inclusive of ends).
        if (start < start1 || start2 < start) {
            if (!(include_poles && start == 0))
                return false;
        }
    }
    double lat, lon;
    GeoEncode::decode(value, lat, lon);
    if (lat < min_lat || lat > max_lat) {
        return false;
    }
    if (lat == -90 || lat == 90) {
        // It's a pole, so the longitude isn't meaningful (will be zero)
        // and we've already checked that the latitude is in range.
        lat_ref = lat;
        lon_ref = 0;
        return true;
    }
    if (discontinuous_longitude_range) {
        if (lon2 < lon && lon < lon1)
            return false;
    } else {
        if (lon < lon1 || lon2 < lon)
            return false;
    }

    lat_ref = lat;
    lon_ref = lon;
    return true;
}