[PostGIS] #5641: ST_TileEnvelope often has wrong XMax

[PostGIS]

#5641: ST_TileEnvelope often has wrong XMax
---------------------+---------------------------
 Reporter:  robe     |      Owner:  pramsey
     Type:  defect   |     Status:  new
 Priority:  medium   |  Milestone:  PostGIS 3.4.2
Component:  postgis  |    Version:  3.4.x
 Keywords:           |
---------------------+---------------------------
 Per velix on matrix/IRC


 {{{
 SELECT ST_TileEnvelope(15, 17017, 10755)
 }}}

 should be:


 {{{
 774154.2224722654
 6883001.523023552
 775377.2149248281
 6884224.515476115
 }}}

 But is :


 {{{
 ┌───────────────────┬───────────────────┬───────────────────┬───────────────────┐
 │      st_xmin      │      st_ymin      │      st_xmax      │      st_ymax
 │
 ├───────────────────┼───────────────────┼───────────────────┼───────────────────┤
 │ 774154.2224722654 │ 6883001.523023551 │ 775377.2149248272 │
 6884224.515476115 │
 └───────────────────┴───────────────────┴───────────────────┴───────────────────┘


 }}}

 Based on plain c code should be computing like this:


 {{{
 #include <stdio.h>
 #include <math.h>

 #define SEMIMAJOR_AXIS 6378137
 #define TILE_SIZE 256
 #define LONGER_PI
 3.1415926535897932384626433832795028841971693993751068316706821L

 void get_tile_envelope(int zoom, int x, int y, long double *envelope) {
     long double earth_circumference = 2.0L * LONGER_PI * SEMIMAJOR_AXIS;
     long double initial_resolution = earth_circumference / TILE_SIZE;

     // Calculate the resolution for the given zoom level
     long double resolution = initial_resolution / powl(2.0L, zoom);

     // Calculate the tile's bounding box coordinates
     long double min_x = -earth_circumference / 2.0L + x * TILE_SIZE *
 resolution;
     long double max_y = +earth_circumference / 2.0L - y * TILE_SIZE *
 resolution;
     long double max_x = min_x + TILE_SIZE * resolution;
     long double min_y = max_y - TILE_SIZE * resolution;

     envelope[0] = min_x;
     envelope[1] = min_y;
     envelope[2] = max_x;
     envelope[3] = max_y;
 }

 int main() {
     long double tile_envelope[4];
     get_tile_envelope(15, 17108, 11130, tile_envelope);

     printf("%.15Lf\n", tile_envelope[0]);
     printf("%.15Lf\n", tile_envelope[1]);
     printf("%.15Lf\n", tile_envelope[2]);
     printf("%.15Lf\n", tile_envelope[3]);

     return 0;
 }
 }}}
-- 
Ticket URL: <https://trac.osgeo.org/postgis/ticket/5641>
PostGIS <http://trac.osgeo.org/postgis/>
The PostGIS Trac is used for bug, enhancement & task tracking, a user and developer wiki, and a view into the subversion code repository of PostGIS project.