-- ****************************************************************************
-- vérification de la validité des géométries (+ "correction")
-- ****************************************************************************

SELECT ST_IsValidReason(geom), *
FROM commune
WHERE NOT ST_IsValid(geom);

SELECT ST_NumGeometries(geom)
FROM commune
WHERE cog = '07091';

UPDATE commune
SET geom = ST_MakeValid(geom)
WHERE NOT ST_IsValid(geom);

-- ****************************************************************************
-- recherche de superposition entre géométries
-- ****************************************************************************

DROP TABLE IF EXISTS erreur_superposition;

CREATE TABLE erreur_superposition AS
SELECT
    c1.cog AS cog_c1, c1.nom AS nom_c1,
    c2.cog AS cog_c2, c2.nom AS nom_c2,
    (ST_Dump(ST_CollectionExtract(ST_Intersection(c1.geom, c2.geom), 3))).geom AS geom
FROM commune c1
INNER JOIN commune c2
ON c1.cog > c2.cog
AND ST_Intersects(c1.geom, c2.geom)
AND NOT ST_Touches(c1.geom, c2.geom);

-- ****************************************************************************
-- recherche de zones blanches entre les géométries
-- ****************************************************************************

DROP TABLE IF EXISTS erreur_zone_blanche;

CREATE TABLE erreur_zone_blanche AS
WITH union_commune AS (
    SELECT ST_Union(geom) AS geom
    FROM commune
    WHERE cog_departement = '30'
)
SELECT gs AS numero_zone, ST_MakePolygon(ST_InteriorRingN(geom, gs)) AS geom
FROM union_commune uc
CROSS JOIN generate_series(1, ST_NumInteriorRings(uc.geom)) gs;

-- ****************************************************************************
-- plus proche voisin textuel
-- ****************************************************************************

CREATE EXTENSION IF NOT EXISTS unaccent;
SELECT unaccent('pgQuaDoGéo - Qualité des données ? Ne vous y trompez pas, adoptez 2 éléphants !');

CREATE EXTENSION IF NOT EXISTS pg_trgm;
SELECT similarity('Bienvenue chez les ch''tis', 'Bienvenue chez les rozes');

SELECT * FROM etablissement_sirene;

DROP TABLE IF EXISTS geocodate_etablissement_sirene;

-- voir :
--  802555,6319899
CREATE TABLE geocodate_etablissement_sirene AS
WITH details_adresse AS (
    SELECT
        a.numero,
        CASE WHEN a.repetition = '' THEN NULL ELSE a.repetition END AS repetition,
        a.voie,
        c.cog AS cog_commune,
        c.nom AS nom_commune,
        a.geom
    FROM adresse a
    INNER JOIN commune c ON a.cog_commune = c.cog
)
SELECT
    siren,
    nic,
    siret,
    concat_ws(' ', e.numerovoieetablissement, e.indicerepetitionetablissement, e.typevoieetablissement, e.libellevoieetablissement) AS adresse_saisie,
    concat_ws(' ', ppv.numero, ppv.repetition, ppv.voie) AS adresse_trouvee,
    ppv.cog_commune,
    ppv.nom_commune,
    ppv.taux_similarite,
    ppv.geom
FROM etablissement_sirene e
LEFT JOIN LATERAL (
    SELECT
        da.*,
        similarity(
            concat_ws(' ', da.numero, da.repetition, da.voie),
            concat_ws(' ', e.numerovoieetablissement, e.indicerepetitionetablissement, e.typevoieetablissement, e.libellevoieetablissement)
        ) AS taux_similarite
    FROM details_adresse da
    WHERE upper(unaccent(da.nom_commune)) = upper(unaccent(e.libellecommuneetablissement))
    ORDER BY taux_similarite DESC
    LIMIT 1
) ppv ON ppv.taux_similarite > 0.3;

SELECT *
FROM geocodate_etablissement_sirene
WHERE geom IS NULL;

-- ****************************************************************************
-- plus proche voisin géographique
-- ****************************************************************************

DROP INDEX IF EXISTS adresse_geom_idx;

CREATE INDEX adresse_geom_idx ON adresse USING gist (geom);

DROP TABLE IF EXISTS ppv_batiment_adresse;

-- voir :
--  810570.8,6292533.2
CREATE TABLE ppv_batiment_adresse AS
SELECT eloignement, b.*, ST_ShortestLine(b.Geom, ppv.Geom) AS geomlien
FROM batiment b
CROSS JOIN LATERAL (
    SELECT *, a.geom <-> b.geom AS eloignement
    FROM adresse a
    ORDER BY eloignement ASC
    LIMIT 1
) ppv
ORDER BY eloignement DESC;

-- ****************************************************************************
-- calcul de l'étendue des voies
-- ****************************************************************************

DROP TABLE IF EXISTS etendue_voie;

-- voir :
--  815507.1,6287858.6
CREATE TABLE etendue_voie AS
SELECT
    cog_commune,
    voie,
    ST_MaxDistance(ST_Collect(geom), ST_Collect(geom)) AS etendue,
    ST_Collect(geom) AS geom
FROM adresse
GROUP BY cog_commune, voie
ORDER BY etendue DESC;

-- ****************************************************************************
-- calcul de l'éloignement des adresses successives d'une même rue
-- ****************************************************************************

DROP TABLE IF EXISTS eloignement_adresses_successives;

CREATE TABLE eloignement_adresses_successives AS
WITH adresse_et_suivante AS (
    SELECT
        *,
        LEAD (numero) OVER (PARTITION BY voie, cog_commune ORDER BY numero, repetition) AS numero_suivant,
        LEAD (repetition) OVER (PARTITION BY voie, cog_commune ORDER BY numero, repetition) AS repetition_suivante,
        LEAD (geom) OVER (PARTITION BY voie, cog_commune ORDER BY numero, repetition) AS geom_suivante
    FROM adresse
)
SELECT
    cog_commune,
    voie,
    concat_ws(' ', numero, repetition, '->', numero_suivant, repetition_suivante) AS succession,
    geom <-> geom_suivante AS eloignement,
    ST_MakeLine(geom, geom_suivante) AS geomlien
FROM adresse_et_suivante
WHERE geom_suivante IS NOT NULL;

-- ****************************************************************************
-- recherche stricte des évolutions entre 2 versions de la BDTopo
-- (basée sur les attributs)
-- ****************************************************************************

-- voir :
--  786123.8,6298145.3
--  789073.82,6340483.62
--  788144.9,6306479.0
--  793156.0,6315399.0
DROP TABLE IF EXISTS evolution_troncon_stricte;

CREATE TABLE evolution_troncon_stricte AS
WITH departementale2021 AS (
    SELECT *
    FROM bdtopo_troncon_de_route_20211215
    WHERE cpx_classement_administratif = 'Départementale'
    AND cpx_gestionnaire = 'Gard'
),
departementale2022 AS (
    SELECT *
    FROM bdtopo_troncon_de_route_20221215
    WHERE cpx_classement_administratif = 'Départementale'
    AND cpx_gestionnaire = 'Gard'
)
SELECT
    CASE
        WHEN v2021.cleabs IS NULL AND v2022.cleabs IS NOT NULL THEN 'ajout'
        WHEN v2022.cleabs IS NULL AND v2021.cleabs IS NOT NULL THEN 'retrait'
        ELSE 'modification'
    END AS action,
    CASE
        WHEN v2022.geom IS NULL THEN v2021.geom
        ELSE v2022.geom
    END AS geom
FROM departementale2021 v2021
FULL JOIN departementale2022 v2022 ON v2021.cleabs = v2022.cleabs
WHERE v2021.cleabs IS NULL
OR v2022.cleabs IS NULL
OR NOT ST_Equals(v2021.geom, v2022.geom);

-- ****************************************************************************
-- recherche lâche des évolutions entre 2 versions de la BDTopo
-- (basée sur la géométrie)
-- ****************************************************************************

-- voir :
--  786123.8,6298145.3
--  789073.82,6340483.62
--  788144.9,6306479.0
--  793156.0,6315399.0
DROP TABLE IF EXISTS evolution_troncon_lache;

CREATE TABLE evolution_troncon_lache AS
WITH departementale2021 AS (
    SELECT cleabs, cpx_numero, geom
    FROM bdtopo_troncon_de_route_20211215
    WHERE cpx_classement_administratif = 'Départementale'
    AND cpx_gestionnaire = 'Gard'
),
departementale2022 AS (
    SELECT cleabs, cpx_numero, geom
    FROM bdtopo_troncon_de_route_20221215
    WHERE cpx_classement_administratif = 'Départementale'
    AND cpx_gestionnaire = 'Gard'
),
buffer2021 AS (
    SELECT cpx_numero, ST_Union(ST_Buffer(geom, 5)) AS geom
    FROM departementale2021
    GROUP BY cpx_numero
),
buffer2022 AS (
    SELECT cpx_numero, ST_Union(ST_Buffer(geom, 5)) AS geom
    FROM departementale2022
    GROUP BY cpx_numero
)
SELECT 'retrait' AS action, d21.*
FROM departementale2021 d21
LEFT JOIN buffer2022 b22 ON b22.cpx_numero = d21.cpx_numero AND ST_Covers(b22.geom, d21.geom)
WHERE b22.geom IS NULL
UNION
SELECT 'ajout' AS action, d22.*
FROM departementale2022 d22
LEFT JOIN buffer2021 b21 ON b21.cpx_numero = d22.cpx_numero AND ST_Covers(b21.geom, d22.geom)
WHERE b21.geom IS NULL;

-- ****************************************************************************
-- IGN vs OSM
-- ****************************************************************************

/*
Export OSM réalisé via https://overpass-turbo.eu/
area
  [admin_level=6]
  [name="Gard"];
way
  [highway]
  (area);
out geom;
*/

-- comparaison IGN vs OSM "sévère" (c-a-d tronçon complet)

-- voir :
--  763135.0,6323781.9
DROP TABLE IF EXISTS comparaison_ign_osm_severe;

CREATE TABLE comparaison_ign_osm_severe AS
WITH departementaleign AS (
    SELECT cleabs AS idsource, cpx_numero AS numero, geom
    FROM bdtopo_troncon_de_route_20221215
    WHERE cpx_classement_administratif = 'Départementale'
    AND cpx_gestionnaire = 'Gard'
    AND nature <> 'Rond-point'
),
departementaleosm AS (
    SELECT id AS idsource, upper(regexp_replace(ref, '^D ', 'D')) AS numero, geom
    FROM osm_highway_20230226
    WHERE ref ~ '^D \d+'
),
bufferign AS (
    SELECT numero, ST_Union(ST_Buffer(geom, 5)) AS geom
    FROM departementaleign
    GROUP BY numero
),
bufferosm AS (
    SELECT numero, ST_Union(ST_Buffer(geom, 5)) AS geom
    FROM departementaleosm
    GROUP BY numero
)
SELECT 'ommission' AS statut, dign.*
FROM departementaleign dign
LEFT JOIN bufferosm bosm ON bosm.numero = dign.numero AND ST_Covers(bosm.geom, dign.geom)
WHERE bosm.geom IS NULL
UNION
SELECT 'excédent' AS statut, dosm.*
FROM departementaleosm dosm
LEFT JOIN bufferign bign ON bign.numero = dosm.numero AND ST_Covers(bign.geom, dosm.geom)
WHERE bign.geom IS NULL;

-- calcul du taux d'exhaustivité sur la base de la comparaison "sévère"

WITH longueur_referentiel AS (
    SELECT sum(ST_Length(geom)) AS longueur
    FROM bdtopo_troncon_de_route_20221215
    WHERE cpx_classement_administratif = 'Départementale'
    AND cpx_gestionnaire = 'Gard'
    AND nature <> 'Rond-point'
),
longueur_excedent AS (
    SELECT sum(ST_Length(geom)) AS longueur
    FROM comparaison_ign_osm_severe
    WHERE statut = 'excédent'
),
longueur_ommission AS (
    SELECT sum(ST_Length(geom)) AS longueur
    FROM comparaison_ign_osm_severe
    WHERE statut = 'ommission'
)
SELECT
    round(r.longueur::numeric / 1000, 2) AS "Longueur référentiel (km)",
    round(e.longueur::numeric / 1000, 2) AS "Longueur excédents (km)",
    round(o.longueur::numeric / 1000, 2) AS "Longueur ommissions (km)",
    round(((1 - (e.longueur + o.longueur) / r.longueur) * 100)::numeric, 2) "Taux exhaustivité (%)"
FROM longueur_referentiel r, longueur_excedent e, longueur_ommission o;

-- comparaison IGN vs OSM "juste" (c-a-d part de tronçon)

DROP TABLE IF EXISTS comparaison_ign_osm_juste;

CREATE TABLE comparaison_ign_osm_juste AS
WITH departementaleign AS (
    SELECT cleabs AS idsource, cpx_numero AS numero, geom
    FROM bdtopo_troncon_de_route_20221215
    WHERE cpx_classement_administratif = 'Départementale'
    AND cpx_gestionnaire = 'Gard'
    AND nature <> 'Rond-point'
),
departementaleosm AS (
    SELECT id AS idsource, upper(regexp_replace(ref, '^D ', 'D')) AS numero, geom
    FROM osm_highway_20230226
    WHERE ref ~ '^D \d+'
),
bufferign AS (
    SELECT numero, ST_Union(ST_Buffer(geom, 5)) AS geom
    FROM departementaleign
    GROUP BY numero
),
bufferosm AS (
    SELECT numero, ST_Union(ST_Buffer(geom, 5)) AS geom
    FROM departementaleosm
    GROUP BY numero
)
SELECT
    'ommission' AS statut,
    dign.idsource,
    dign.numero,
    CASE
        WHEN bosm.geom IS NULL THEN dign.geom
        ELSE ST_Difference(dign.geom, bosm.geom)
    END AS geom
FROM departementaleign dign
LEFT JOIN bufferosm bosm ON bosm.numero = dign.numero
WHERE bosm.geom IS NULL OR NOT ST_Covers(bosm.geom, dign.geom)
UNION
SELECT
    'excédent' AS statut,
    dosm.idsource,
    dosm.numero,
    CASE
        WHEN bign.geom IS NULL THEN dosm.geom
        ELSE ST_Difference(dosm.geom, bign.geom)
    END AS geom
FROM departementaleosm dosm
LEFT JOIN bufferign bign ON bign.numero = dosm.numero
WHERE bign.geom IS NULL OR NOT ST_Covers(bign.geom, dosm.geom);

-- calcul du taux d'exhaustivité sur la base de la comparaison "juste"

WITH longueur_referentiel AS (
    SELECT sum(ST_Length(geom)) AS longueur
    FROM bdtopo_troncon_de_route_20221215
    WHERE cpx_classement_administratif = 'Départementale'
    AND cpx_gestionnaire = 'Gard'
    AND nature <> 'Rond-point'
),
longueur_excedent AS (
    SELECT sum(ST_Length(geom)) AS longueur
    FROM comparaison_ign_osm_juste
    WHERE statut = 'excédent'
),
longueur_ommission AS (
    SELECT sum(ST_Length(geom)) AS longueur
    FROM comparaison_ign_osm_juste
    WHERE statut = 'ommission'
)
SELECT
    round(r.longueur::numeric / 1000, 2) AS "Longueur référentiel (km)",
    round(e.longueur::numeric / 1000, 2) AS "Longueur excédents (km)",
    round(o.longueur::numeric / 1000, 2) AS "Longueur ommissions (km)",
    round(((1 - (e.longueur + o.longueur) / r.longueur) * 100)::numeric, 2) "Taux exhaustivité (%)"
FROM longueur_referentiel r, longueur_excedent e, longueur_ommission o;

-- ****************************************************************************
-- recherche des giratoires OSM non fermés
-- ****************************************************************************

DROP TABLE IF EXISTS giratoire_osm_non_ferme;

CREATE TABLE giratoire_osm_non_ferme AS
WITH giratoire AS (
    SELECT ST_LineMerge(ST_CollectionExtract(unnest(ST_ClusterIntersecting(geom)), 2)) as geom
    FROM osm_highway_20230226
    WHERE junction = 'roundabout'
)
SELECT row_number() over() as id, geom
FROM giratoire
WHERE NOT ST_IsClosed(geom);

-- ****************************************************************************
-- exemple de requêtes récursives
-- ****************************************************************************

DROP TABLE IF EXISTS orienter_troncon;

CREATE TABLE orienter_troncon AS
WITH RECURSIVE orienter AS (
    SELECT
        cpx_numero,
        cleabs::varchar AS chemin,
        cleabs,
        liens_vers_route_nommee,
        0::numeric as distancedebut,
        ST_Length(geom)::numeric as distancefin,
        geom::geometry
    FROM bdtopo_troncon_de_route_20221215
    WHERE cleabs = 'TRONROUT0000000027900538' -- premier tronçon de la D1
    
    UNION
    
    SELECT
        ppv.cpx_numero,
        concat_ws('|', o.chemin, ppv.cleabs) AS chemin,
        ppv.cleabs,
        ppv.liens_vers_route_nommee,
        o.distancefin as distancedebut,
        o.distancefin + ST_Length(ppv.geom)::numeric as distancefin,
        ppv.geom
    FROM orienter o
    CROSS JOIN LATERAL (
        SELECT
            cpx_numero,
            cleabs,
            liens_vers_route_nommee,
            CASE
                WHEN t.sens_de_circulation = 'Sens inverse' THEN ST_Reverse(t.geom)
                WHEN t.sens_de_circulation = 'Sens direct' THEN t.geom
                WHEN ST_EndPoint(o.geom) <-> ST_EndPoint(t.geom) < ST_EndPoint(o.geom) <-> ST_StartPoint(t.geom) THEN ST_Reverse(t.geom)
                ELSE t.geom
            END AS geom
        FROM bdtopo_troncon_de_route_20221215 t
        WHERE o.liens_vers_route_nommee = t.liens_vers_route_nommee -- le tronçon doit être de la même RD
        AND o.chemin not like '%' || t.cleabs || '%' -- et le tronçon ne doit pas déjà participer au chemin (pour éviter la récursivité infinie)
        ORDER BY ST_EndPoint(o.geom) <-> t.geom
        LIMIT 1
    ) ppv
)
SELECT cpx_numero, cleabs, liens_vers_route_nommee, distancedebut, distancefin, geom, ST_AddMeasure(geom, distancedebut, distancefin) as geom_avec_mesures
FROM orienter;

DROP TABLE IF EXISTS point24500;

CREATE TABLE point24500 AS
SELECT ST_Force2D(ST_LocateAlong(geom_avec_mesures, 24500)) as geom
FROM orienter_troncon
WHERE 24500 BETWEEN ST_M(ST_StartPoint(geom_avec_mesures)) AND ST_M(ST_EndPoint(geom_avec_mesures));