xapian-core/geospatial/geoencode.cc

   1 /** @file geoencode.cc
   2  * @brief Encodings for geospatial coordinates.
   3  */
   4 /* Copyright (C) 2011 Richard Boulton
   5  * Based closely on a python version, copyright (C) 2010 Olly Betts
   6  *
   7  * Permission is hereby granted, free of charge, to any person obtaining a copy
   8  * of this software and associated documentation files (the "Software"), to
   9  * deal in the Software without restriction, including without limitation the
  10  * rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
  11  * sell copies of the Software, and to permit persons to whom the Software is
  12  * furnished to do so, subject to the following conditions:
  13  *
  14  * The above copyright notice and this permission notice shall be included in
  15  * all copies or substantial portions of the Software.
  16  *
  17  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
  18  * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
  19  * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
  20  * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
  21  * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
  22  * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS
  23  * IN THE SOFTWARE.
  24  */
  25
  26 #include <config.h>
  27 #include "geoencode.h"
  28
  29 #include <cmath>
  30
  31 using namespace std;
  32
  33 /** Angles, split into degrees, minutes and seconds.
  34  *
  35  *  Only designed to work with positive angles.
  36  */
  37 struct DegreesMinutesSeconds {
  38     /** Number of degrees.
  39      *
  40      *  Range 0 <= degrees <= 180 for latitude, 0 <= degrees < 360 for
  41      *  longitude.
  42      */
  43     int degrees;
  44
  45     /** Number of minutes: 0 to 59 */
  46     int minutes;
  47
  48     /** Number of seconds: 0 to 59 */
  49     int seconds;
  50
  51     /** Number of 16ths of a second: 0 to 15 */
  52     int sec16ths;
  53
  54     /** Initialise with a (positive) angle, as an integer representing the
  55      *  number of 16ths of a second, rounding to nearest.
  56      *
  57      *  The range of valid angles is assumed to be 0 <= angle in degrees < 360,
  58      *  so range of angle_16th_secs is 0..20735999, which fits easily into a 32
  59      *  bit int.  (Latitudes are represented in the range 0 <= angle <= 180,
  60      *  where 0 is the south pole.)
  61      */
  62     explicit DegreesMinutesSeconds(int angle_16th_secs) {
  63         degrees = angle_16th_secs / (3600 * 16);
  64         angle_16th_secs = angle_16th_secs % (3600 * 16);
  65         minutes = angle_16th_secs / (60 * 16);
  66         angle_16th_secs = angle_16th_secs % (60 * 16);
  67         seconds = angle_16th_secs / 16;
  68         sec16ths = angle_16th_secs % 16;
  69     }
  70 };
  71
  72 bool
  73 GeoEncode::encode(double lat, double lon, string & result)
  74 {
  75     // Check range of latitude.
  76     if (rare(lat < -90.0 || lat > 90.0)) {
  77         return false;
  78     }
  79
  80     // Wrap longitude to range [0,360).
  81     lon = fmod(lon, 360.0);
  82     if (lon < 0) {
  83         lon += 360;
  84     }
  85
  86     int lat_16ths, lon_16ths;
  87     lat_16ths = lround((lat + 90.0) * 57600.0);
  88     if (lat_16ths == 0 || lat_16ths == 57600 * 180) {
  89         lon_16ths = 0;
  90     } else {
  91         lon_16ths = lround(lon * 57600.0);
  92         if (lon_16ths == 57600 * 360) {
  93             lon_16ths = 0;
  94         }
  95     }
  96
  97     DegreesMinutesSeconds lat_dms(lat_16ths);
  98     DegreesMinutesSeconds lon_dms(lon_16ths);
  99
 100     size_t old_len = result.size();
 101     result.resize(old_len + 6);
 102
 103     // Add degrees parts as first two bytes.
 104     unsigned dd = lat_dms.degrees + lon_dms.degrees * 181;
 105     // dd is in range 0..180*360+359 = 0..65159
 106     result[old_len] = char(dd >> 8);
 107     result[old_len + 1] = char(dd & 0xff);
 108
 109     // Add minutes next; 4 bits from each in the first byte.
 110     result[old_len + 2] = char(((lat_dms.minutes / 4) << 4) |
 111                                (lon_dms.minutes / 4)
 112                               );
 113
 114     result[old_len + 3] = char(
 115                                ((lat_dms.minutes % 4) << 6) |
 116                                ((lon_dms.minutes % 4) << 4) |
 117                                ((lat_dms.seconds / 15) << 2) |
 118                                (lon_dms.seconds / 15)
 119                               );
 120
 121     result[old_len + 4] = char(
 122                                ((lat_dms.seconds % 15) << 4) |
 123                                (lon_dms.seconds % 15)
 124                               );
 125
 126     result[old_len + 5] = char(
 127                                (lat_dms.sec16ths << 4) |
 128                                lon_dms.sec16ths
 129                               );
 130
 131     return true;
 132 }
 133
 134 void
 135 GeoEncode::decode(const char * value, size_t len,
 136                   double & lat_ref, double & lon_ref)
 137 {
 138     const unsigned char * ptr
 139             = reinterpret_cast<const unsigned char *>(value);
 140     unsigned tmp = (ptr[0] & 0xff) << 8 | (ptr[1] & 0xff);
 141     lat_ref = tmp % 181;
 142     lon_ref = tmp / 181;
 143     if (len > 2) {
 144         tmp = ptr[2];
 145         double lat_m = (tmp >> 4) * 4;
 146         double lon_m = (tmp & 0xf) * 4;
 147
 148         if (len > 3) {
 149             tmp = ptr[3];
 150             lat_m += (tmp >> 6) & 3;
 151             lon_m += (tmp >> 4) & 3;
 152             double lat_s = ((tmp >> 2) & 3) * 15;
 153             double lon_s = (tmp & 3) * 15;
 154
 155             if (len > 4) {
 156                 tmp = ptr[4];
 157                 lat_s += (tmp >> 4) & 0xf;
 158                 lon_s += tmp & 0xf;
 159
 160                 if (len > 5) {
 161                     tmp = ptr[5];
 162                     lat_s += ((tmp >> 4) / 16.0);
 163                     lon_s += ((tmp & 0xf) / 16.0);
 164                 }
 165             }
 166
 167             lat_m += lat_s / 60.0;
 168             lon_m += lon_s / 60.0;
 169         }
 170
 171         lat_ref += lat_m / 60.0;
 172         lon_ref += lon_m / 60.0;
 173     }
 174
 175     lat_ref -= 90.0;
 176 }
 177
 178 /// Calc latitude and longitude in integral number of 16ths of a second
 179 static void
 180 calc_latlon_16ths(double lat, double lon, int & lat_16ths, int & lon_16ths)
 181 {
 182     lat_16ths = lround((lat + 90.0) * 57600.0);
 183     lon_16ths = lround(lon * 57600.0);
 184     if (lon_16ths == 57600 * 360) {
 185         lon_16ths = 0;
 186     }
 187 }
 188
 189 GeoEncode::DecoderWithBoundingBox::DecoderWithBoundingBox
 190 (double lat1, double lon1_, double lat2, double lon2_)
 191         : lon1(lon1_), lon2(lon2_),
 192           min_lat(lat1), max_lat(lat2),
 193           include_poles(false)
 194 {
 195     // Wrap longitudes to range [0,360).
 196     lon1 = fmod(lon1, 360.0);
 197     if (lon1 < 0) {
 198         lon1 += 360;
 199     }
 200
 201     lon2 = fmod(lon2, 360.0);
 202     if (lon2 < 0) {
 203         lon2 += 360;
 204     }
 205
 206     // Calculate start1
 207     int lat_16ths, lon_16ths;
 208     calc_latlon_16ths(lat1, lon1, lat_16ths, lon_16ths);
 209     if (lat_16ths == 0 || lat_16ths == 57600 * 180) {
 210         include_poles = true;
 211     }
 212     unsigned dd = lat_16ths / (3600 * 16) + (lon_16ths / (3600 * 16)) * 181;
 213     start1 = char(dd >> 8);
 214
 215     calc_latlon_16ths(lat2, lon2, lat_16ths, lon_16ths);
 216     if (lat_16ths == 0 || lat_16ths == 57600 * 180) {
 217         include_poles = true;
 218     }
 219     dd = lat_16ths / (3600 * 16) + (lon_16ths / (3600 * 16)) * 181;
 220     start2 = char(dd >> 8);
 221
 222     discontinuous_longitude_range = (lon1 > lon2);
 223 }
 224
 225 bool
 226 GeoEncode::DecoderWithBoundingBox::decode(const std::string & value,
 227                                           double & lat_ref,
 228                                           double & lon_ref) const
 229 {
 230     unsigned char start = value[0];
 231     if (discontinuous_longitude_range) {
 232         // start must be outside range of (start2..start1)
 233         // (start2 will be > start1)
 234         if (start2 < start && start < start1) {
 235             if (!(include_poles && start == 0))
 236                 return false;
 237         }
 238     } else {
 239         // start must be inside range of [start1..start2] (inclusive of ends).
 240         if (start < start1 || start2 < start) {
 241             if (!(include_poles && start == 0))
 242                 return false;
 243         }
 244     }
 245     double lat, lon;
 246     GeoEncode::decode(value, lat, lon);
 247     if (lat < min_lat || lat > max_lat) {
 248         return false;
 249     }
 250     if (lat == -90 || lat == 90) {
 251         // It's a pole, so the longitude isn't meaningful (will be zero)
 252         // and we've already checked that the latitude is in range.
 253         lat_ref = lat;
 254         lon_ref = 0;
 255         return true;
 256     }
 257     if (discontinuous_longitude_range) {
 258         if (lon2 < lon && lon < lon1)
 259             return false;
 260     } else {
 261         if (lon < lon1 || lon2 < lon)
 262             return false;
 263     }
 264
 265     lat_ref = lat;
 266     lon_ref = lon;
 267     return true;
 268 }