Opened 7 months ago

Closed 6 months ago

#5641 closed defect (invalid)

ST_TileEnvelope often has wrong XMax

Reported by: robe Owned by: pramsey
Priority: medium Milestone: PostGIS 3.4.2
Component: postgis Version: 3.4.x
Keywords: Cc:

Description

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;
}

Change History (1)

comment:1 by pramsey, 6 months ago

Resolution: invalid
Status: newclosed

I'm not convinced I'm seeing anything here other than an order of operations difference. The demo code divides by and then multiplies by tile_size, for example.

Note: See TracTickets for help on using tickets.