Générer une version de données où les EPCI des DROM sont regroupés à proximité du territoire métropolitain

Cette demande fait suite à une demande pour avoir un résultat similaire à "Contours des communes de France simplifié, avec régions et département d'outre-mer rapprochés" https://www.data.gouv.fr/fr/datasets/contours-des-communes-de-france-simplifie-avec-regions-et-departement-doutre-mer-rapproches/

Nous avions vu passer le billet de blog d'Eric Mauvière à ce sujet https://www.icem7.fr/cartographie/un-fond-de-carte-france-par-commune-optimise-pour-le-web-et-lanalyse-statistique/

Méthodologie adoptée

Il passait par les codes départements pour déterminer quelles entités devaient être translatés/redimensionnées sur la métropole. Nous avons fait le choix de partir des rectangles englobants ("bounding box") pour pouvoir assigner les mêmes codes et réutiliser une partie de la recette du billet de blog. Ce fonctionnement présente en outre l'avantage d'être générique à toutes les unités administratives indépendamment que leurs colonnes originales comporte une référence au code département/drom/com.

Nos choix pour réassigner des codes sont les suivants:

Le code '000' regroupe la zone métropolitaine en incluant la Corse. Elle n'a rien d'officielle mais nous permet de faire notre tri dans les objets à translater/redimensionner

Les codes directement attribués sont:

Mayotte (976)
La Réunion (974)
Guyane (973)
Martinique (972)
Guadeloupe (971)
Saint-Pierre-et-Miquelon (975)
Saint-Barthélemy (977)
Saint-Martin (978)
Terres australes et antarctiques françaises (984)
Wallis et Futuna (986)
Polynésie française (987)
Nouvelle-Calédonie (988)
Île de Clipperton (989)

On exclut les codes 975, 977, 978, 984, 986, 987, 988, 989 car on sait qu'ils n'ont pas d'EPCIs et surtout ultérieurement on n'a pas prévue de translation pour ces unités administratives.

Traitement des données en pratique

On récupère la donnée

wget http://etalab-datasets.geo.data.gouv.fr/contours-administratifs/2024/geojson/epci-1000m.geojson

On debugge nos if/else lors de nos essais initiaux. Attention, la syntaxe utilisée dans la requête SQL CASE WHEN condition1 THEN result1 WHEN condition2 THEN result2 ELSE alternative END ne fonctionne que si vous avez une version de GDAL récente car associée à une version récente de SQLite.

ogr2ogr -f CSV -dialect SQLite -sql "SELECT * FROM (SELECT *, CASE WHEN ST_WITHIN(geometry, ST_Envelope(GeomFromText('MULTIPOINT(44.8 -13.1,45.5 -12.6)', 4326))) THEN '976' WHEN ST_WITHIN(geometry, ST_Envelope(GeomFromText('MULTIPOINT(55.1 -21.5,56.1 -20.8)', 4326))) THEN '974' WHEN ST_WITHIN(geometry, ST_Envelope(GeomFromText('MULTIPOINT(-57.1 1.3,-49.9 6.2)', 4326))) THEN '973' WHEN ST_WITHIN(geometry, ST_Envelope(GeomFromText('MULTIPOINT(-61.5 14.3,-60.6 14.9)', 4326))) THEN '972' WHEN ST_WITHIN(geometry, ST_Envelope(GeomFromText('MULTIPOINT(-62.2 15.6,-60.6 16.7)', 4326))) THEN '971' WHEN ST_WITHIN(geometry, ST_Envelope(GeomFromText('MULTIPOINT(-5.7 40.8,11.1 52.2)', 4326))) THEN '000' WHEN ST_WITHIN(geometry, ST_Envelope(GeomFromText('MULTIPOINT(-56.4556199364351 46.7007782234422,-56.0764709868039 46.7007782234422)', 4326))) THEN '975' WHEN ST_WITHIN(geometry, ST_Envelope(GeomFromText('MULTIPOINT(-62.960429718117 17.8218601145963,-62.7390391961512 17.8218601145963)', 4326))) THEN '977' WHEN ST_WITHIN(geometry, ST_Envelope(GeomFromText('MULTIPOINT(-63.2033251851506 17.9965904153951,-62.921750022029 17.9965904153951)', 4326))) THEN '978' WHEN ST_WITHIN(geometry, ST_Envelope(GeomFromText('MULTIPOINT(-178.23161300831 -14.4121737535224,-176.112320001801 -14.4121737535224)', 4326))) THEN '986' WHEN ST_WITHIN(geometry, ST_Envelope(GeomFromText('MULTIPOINT(-154.772612565755 -27.7024066941686,-134.402483042045 -27.7024066941686)', 4326))) THEN '987' WHEN ST_WITHIN(geometry, ST_Envelope(GeomFromText('MULTIPOINT(163.519710180826 -22.8976911109928,168.183654941539 -22.8976911109928)', 4326))) THEN '988' WHEN ST_WITHIN(geometry, ST_Envelope(GeomFromText('MULTIPOINT(-109.284589731653 10.237159563632,-109.149831814828 10.237159563632)', 4326))) THEN '989' WHEN ST_WITHIN(geometry, ST_Envelope(GeomFromText('MULTIPOINT(39.5914244551481 -85.1,142.05 -85.1)', 4326))) THEN '984' ELSE '999' END AS dpt_simpl FROM \"epci-1000m\") AS foo WHERE dpt_simpl NOT IN ('975', '977', '978', '984', '986', '987', '988', '989', '999')" /vsistdout/ epci-1000m.geojson

On produit le GeoJSON avec les codes dept/drom/com

ogr2ogr -f GeoJSON -dialect SQLite -sql "SELECT * FROM (SELECT *, CASE WHEN ST_WITHIN(geometry, ST_Envelope(GeomFromText('MULTIPOINT(44.8 -13.1,45.5 -12.6)', 4326))) THEN '976' WHEN ST_WITHIN(geometry, ST_Envelope(GeomFromText('MULTIPOINT(55.1 -21.5,56.1 -20.8)', 4326))) THEN '974' WHEN ST_WITHIN(geometry, ST_Envelope(GeomFromText('MULTIPOINT(-57.1 1.3,-49.9 6.2)', 4326))) THEN '973' WHEN ST_WITHIN(geometry, ST_Envelope(GeomFromText('MULTIPOINT(-61.5 14.3,-60.6 14.9)', 4326))) THEN '972' WHEN ST_WITHIN(geometry, ST_Envelope(GeomFromText('MULTIPOINT(-62.2 15.6,-60.6 16.7)', 4326))) THEN '971' WHEN ST_WITHIN(geometry, ST_Envelope(GeomFromText('MULTIPOINT(-5.7 40.8,11.1 52.2)', 4326))) THEN '000' WHEN ST_WITHIN(geometry, ST_Envelope(GeomFromText('MULTIPOINT(-56.4556199364351 46.7007782234422,-56.0764709868039 46.7007782234422)', 4326))) THEN '975' WHEN ST_WITHIN(geometry, ST_Envelope(GeomFromText('MULTIPOINT(-62.960429718117 17.8218601145963,-62.7390391961512 17.8218601145963)', 4326))) THEN '977' WHEN ST_WITHIN(geometry, ST_Envelope(GeomFromText('MULTIPOINT(-63.2033251851506 17.9965904153951,-62.921750022029 17.9965904153951)', 4326))) THEN '978' WHEN ST_WITHIN(geometry, ST_Envelope(GeomFromText('MULTIPOINT(-178.23161300831 -14.4121737535224,-176.112320001801 -14.4121737535224)', 4326))) THEN '986' WHEN ST_WITHIN(geometry, ST_Envelope(GeomFromText('MULTIPOINT(-154.772612565755 -27.7024066941686,-134.402483042045 -27.7024066941686)', 4326))) THEN '987' WHEN ST_WITHIN(geometry, ST_Envelope(GeomFromText('MULTIPOINT(163.519710180826 -22.8976911109928,168.183654941539 -22.8976911109928)', 4326))) THEN '988' WHEN ST_WITHIN(geometry, ST_Envelope(GeomFromText('MULTIPOINT(-109.284589731653 10.237159563632,-109.149831814828 10.237159563632)', 4326))) THEN '989' WHEN ST_WITHIN(geometry, ST_Envelope(GeomFromText('MULTIPOINT(39.5914244551481 -85.1,142.05 -85.1)', 4326))) THEN '984' ELSE '999' END AS dpt_simpl FROM \"epci-1000m\") AS foo WHERE dpt_simpl NOT IN ('975', '977', '978', '984', '986', '987', '988', '989')" -lco WRITE_NAME=NO epci-1000m-dept-simpl.geojson epci-1000m.geojson

Translater et redimensionner

En reprenant la recette de Éric Mauvière sur https://www.icem7.fr/cartographie/un-fond-de-carte-france-par-commune-optimise-pour-le-web-et-lanalyse-statistique/, je peux maintenant appliquer la redimension et translation des coordonnées. J'ai exclu les COM (Collectivités d'Outre Mer) car je ne voulais pas dans l'immédiat me préoccuper de trouver la "bonne" translation et redimension pour ces entités. Les EPCIs ne sont pas concernés par ce cas.

Version horizontale

mapshaper -i epci-1000m-dept-simpl.geojson -proj webmercator \
          -affine where="dpt_simpl.indexOf('971')==0" shift=6355000,3330000 scale=1.5 \
          -affine where="dpt_simpl.indexOf('972')==0" shift=6480000,3505000 scale=1.5 \
          -affine where="dpt_simpl.indexOf('973')==0" shift=5760000,4720000 scale=0.35 \
          -affine where="dpt_simpl.indexOf('974')==0" shift=-6170000,7560000 scale=1.5 \
          -affine where="dpt_simpl.indexOf('976')==0" shift=-4885000,6590000 scale=1.5 \
          -filter-fields code,nom \
          -proj wgs84 \
          -o precision=0.0001 epci-1000m-dept-simpl-displaced-horizontal.geojson

Version verticale

mapshaper -i epci-1000m-dept-simpl.geojson -proj webmercator \
          -affine where="dpt_simpl.indexOf('971')==0" shift=6455000,4030000 scale=1.5 \
          -affine where="dpt_simpl.indexOf('972')==0" shift=6400000,4075000 scale=1.5 \
          -affine where="dpt_simpl.indexOf('973')==0" shift=5540000,5140000 scale=0.35 \
          -affine where="dpt_simpl.indexOf('974')==0" shift=-6540000,7830000 scale=1.5 \
          -affine where="dpt_simpl.indexOf('976')==0" shift=-5375000,6740000 scale=1.5 \
          -filter-fields code,nom \
          -proj wgs84 \
          -o precision=0.0001 epci-1000m-dept-simpl-displaced-vertical.geojson

Rendus

Aperçu du rendu horizontal

mapshaper -i epci-1000m-dept-simpl.geojson -proj webmercator \
          -affine where="dpt_simpl.indexOf('971')==0" shift=6355000,3330000 scale=1.5 \
          -affine where="dpt_simpl.indexOf('972')==0" shift=6480000,3505000 scale=1.5 \
          -affine where="dpt_simpl.indexOf('973')==0" shift=5760000,4720000 scale=0.35 \
          -affine where="dpt_simpl.indexOf('974')==0" shift=-6170000,7560000 scale=1.5 \
          -affine where="dpt_simpl.indexOf('976')==0" shift=-4885000,6590000 scale=1.5 \
          -proj EPSG:2154 \
          -each 'idwpref="e" + code' \
          -each 'idreg="r" + dpt_simpl' \
          -colorizer name=calcFill colors='rgba(255, 0, 0, 0.1),rgba(0, 0, 255, 0.1),rgba(0, 128, 0, 0.1),rgba(255, 0, 255, 0.1),rgba(255, 255, 0, 0.1),rgba(0, 255, 255, 0.1)' categories='000,971,972,973,974,976' \
          -style fill='calcFill(dpt_simpl)' stroke='#aaa' stroke-width=0.6 \
          -o svg-data=idwpref,idreg \
          epci-1000m-dept-simpl-displaced-2154-horizontal.svg

Le résultat est le suivant:

Aperçu du rendu vertical

mapshaper -i epci-1000m-dept-simpl.geojson -proj webmercator \
          -affine where="dpt_simpl.indexOf('971')==0" shift=6455000,4030000 scale=1.5 \
          -affine where="dpt_simpl.indexOf('972')==0" shift=6400000,4075000 scale=1.5 \
          -affine where="dpt_simpl.indexOf('973')==0" shift=5540000,5140000 scale=0.35 \
          -affine where="dpt_simpl.indexOf('974')==0" shift=-6540000,7830000 scale=1.5 \
          -affine where="dpt_simpl.indexOf('976')==0" shift=-5375000,6740000 scale=1.5 \
          -proj EPSG:2154 \
          -each 'idwpref="e" + code' \
          -each 'idreg="r" + dpt_simpl' \
          -colorizer name=calcFill colors='rgba(255, 0, 0, 0.1),rgba(0, 0, 255, 0.1),rgba(0, 128, 0, 0.1),rgba(255, 0, 255, 0.1),rgba(255, 255, 0, 0.1),rgba(0, 255, 255, 0.1)' categories='000,971,972,973,974,976' \
          -style fill='calcFill(dpt_simpl)' stroke='#aaa' stroke-width=0.6 \
          -o svg-data=idwpref,idreg \
          epci-1000m-dept-simpl-displaced-2154-vertical.svg