Missende BGT records na BGT Lean SQL script generatie

[justb4]

Bijv Wegdeel voor geometrie_vlak in: lokaalid = 'G0307.f86336b663854a3b9d00dd4d89a49b56' (fietspad in Amersfoort).

BGT Lean verwerpt dit record omdat ST_IsValid(geometrie_vlak) == False van origineel BGT Record uit wegdeelactueelbestaand VIEW in bgt.latest schema. Misschien door toleranties. In ieder geval is het een CURVE geometrie, als CSV WKT aangehecht. Deze kan ok in QGIS gesleept (in EPSG:28992 project).

G0307_f86336b663854a3b9d00dd4d89a49b56.csv

Dus is raar. Dit is wat bgt-lean.sql doet voor Wegdeel:

BEGIN;
-- Helper VIEW to select columns for final tables once
DROP VIEW IF EXISTS bgt_lean.wegdeel_cols CASCADE;
CREATE VIEW bgt_lean.wegdeel_cols AS
SELECT
    ogc_fid,
    lokaalid,
    relatievehoogteligging,
    bgt_functie,
    plus_functie,
    bgt_fysiekvoorkomen,
    plus_fysiekvoorkomen,
    wegdeeloptalud
FROM latest.wegdeelactueelbestaand;

DROP TABLE IF EXISTS bgt_lean.wegdeel_vlak CASCADE;
CREATE TABLE bgt_lean.wegdeel_vlak AS
SELECT
    a.*,
    ST_MakeValid(ST_CurveToLine(b.geometrie_vlak))::geometry(Polygon, 28992) as geometrie_vlak
FROM bgt_lean.wegdeel_cols a, latest.wegdeelactueelbestaand b
WHERE a.ogc_fid = b.ogc_fid AND b.geometrie_vlak is not null AND ST_IsValid(b.geometrie_vlak);
CREATE INDEX wegdeel_vlak_geometrie_vlak_geom_idx ON bgt_lean.wegdeel_vlak USING gist(geometrie_vlak);
-- CREATE INDEX wegdeel_vlak_lokaalid_idx ON bgt_lean.wegdeel_vlak USING btree (lokaalid);
ALTER TABLE ONLY bgt_lean.wegdeel_vlak ADD CONSTRAINT wegdeel_vlak_pkey PRIMARY KEY (ogc_fid);

DROP TABLE IF EXISTS bgt_lean.wegdeel_kruinlijn CASCADE;
CREATE TABLE bgt_lean.wegdeel_kruinlijn AS
SELECT
    a.*,
    ST_MakeValid(ST_CurveToLine(geometrie_kruinlijn))::geometry(LineString, 28992) as geometrie_kruinlijn
FROM bgt_lean.wegdeel_cols a, latest.wegdeelactueelbestaand b
WHERE a.ogc_fid = b.ogc_fid AND b.geometrie_kruinlijn is not null AND ST_IsValid(b.geometrie_kruinlijn);
CREATE INDEX wegdeel_kruinlijn_geometrie_kruinlijn_geom_idx ON bgt_lean.wegdeel_kruinlijn USING gist(geometrie_kruinlijn);
-- CREATE INDEX wegdeel_kruinlijn_lokaalid_idx ON bgt_lean.wegdeel_kruinlijn USING btree (lokaalid);
ALTER TABLE ONLY bgt_lean.wegdeel_kruinlijn ADD CONSTRAINT wegdeel_kruinlijn_pkey PRIMARY KEY (ogc_fid);

DROP VIEW bgt_lean.wegdeel_cols;
COMMIT;

Misschien maar ST_IsValid() conditie weglaten. En ST_MakeValid(ST_CurveToLine(b.geometrie_vlak))::geometry(Polygon, 28992) as geometrie_vlak -deel anders want levert een MULTIPOLYGON en ook ST_IsValid() False...

[justb4]

Uiteindelijk zijn er 851 records met invalid geometrie_vlak in Wegdeelactueelbestaand, althans dat komt uit deze query:

SELECT count(geometrie_vlak) FROM latest.wegdeelactueelbestaand where st_isvalid(geometrie_vlak) is false;

Allemaal 'self-intersecting' Polygons m.i.

Dus relatief heel weinig op 7.7 miljoen Wegdeel records. Misschien meldingen maken? Of is het toch iets in PostGIS (de GEOS lib). Er kan ook geen vlag of tolerantie gezet.

[justb4]

Voor Wegdeel nu: ST_CurveToLine(b.geometrie_vlak)::geometry(Polygon, 28992) dus zonder ST_IsValid() selectie van origineel en zonder ST_MakeValid() op ST_CurveToLine() resultaat. De laatste maakt er een MULTIPOLYGON van, wil je ook niet. Nu komt bewuste Fietspad wegdeel wel goed door, althans QGIS en Mapnik klagen niet:

[justb4]

Melding BGT00140630 gemaakt: https://www.verbeterdekaart.nl/#?geometry.x=154987.2&geometry.y=463328.39&zoomlevel=16

[willemhoffmans]

Ha Just, je kan er ook voor kiezen om ST_Buffer te gebruiken i.p.v. ST_MakeValid.

Dus: ST_Buffer(ST_CurveToLine(geometrie_vlak),0)

Dit maakt de geometrie wel valide, maar maakt er géén multipolygon van.

[justb4]

Bedankt @willemhoffmans . Ik was net bezig uitproberen tolerantie parameters van ST_CurveToLine(). Zo wordt resultaat ook een valide POLYGON:

select ST_CurveToLine(geometrie_vlak, 3, 1, 1) from latest.wegdeelactueelbestaand where lokaalid = 'G0307.f86336b663854a3b9d00dd4d89a49b56'

Tolerantie getal maakt niet eens uit 3 of 10 of 20, etc .Lastig in te schatten voor mij waarom nu wel goed...De ST_Buffer(0) methode lijkt simpeler.

[justb4]

Als ik bovenstaande toepas op alle vlakken in hele BGT zie ik dat zowel ST_Buffer(0) als tolerantie parameters in ST_CurveToLine() probleemgevallen opleveren:

• ST_Buffer(0) levert, lijkt het, altijd een valide geometrie, maar vaak een MultiPolygon waar Polygon gewenst, bijv in kunstwerkdeel vlak.
• er zijn geen ST_CurveToLine() tolerantie-parameters die in alle gevallen een valide geometrie opleveren

De vraag is dus: hoe kunnen we altijd een valide geometrie opleveren? Of dan maar MultiPolygons accepteren na ST_Buffer(0)?

Er is zelfs een stuk wetenschap hiervoor, door, en met literatuurverwijzingen naar, niet de minste auteurs: https://r-spatial.org/r/2017/03/19/invalid.html ...

Ik zie wel iets in een oplossing via een dedicated PG Function met daarin heuristiek om enkelvlakken POLYGON te maken. Mogelijk alleen voor de probleemgevallen ivm performance.

[justb4]

Voor nu misschien niet de meest elegante oplossing in commit 43a11e6, maar wel een oplossing die de polygonen zoveel mogelijk behoudt/repareert. Alle BGT vlakken die eerder Polygon waren zijn nu MultiPolygon. Voor toepassingen zoals mapping/kaarten en ruimtelijke analyses moet dit niet uitmaken en missen tenminste geen elementen. De typische manipulatie van een geometrie_vlak is nu:

ST_Multi(ST_Buffer(ST_CurveToLine(b.geometrie_vlak), 0))::geometry(MultiPolygon, 28992) as geometrie_vlak

Rationale:

• ST_Buffer(0) geeft regelmatig een MULTIPOLYGON. Dan maken we alle vlakken maar van dit type, middels ST_Multi().
• Geen foutmeldingen. Ook niet voor de cast ::geometry(MultiPolygon, 28992) (Zou m.i. zelfs overbodig zijn...)
• alle resulterende geometrieen zijn valide (ST_IsValid() True)
• belangrijkst: er gaat geen geometrie-vlak of deel daarvan verloren

ST_Buffer(0) +St_Multi() is m.i. in lijn eerdere oplossingen van @willemhoffmans (?).

[janwinsemius]

Ik ben toevallig zelf bgt aan het parsen, doe dat direct uit gml in rust. Dus daar kom ik ook die curves tegen. Ik schrijf de segmentering dus eigenlijk zelf en gooi het pas na parsen in postgis. Uit interesse ben ik eens gaan kijken wat ST_CurveToLine(b.geometrie_vlak) nu produceert in postgis, om erachter te komen wat nu de standaardinstelling van postgis is.

Als ik in de documentatie kijk, begrijp ik dat er standaard 32 segmenten per kwart worden gemaakt. De output van postgis lijkt nogal onregelmatig. I het gegeven voorbeeld (lokaalid = 'G0307.f86336b663854a3b9d00dd4d89a49b56') krijg je de ene keer voor een curve met een lengte van 3 meter 4 punten. De andere keer 30. Ik kan er niet echt chocolade van maken. Als ik kijk in pdok, dan geeft diezelfde lokaalid voor de gecurvde delen zo te zien steeds 100 segmenten.

Zoals ik het nu zie zou je wel aan de knoppen van de ST_CurveToLine moeten draaien om wat helders te krijgen als (multi)polygon, maar ik kan het natuurlijk weer helemaal mis hebben. Het gaat uiteindelijk om de tolerantie. Dus

Om te weten of we het over hetzelfde hebben zou het goed zijn om ook je output van de ST_Multi(ST_Buffer(ST_CurveToLine(b.geometrie_vlak), 0))::geometry(MultiPolygon, 28992) as geometrie_vlak query er eens in te gooien. Steeds voor dat populaire fietspad in Amersfoort natuurlijk.

In rust heb ik die ellende niet. Ik pak gewoon de lengte van de arc en neem een aantal segmenten op basis van een tolerantie die past bij de bgt. Elke arc wordt dan omgezet in punten die ongeveer op dezelfde afstand van elkaar liggen. ST_CurveToLine(b.geometrie_vlak),0.001,1,1) zorgt er voor dat onze polygon maximaal 0.001 van de cirkel afligt, waarbij de laatste 1 voor een gelijke verdeling van segmenten zou moeten zorgen, maar die heb ik nog niet nagekeken.

[willemhoffmans]

ST_Multi(ST_Buffer(ST_CurveToLine(b.geometrie_vlak), 0))::geometry(MultiPolygon, 28992) as geometrie_vlak De ST_CurveToLine is hierbij m.i. overbodig, evenals de casting naar ::geometry(MultiPolygon, 28992)

ST_Buffer maakt vanzelf al een geometrie zonder curves. Overigens constateer ik dat ST_Buffer(geometrie_vlak, 0)32 segmenten per kwart maakt, terwijl de in de documentatie aangegeven default 8 segmenten is. Ik heb dat uitgeprobeerd op een cirkel. Bufferen met verschillende toleranties levert dan steeds 32 segmenten op, het maakt dus niks uit welke je kiest. Bijvoorbeeld: ST_Buffer(geometrie_vlak, 0, 1) of ST_Buffer(geometrie_Vlak, 0, 100)

Steeds 32 segmenten!

[justb4]

@janwinsemius denk inderdaad dat de parameters voor ST_CurveToLine() 'op maat' moeten zijn. Ik had ook verschillende waarden uitgeprobeerd, maar op complete BGT dook er altijd wel weer een fout op. Per geometrie kunnen er meerdere 'bogen', bijv CIRCULARSTRINGs voorkomen, met ieder eigen afstanden. Weet niet of PostGIS dergelijke functies bezit om dat te analyseren.

@willemhoffmans ja is dubbel-dubbelop. Dan zal denk ik ST_Buffer intern een ST_CurveToLine doen (?), maar voor begrip wel helder die erin te houden. De cast naar MultiPolygon is nog erfenis, voordat ST_Multi gebruikt werd, omdat ST_CurveToLine in vele/meeste gevallen Polygon opleverde

Maar belangrijkste, althans voor mij, is dat er werkbare geometrie is, en er niets verloren gaat. De aansluiting zal niet perfect zijn, maar verwacht voor 'mapping' toepassing bijv voldoende.

[justb4]

Ik sluit deze, er is een werkbare situatie waarin m.i. geen geometrieën verloren gaan. Als er toch problemen of verbeteringen nodig zijn, gaarne nieuw issue.