Skip to content

T-Rex doesn't handle invalid geometry results correctly when transforming between coordinate systems #315

Description

@hamarituc

I encountered an issue where invalid point geometries are served.

Steps to reproduce

  1. Consider the following data in a PostGIS database
BEGIN;
CREATE TABLE "public"."test" ( "ogc_fid" SERIAL, CONSTRAINT "test_pk" PRIMARY KEY ("ogc_fid") );
SELECT AddGeometryColumn('public','test','wkb_geometry',25833,'POLYGON',2);
CREATE INDEX "test_wkb_geometry_geom_idx" ON "public"."test" USING GIST ("wkb_geometry");
INSERT INTO "public"."test" ("wkb_geometry" , "ogc_fid" ) VALUES ('0103000020E964000001000000050000009AA6C935869C154110164AAEEB7D55410E92AED1A09C15410E9FB145EE7D554178EADB7FCD9C15413B2A246DEB7D55418E44E7E2989C15411D3A977BEA7D55419AA6C935869C154110164AAEEB7D5541', 1);
INSERT INTO "public"."test" ("wkb_geometry" , "ogc_fid" ) VALUES ('0103000020E9640000010000000500000092E008F4179D1541534A8CE8EF7D5541BA518FE92F9D1541A47BCA5FED7D5541341B6238FA9C154171581624EC7D5541C3CA4C05DF9C15413675F754F07D554192E008F4179D1541534A8CE8EF7D5541', 2);
INSERT INTO "public"."test" ("wkb_geometry" , "ogc_fid" ) VALUES ('0103000020E96400000100000006000000023A12A91A9D154167CE46A0F37D55413236D5033A9D15413FFE54E8F27D5541DCBE1A20359D1541CABD12D9F17D55413E4F782A419D1541078DD709F07D554118FCCDC60E9D154136BDE349F07D5541023A12A91A9D154167CE46A0F37D5541', 3);
INSERT INTO "public"."test" ("wkb_geometry" , "ogc_fid" ) VALUES ('0103000020E9640000010000000400000023173BB3B89C1541FD196035EA7D55415954A095C79C1541A63D3AA9E87D5541AB142256D79C154159990D34EA7D554123173BB3B89C1541FD196035EA7D5541', 4);
INSERT INTO "public"."test" ("wkb_geometry" , "ogc_fid" ) VALUES ('0103000020E964000001000000050000002EE515095A9C15417FC68947F47D5541704A76218B9C1541618A4738F77D5541DFE97506CC9C1541E1C4BEF9F17D5541133529E45A9C15413048178FF17D55412EE515095A9C15417FC68947F47D5541', 5);
INSERT INTO "public"."test" ("wkb_geometry" , "ogc_fid" ) VALUES ('0103000020E964000001000000050000005C84D2275C9C15418E1EA353F07D554176BCCBA6819C15411D415EFCEF7D55410C1A04257F9C1541F43B6077EE7D55414C1200FB4C9C1541E9A6B3E9EE7D55415C84D2275C9C15418E1EA353F07D5541', 6);
INSERT INTO "public"."test" ("wkb_geometry" , "ogc_fid" ) VALUES ('0103000020E9640000010000000500000055487DCD589C154156222501EE7D5541C714C0706E9C154184D3A3C3ED7D554156540F1E6C9C1541D7E2667EEC7D55412966C7475C9C15415A164C7FEC7D554155487DCD589C154156222501EE7D5541', 7);
INSERT INTO "public"."test" ("wkb_geometry" , "ogc_fid" ) VALUES ('0103000020E96400000100000005000000147B7306D69C154124585008F67D55413EC200D2D99C1541356B1CEBF27D55413D5F9FE8FC9C1541397C7E6BF27D554106D8B1FEFA9C1541043E24ADF57D5541147B7306D69C154124585008F67D5541', 8);
COMMIT;
  1. Configure T-Rex to serve polygon centroids from this database
[service.mvt]
viewer = false

[[datasource]]
name = "test"
dbconn = "as required"

[grid]
predefined = "web_mercator"

[[tileset]]
name = "test"
#minzoom = 15
#maxzoom = 22
extent = { minx = 12.9053, maxx = 12.9609, miny = 50.7920, maxy = 50.8538 }
center = [ 12.92804, 50.83904 ]
start_zoom = 17

[[tileset.layer]]
name = "centroid"
datasource = "test"
geometry_field = "geom"
geometry_type = "POINT"
srid = 25833
fid_field = "ogc_fid"
# Clip polygons with a buffer
buffer_size = 200
simplify = false

  [[tileset.layer.query]]
  minzoom = 17
  maxzoom = 24
  sql = "SELECT ogc_fid, ST_PointOnSurface(wkb_geometry) AS geom FROM test WHERE ST_PointOnSurface(wkb_geometry) && !bbox!"

[cache.file]
base = "/tmp/trex/cache"

[webserver]
bind = "localhost"
port = 6767
  1. Fetch a tile from /vector/test/19/280971/175891.pbf

Observed Behavior

Loading 175891.pbf into QGis reveals there is a single point at the lower left corner of the tile not present in the database:

grafik

Expected Behavior

The served tile should not contain invalid geometries.

I suppose this invalid geometries arise when the databases uses a different coordinate reference system than Web mercator. According to the docs the &&-operator should be used inside the WHERE-clause to filter for all relevant geometries for the current tile. This operator compares the bounding boxes of both of its operands. Later on T-Rex clips the result of the query by means of the ST_Intersection()-function to the exact tile extent.

SELECT
  ST_Transform(
    ST_Intersection(
      geom,
      ST_Transform(
        ST_Segmentize(
          ST_MakeEnvelope(
            1439079.9315250143-12.5*0.2985821417389698::FLOAT8,
            6592846.563653022-12.5*0.2985821417389698::FLOAT8,
            1439156.3685532995+12.5*0.2985821417389698::FLOAT8,
            6592923.000681307+12.5*0.2985821417389698::FLOAT8,
            3857
          ),
          (1439156.3685532995-1439079.9315250143)/512
        ),
        25833
      )
    ),
    3857
  )
AS geom,"ogc_fid"
FROM
(
  SELECT ogc_fid, ST_PointOnSurface(wkb_geometry) AS geom
  FROM test
  WHERE ST_PointOnSurface(wkb_geometry) &&
    ST_Transform(
      ST_Segmentize(
        ST_MakeEnvelope(
          1439079.9315250143-12.5*0.2985821417389698::FLOAT8,
          6592846.563653022-12.5*0.2985821417389698::FLOAT8,
          1439156.3685532995+12.5*0.2985821417389698::FLOAT8,
          6592923.000681307+12.5*0.2985821417389698::FLOAT8,
          3857
        ),
        (1439156.3685532995-1439079.9315250143)/512
      ),
      25833
    )
) AS _q

The tile area (including a buffer, if configured) itself is transformed into the database native CRS (light-red polygon). The dark-red box shows the bounding box of the tile area. The green centroid therefore fulfills the WHERE-clause, but is subsequently clipped in the SELECT clause which yields an POINT EMPTY value for this row.

grafik

T-Rex doesn't handle that case correctly. It should ignore empty point geometries in this case.

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Type

    No type

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions