/**********************************************************************

 *

 * PostGIS - Spatial Types for PostgreSQL

 * http://postgis.net

 *

 * PostGIS is free software: you can redistribute it and/or modify

 * it under the terms of the GNU General Public License as published by

 * the Free Software Foundation, either version 2 of the License, or

 * (at your option) any later version.

 *

 * PostGIS is distributed in the hope that it will be useful,

 * but WITHOUT ANY WARRANTY; without even the implied warranty of

 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the

 * GNU General Public License for more details.

 *

 * You should have received a copy of the GNU General Public License

 * along with PostGIS.  If not, see <http://www.gnu.org/licenses/>.

 *

 **********************************************************************

 *

 * Copyright 2008 Paul Ramsey

 *

 **********************************************************************/

#include "liblwgeom_internal.h"

#include "lwgeom_log.h"

#include <ctype.h> /* for tolower */

#include <stdbool.h>

int

p4d_same(const POINT4D *p1, const POINT4D *p2)

{

  return p1->x == p2->x

      && p1->y == p2->y

      && p1->z == p2->z

      && p1->m == p2->m;

}

int

p3d_same(const POINT3D *p1, const POINT3D *p2)

{

  return p1->x == p2->x

      && p1->y == p2->y

      && p1->z == p2->z;

}

int

p3dz_same(const POINT3DZ *p1, const POINT3DZ *p2)

{

  return p1->x == p2->x

      && p1->y == p2->y

      && p1->z == p2->z;

}

int

p2d_same(const POINT2D *p1, const POINT2D *p2)

{

  return p1->x == p2->x

      && p1->y == p2->y;

}

/**

* lw_segment_side()

*

* Return -1  if point Q is left of segment P

* Return  1  if point Q is right of segment P

* Return  0  if point Q in on segment P

*/

int lw_segment_side(const POINT2D *p1, const POINT2D *p2, const POINT2D *q)

{

  double side = ( (q->x - p1->x) * (p2->y - p1->y) - (p2->x - p1->x) * (q->y - p1->y) );

  return SIGNUM(side);

}

/**

* Returns the length of a linear segment

*/

double

lw_seg_length(const POINT2D *A1, const POINT2D *A2)

{

  return sqrt((A1->x-A2->x)*(A1->x-A2->x)+(A1->y-A2->y)*(A1->y-A2->y));

}

/**

* Returns true if P is on the same side of the plane partition

* defined by A1/A3 as A2 is. Only makes sense if P has already been

* determined to be on the circle defined by A1/A2/A3.

*/

int

lw_pt_in_arc(const POINT2D *P, const POINT2D *A1, const POINT2D *A2, const POINT2D *A3)

{

  int pt_side = lw_segment_side(A1, A3, P);

  int arc_side = lw_segment_side(A1, A3, A2);

  return pt_side == 0 || pt_side == arc_side;

}

/**

* Returns true if P is between A1/A2. Only makes sense if P has already been

* determined to be on the line defined by A1/A2.

*/

int

lw_pt_in_seg(const POINT2D *P, const POINT2D *A1, const POINT2D *A2)

{

  return ((A1->x <= P->x && P->x < A2->x) || (A1->x >= P->x && P->x > A2->x)) ||

         ((A1->y <= P->y && P->y < A2->y) || (A1->y >= P->y && P->y > A2->y));

}

int

lw_pt_on_segment(const POINT2D* A1, const POINT2D* A2, const POINT2D* P)

{

  int side = lw_segment_side(A1, A2, P);

  if (side != 0) return LW_FALSE;

  return lw_pt_in_seg(P, A1, A2);

}

/**

* Returns true if arc A is actually a point (all vertices are the same) .

*/

int

lw_arc_is_pt(const POINT2D *A1, const POINT2D *A2, const POINT2D *A3)

{

  if ( A1->x == A2->x && A2->x == A3->x &&

       A1->y == A2->y && A2->y == A3->y )

    return LW_TRUE;

  else

    return LW_FALSE;

}

/**

* Returns the length of a circular arc segment

*/

double

lw_arc_length(const POINT2D *A1, const POINT2D *A2, const POINT2D *A3)

{

  POINT2D C;

  double radius_A, circumference_A;

  int a2_side, clockwise;

  double a1, a3;

  double angle;

  if ( lw_arc_is_pt(A1, A2, A3) )

    return 0.0;

  radius_A = lw_arc_center(A1, A2, A3, &C);

  /* Co-linear! Return linear distance! */

  if ( radius_A < 0 )

  {

        double dx = A1->x - A3->x;

        double dy = A1->y - A3->y;

    return sqrt(dx*dx + dy*dy);

  }

  /* Closed circle! Return the circumference! */

  circumference_A = M_PI * 2 * radius_A;

  if ( p2d_same(A1, A3) )

    return circumference_A;

  /* Determine the orientation of the arc */

  a2_side = lw_segment_side(A1, A3, A2);

  /* The side of the A1/A3 line that A2 falls on dictates the sweep

     direction from A1 to A3. */

  if ( a2_side == -1 )

    clockwise = LW_TRUE;

  else

    clockwise = LW_FALSE;

  /* Angles of each point that defines the arc section */

  a1 = atan2(A1->y - C.y, A1->x - C.x);

  a3 = atan2(A3->y - C.y, A3->x - C.x);

  /* What's the sweep from A1 to A3? */

  if ( clockwise )

  {

    if ( a1 > a3 )

      angle = a1 - a3;

    else

      angle = 2*M_PI + a1 - a3;

  }

  else

  {

    if ( a3 > a1 )

      angle = a3 - a1;

    else

      angle = 2*M_PI + a3 - a1;

  }

  /* Length as proportion of circumference */

  return circumference_A * (angle / (2*M_PI));

}

int lw_arc_side(const POINT2D *A1, const POINT2D *A2, const POINT2D *A3, const POINT2D *Q)

{

  POINT2D C;

  double radius_A;

  double side_Q, side_A2;

  double d;

  side_Q = lw_segment_side(A1, A3, Q);

  radius_A = lw_arc_center(A1, A2, A3, &C);

  side_A2 = lw_segment_side(A1, A3, A2);

  /* Linear case */

  if ( radius_A < 0 )

    return side_Q;

  d = distance2d_pt_pt(Q, &C);

  /* Q is on the arc boundary */

  if ( d == radius_A && side_Q == side_A2 )

  {

    return 0;

  }

  /* Q on A1-A3 line, so its on opposite side to A2 */

  if ( side_Q == 0 )

  {

    return -1 * side_A2;

  }

  /*

  * Q is inside the arc boundary, so it's not on the side we

  * might think from examining only the end points

  */

  if ( d < radius_A && side_Q == side_A2 )

  {

    side_Q *= -1;

  }

  return side_Q;

}

/**

* Determines the center of the circle defined by the three given points.

* In the event the circle is complete, the midpoint of the segment defined

* by the first and second points is returned.  If the points are collinear,

* as determined by equal slopes, then -1.0 is returned.  If the interior

* point is coincident with either end point, they are taken as collinear.

* For non-collinear cases, arc radious is returned.

*/

double

lw_arc_center(const POINT2D *p1, const POINT2D *p2, const POINT2D *p3, POINT2D *result)

{

  POINT2D c;

  double cx, cy, cr;

  double dx21, dy21, dx31, dy31, h21, h31, d;

  c.x = c.y = 0.0;

  LWDEBUGF(2, "lw_arc_center called (%.16f,%.16f), (%.16f,%.16f), (%.16f,%.16f).", p1->x, p1->y, p2->x, p2->y, p3->x, p3->y);

  /* Closed circle */

  if (fabs(p1->x - p3->x) < EPSILON_SQLMM &&

      fabs(p1->y - p3->y) < EPSILON_SQLMM)

  {

    cx = p1->x + (p2->x - p1->x) / 2.0;

    cy = p1->y + (p2->y - p1->y) / 2.0;

    c.x = cx;

    c.y = cy;

    *result = c;

    cr = sqrt(pow(cx - p1->x, 2.0) + pow(cy - p1->y, 2.0));

    return cr;

  }

  /* Using cartesian eguations from page https://en.wikipedia.org/wiki/Circumscribed_circle */

  dx21 = p2->x - p1->x;

  dy21 = p2->y - p1->y;

  dx31 = p3->x - p1->x;

  dy31 = p3->y - p1->y;

  h21 = pow(dx21, 2.0) + pow(dy21, 2.0);

  h31 = pow(dx31, 2.0) + pow(dy31, 2.0);

  /* 2 * |Cross product|, d<0 means clockwise and d>0 counterclockwise sweeping angle */

  d = 2 * (dx21 * dy31 - dx31 * dy21);

  /* Check colinearity, |Cross product| = 0 */

  if (fabs(d) < EPSILON_SQLMM)

    return -1.0;

  /* Calculate centroid coordinates and radius */

  cx = p1->x + (h21 * dy31 - h31 * dy21) / d;

  cy = p1->y - (h21 * dx31 - h31 * dx21) / d;

  c.x = cx;

  c.y = cy;

  *result = c;

  cr = sqrt(pow(cx - p1->x, 2) + pow(cy - p1->y, 2));

  LWDEBUGF(2, "lw_arc_center center is (%.16f,%.16f)", result->x, result->y);

  return cr;

}

static int

lw_seg_interact(const POINT2D *p1, const POINT2D *p2, const POINT2D *q1, const POINT2D *q2)

{

  double minq=FP_MIN(q1->x,q2->x);

  double maxq=FP_MAX(q1->x,q2->x);

  double minp=FP_MIN(p1->x,p2->x);

  double maxp=FP_MAX(p1->x,p2->x);

  if (FP_GT(minp,maxq) || FP_LT(maxp,minq))

    return LW_FALSE;

  minq=FP_MIN(q1->y,q2->y);

  maxq=FP_MAX(q1->y,q2->y);

  minp=FP_MIN(p1->y,p2->y);

  maxp=FP_MAX(p1->y,p2->y);

  if (FP_GT(minp,maxq) || FP_LT(maxp,minq))

    return LW_FALSE;

  return LW_TRUE;

}

/**

** @brief returns the kind of #CG_SEGMENT_INTERSECTION_TYPE  behavior of lineseg 1 (constructed from p1 and p2) and lineseg 2 (constructed from q1 and q2)

**  @param p1 start point of first straight linesegment

**  @param p2 end point of first straight linesegment

**  @param q1 start point of second line segment

**  @param q2 end point of second line segment

**  @return a #CG_SEGMENT_INTERSECTION_TYPE

**  Returns one of

**    SEG_ERROR = -1,

**    SEG_NO_INTERSECTION = 0,

**    SEG_COLINEAR = 1,

**    SEG_CROSS_LEFT = 2,

**    SEG_CROSS_RIGHT = 3,

*/

int lw_segment_intersects(const POINT2D *p1, const POINT2D *p2, const POINT2D *q1, const POINT2D *q2)

{

  int pq1, pq2, qp1, qp2;

  /* No envelope interaction => we are done. */

  if (!lw_seg_interact(p1, p2, q1, p2))

  {

    return SEG_NO_INTERSECTION;

  }

  /* Are the start and end points of q on the same side of p? */

  pq1=lw_segment_side(p1,p2,q1);

  pq2=lw_segment_side(p1,p2,q2);

  if ((pq1>0 && pq2>0) || (pq1<0 && pq2<0))

  {

    return SEG_NO_INTERSECTION;

  }

  /* Are the start and end points of p on the same side of q? */

  qp1=lw_segment_side(q1,q2,p1);

  qp2=lw_segment_side(q1,q2,p2);

  if ( (qp1 > 0.0 && qp2 > 0.0) || (qp1 < 0.0 && qp2 < 0.0) )

  {

    return SEG_NO_INTERSECTION;

  }

  /* Nobody is on one side or another? Must be colinear. */

  if ( pq1 == 0.0 && pq2 == 0.0 && qp1 == 0.0 && qp2 == 0.0 )

  {

    return SEG_COLINEAR;

  }

  /*

  ** When one end-point touches, the sidedness is determined by the

  ** location of the other end-point. Only touches by the first point

  ** will be considered "real" to avoid double counting.

  */

  LWDEBUGF(4, "pq1=%d pq2=%d", pq1, pq2);

  LWDEBUGF(4, "qp1=%d qp2=%d", qp1, qp2);

  /* Second point of p or q touches, it's not a crossing. */

  if ( pq2 == 0 || qp2 == 0 )

  {

    return SEG_NO_INTERSECTION;

  }

  /* First point of p touches, it's a "crossing". */

  if ( pq1 == 0 )

  {

    if ( pq2 > 0 )

      return SEG_CROSS_RIGHT;

    else

      return SEG_CROSS_LEFT;

  }

  /* First point of q touches, it's a crossing. */

  if ( qp1 == 0 )

  {

    if ( pq1 < pq2 )

      return SEG_CROSS_RIGHT;

    else

      return SEG_CROSS_LEFT;

  }

  /* The segments cross, what direction is the crossing? */

  if ( pq1 < pq2 )

    return SEG_CROSS_RIGHT;

  else

    return SEG_CROSS_LEFT;

  /* This should never happen! */

  return SEG_ERROR;

}

/**

** @brief lwline_crossing_direction: returns the kind of #CG_LINE_CROSS_TYPE behavior  of 2 linestrings

** @param l1 first line string

** @param l2 second line string

** @return a #CG_LINE_CROSS_TYPE

**   LINE_NO_CROSS = 0

**   LINE_CROSS_LEFT = -1

**   LINE_CROSS_RIGHT = 1

**   LINE_MULTICROSS_END_LEFT = -2

**   LINE_MULTICROSS_END_RIGHT = 2

**   LINE_MULTICROSS_END_SAME_FIRST_LEFT = -3

**   LINE_MULTICROSS_END_SAME_FIRST_RIGHT = 3

**

*/

int lwline_crossing_direction(const LWLINE *l1, const LWLINE *l2)

{

  uint32_t i = 0, j = 0;

  const POINT2D *p1, *p2, *q1, *q2;

  POINTARRAY *pa1 = NULL, *pa2 = NULL;

  int cross_left = 0;

  int cross_right = 0;

  int first_cross = 0;

  int this_cross = 0;

#if POSTGIS_DEBUG_LEVEL >= 4

  char *geom_ewkt;

#endif

  pa1 = (POINTARRAY*)l1->points;

  pa2 = (POINTARRAY*)l2->points;

  /* One-point lines can't intersect (and shouldn't exist). */

  if ( pa1->npoints < 2 || pa2->npoints < 2 )

    return LINE_NO_CROSS;

  /* Zero length lines don't have a side. */

  if ( ptarray_length_2d(pa1) == 0 || ptarray_length_2d(pa2) == 0 )

    return LINE_NO_CROSS;

#if POSTGIS_DEBUG_LEVEL >= 4

  geom_ewkt = lwgeom_to_ewkt((LWGEOM*)l1);

  LWDEBUGF(4, "l1 = %s", geom_ewkt);

  lwfree(geom_ewkt);

  geom_ewkt = lwgeom_to_ewkt((LWGEOM*)l2);

  LWDEBUGF(4, "l2 = %s", geom_ewkt);

  lwfree(geom_ewkt);

#endif

  /* Initialize first point of q */

  q1 = getPoint2d_cp(pa2, 0);

  for ( i = 1; i < pa2->npoints; i++ )

  {

    /* Update second point of q to next value */

    q2 = getPoint2d_cp(pa2, i);

    /* Initialize first point of p */

    p1 = getPoint2d_cp(pa1, 0);

    for ( j = 1; j < pa1->npoints; j++ )

    {

      /* Update second point of p to next value */

      p2 = getPoint2d_cp(pa1, j);

      this_cross = lw_segment_intersects(p1, p2, q1, q2);

      LWDEBUGF(4, "i=%d, j=%d (%.8g %.8g, %.8g %.8g)", i, j, p1->x, p1->y, p2->x, p2->y);

      if ( this_cross == SEG_CROSS_LEFT )

      {

        LWDEBUG(4,"this_cross == SEG_CROSS_LEFT");

        cross_left++;

        if ( ! first_cross )

          first_cross = SEG_CROSS_LEFT;

      }

      if ( this_cross == SEG_CROSS_RIGHT )

      {

        LWDEBUG(4,"this_cross == SEG_CROSS_RIGHT");

        cross_right++;

        if ( ! first_cross )

          first_cross = SEG_CROSS_RIGHT;

      }

      /*

      ** Crossing at a co-linearity can be turned handled by extending

      ** segment to next vertex and seeing if the end points straddle

      ** the co-linear segment.

      */

      if ( this_cross == SEG_COLINEAR )

      {

        LWDEBUG(4,"this_cross == SEG_COLINEAR");

        /* TODO: Add logic here and in segment_intersects()

        continue;

        */

      }

      LWDEBUG(4,"this_cross == SEG_NO_INTERSECTION");

      /* Turn second point of p into first point */

      p1 = p2;

    }

    /* Turn second point of q into first point */

    q1 = q2;

  }

  LWDEBUGF(4, "first_cross=%d, cross_left=%d, cross_right=%d", first_cross, cross_left, cross_right);

  if ( !cross_left && !cross_right )

    return LINE_NO_CROSS;

  if ( !cross_left && cross_right == 1 )

    return LINE_CROSS_RIGHT;

  if ( !cross_right && cross_left == 1 )

    return LINE_CROSS_LEFT;

  if ( cross_left - cross_right == 1 )

    return LINE_MULTICROSS_END_LEFT;

  if ( cross_left - cross_right == -1 )

    return LINE_MULTICROSS_END_RIGHT;

  if ( cross_left - cross_right == 0 && first_cross == SEG_CROSS_LEFT )

    return LINE_MULTICROSS_END_SAME_FIRST_LEFT;

  if ( cross_left - cross_right == 0 && first_cross == SEG_CROSS_RIGHT )

    return LINE_MULTICROSS_END_SAME_FIRST_RIGHT;

  return LINE_NO_CROSS;

}

static char *base32 = "0123456789bcdefghjkmnpqrstuvwxyz";

#define GEOHASH_MAX_DOUBLE_PRECISION_CHARS 20

/*

** Calculate the geohash, iterating downwards and gaining precision.

** From geohash-native.c, (c) 2008 David Troy <dave@roundhousetech.com>

** Released under the MIT License.

*/

lwvarlena_t *

geohash_point(double longitude, double latitude, int precision)

{

  int is_even=1, i=0;

  double lat[2], lon[2], mid;

  char bits[] = {16,8,4,2,1};

  int bit=0, ch=0;

  lwvarlena_t *v = lwalloc(precision + LWVARHDRSZ);

  LWSIZE_SET(v->size, precision + LWVARHDRSZ);

  char *geohash = v->data;

  lat[0] = -90.0;

  lat[1] = 90.0;

  lon[0] = -180.0;

  lon[1] = 180.0;

  while (i < precision)

  {

    if (is_even)

    {

      mid = (lon[0] + lon[1]) / 2;

      if (longitude >= mid)

      {

        ch |= bits[bit];

        lon[0] = mid;

      }

      else

      {

        lon[1] = mid;

      }

    }

    else

    {

      mid = (lat[0] + lat[1]) / 2;

      if (latitude >= mid)

      {

        ch |= bits[bit];

        lat[0] = mid;

      }

      else

      {

        lat[1] = mid;

      }

    }

    is_even = !is_even;

    if (bit < 4)

    {

      bit++;

    }

    else

    {

      geohash[i++] = base32[ch];

      bit = 0;

      ch = 0;

    }

  }

  return v;

}

/*

** Calculate the geohash, iterating downwards and gaining precision.

** From geohash-native.c, (c) 2008 David Troy <dave@roundhousetech.com>

** Released under the MIT License.

*/

unsigned int geohash_point_as_int(POINT2D *pt)

{

  int is_even=1;

  double lat[2], lon[2], mid;

  int bit=32;

  unsigned int ch = 0;

  double longitude = pt->x;

  double latitude = pt->y;

  lat[0] = -90.0;

  lat[1] = 90.0;

  lon[0] = -180.0;

  lon[1] = 180.0;

  while (--bit >= 0)

  {

    if (is_even)

    {

      mid = (lon[0] + lon[1]) / 2;

      if (longitude > mid)

      {

        ch |= 0x0001u << bit;

        lon[0] = mid;

      }

      else

      {

        lon[1] = mid;

      }

    }

    else

    {

      mid = (lat[0] + lat[1]) / 2;

      if (latitude > mid)

      {

        ch |= 0x0001 << bit;

        lat[0] = mid;

      }

      else

      {

        lat[1] = mid;

      }

    }

    is_even = !is_even;

  }

  return ch;

}

/*

** Decode a GeoHash into a bounding box. The lat and lon arguments should

** both be passed as double arrays of length 2 at a minimum where the values

** set in them will be the southwest and northeast coordinates of the bounding

** box accordingly. A precision less than 0 indicates that the entire length

** of the GeoHash should be used.

** It will call `lwerror` if an invalid character is found

*/

void decode_geohash_bbox(char *geohash, double *lat, double *lon, int precision)

{

  bool is_even = 1;

  lat[0] = -90.0;

  lat[1] = 90.0;

  lon[0] = -180.0;

  lon[1] = 180.0;

  size_t hashlen = strlen(geohash);

  if (precision < 0 || (size_t)precision > hashlen)

  {

    precision = (int)hashlen;

  }

  for (int i = 0; i < precision; i++)

  {

    char c = tolower(geohash[i]);

    /* Valid characters are all digits in base32 */

    char *base32_pos = strchr(base32, c);

    if (!base32_pos)

    {

      lwerror("%s: Invalid character '%c'", __func__, geohash[i]);

      return;

    }

    char cd = base32_pos - base32;

    for (size_t j = 0; j < 5; j++)

    {

      const char bits[] = {16, 8, 4, 2, 1};

      char mask = bits[j];

      if (is_even)

      {

        lon[!(cd & mask)] = (lon[0] + lon[1]) / 2;

      }

      else

      {

        lat[!(cd & mask)] = (lat[0] + lat[1]) / 2;

      }

      is_even = !is_even;

    }

  }

}

int lwgeom_geohash_precision(GBOX bbox, GBOX *bounds)

{

  double minx, miny, maxx, maxy;

  double latmax, latmin, lonmax, lonmin;

  double lonwidth, latwidth;

  double latmaxadjust, lonmaxadjust, latminadjust, lonminadjust;

  int precision = 0;

  /* Get the bounding box, return error if things don't work out. */

  minx = bbox.xmin;

  miny = bbox.ymin;

  maxx = bbox.xmax;

  maxy = bbox.ymax;

  if ( minx == maxx && miny == maxy )

  {

    /* It's a point. Doubles have 51 bits of precision.

    ** 2 * 51 / 5 == 20 */

    return GEOHASH_MAX_DOUBLE_PRECISION_CHARS;

  }

  lonmin = -180.0;

  latmin = -90.0;

  lonmax = 180.0;

  latmax = 90.0;

  /* Shrink a world bounding box until one of the edges interferes with the

  ** bounds of our rectangle. */

  while ( 1 )

  {

    double old_lonmin = lonmin;

    double old_lonmax = lonmax;

    double old_latmin = latmin;

    double old_latmax = latmax;

    lonwidth = lonmax - lonmin;

    latwidth = latmax - latmin;

    latmaxadjust = lonmaxadjust = latminadjust = lonminadjust = 0.0;

    if ( minx >= lonmin + lonwidth / 2.0 )

    {

      lonminadjust = lonwidth / 2.0;

    }

    else if ( maxx <= lonmax - lonwidth / 2.0 )

    {

      lonmaxadjust = -1 * lonwidth / 2.0;

    }

    if ( lonminadjust || lonmaxadjust )

    {

      lonmin += lonminadjust;

      lonmax += lonmaxadjust;

      /* Stop before counting bits that no longer move the double bounds. */

      if (lonmin == old_lonmin && lonmax == old_lonmax)

        break;

      /* Each adjustment cycle corresponds to 2 bits of storage in the

      ** geohash. */

      precision++;

    }

    else

    {

      break;

    }

    if ( miny >= latmin + latwidth / 2.0 )

    {

      latminadjust = latwidth / 2.0;

    }

    else if (maxy <= latmax - latwidth / 2.0 )

    {

      latmaxadjust = -1 * latwidth / 2.0;

    }

    /* Only adjust if adjustments are legal (we haven't crossed any edges). */

    if ( latminadjust || latmaxadjust )

    {

      latmin += latminadjust;

      latmax += latmaxadjust;

      /* Stop before counting bits that no longer move the double bounds. */

      if (latmin == old_latmin && latmax == old_latmax)

        break;

      /* Each adjustment cycle corresponds to 2 bits of storage in the

      ** geohash. */

      precision++;

    }

    else

    {

      break;

    }

  }

  /* Save the edges of our bounds, in case someone cares later. */

  bounds->xmin = lonmin;

  bounds->xmax = lonmax;

  bounds->ymin = latmin;

  bounds->ymax = latmax;

  /* Each geohash character (base32) can contain 5 bits of information.

  ** We are returning the precision in characters, so here we divide. */

  precision /= 5;

  return precision > GEOHASH_MAX_DOUBLE_PRECISION_CHARS ? GEOHASH_MAX_DOUBLE_PRECISION_CHARS : precision;

}

/*

** Return a geohash string for the geometry. <http://geohash.org>

** Where the precision is non-positive, calculate a precision based on the

** bounds of the feature. Big features have loose precision.

** Small features have tight precision.

*/

lwvarlena_t *

lwgeom_geohash(const LWGEOM *lwgeom, int precision)

{

  GBOX gbox = {0};

  GBOX gbox_bounds = {0};

  double lat, lon;

  int result;

  gbox_init(&gbox);

  gbox_init(&gbox_bounds);

  result = lwgeom_calculate_gbox_cartesian(lwgeom, &gbox);

  if ( result == LW_FAILURE ) return NULL;

  /* Return error if we are being fed something outside our working bounds */

  if ( gbox.xmin < -180 || gbox.ymin < -90 || gbox.xmax > 180 || gbox.ymax > 90 )

  {

    lwerror("Geohash requires inputs in decimal degrees, got (%g %g, %g %g).",

       gbox.xmin, gbox.ymin,

       gbox.xmax, gbox.ymax);

    return NULL;

  }

  /* What is the center of our geometry bounds? We'll use that to

  ** approximate location. */

  lon = gbox.xmin + (gbox.xmax - gbox.xmin) / 2;

  lat = gbox.ymin + (gbox.ymax - gbox.ymin) / 2;

  if ( precision <= 0 )

  {

    precision = lwgeom_geohash_precision(gbox, &gbox_bounds);

  }

  /*

  ** Return the geohash of the center, with a precision determined by the

  ** extent of the bounds.

  ** Possible change: return the point at the center of the precision bounds?

  */

  return geohash_point(lon, lat, precision);

}