Correct polygon winding directions

camdecoster · camdecoster · commit 15dd77a6724d · 2025-06-25T14:21:19.000-06:00
diff --git a/tasks/topojson/process_geodata.mjs b/tasks/topojson/process_geodata.mjs
@@ -1,12 +1,10 @@
+import rewind from '@mapbox/geojson-rewind'
 import { geoIdentity, geoPath } from 'd3-geo';
-import { geoStitch } from 'd3-geo-projection'
 import fs from 'fs';
 import mapshaper from 'mapshaper';
 import path from 'path';
-import config, { getNEFilename } from './config.mjs';
-import { topology } from 'topojson-server';
 import topojsonLib from 'topojson';
-import rewind from '@mapbox/geojson-rewind'
+import config, { getNEFilename } from './config.mjs';
 
 const { filters, inputDir, layers, resolutions, scopes, unFilename, vectors } = config;
 
@@ -32,19 +30,6 @@ function getJsonFile(filename) {
     }
 }
 
-async function createCountriesLayer({ bounds, filter, name, resolution, source }) {
-    const inputFilePath = `${outputDirGeojson}/${unFilename}_${resolution}m/${source}.geojson`;
-    const outputFilePath = `${outputDirGeojson}/${name}_${resolution}m/countries.geojson`;
-    const commands = [
-        inputFilePath,
-        bounds.length ? `-clip bbox=${bounds.join(',')}` : '',
-        filter ? `-filter '${filter}'` : '',
-        `-o ${outputFilePath}`
-    ].join(' ');
-    await mapshaper.runCommands(commands);
-    addCentroidsToGeojson(outputFilePath);
-}
-
 function addCentroidsToGeojson(geojsonPath) {
     const geojson = getJsonFile(geojsonPath);
     if (!geojson.features) return;
@@ -59,6 +44,100 @@ function addCentroidsToGeojson(geojsonPath) {
     fs.writeFileSync(geojsonPath, JSON.stringify({ ...geojson, features }));
 }
 
+// Wind the polygon rings in the correct direction to indicate what is solid and what is whole
+const rewindGeojson = (geojson, clockwise = true) => rewind(geojson, clockwise)
+
+// Snap x-coordinates that are close to be on the antimeridian
+function snapToAntimeridian(inputFilepath, outputFilepath) {
+    outputFilepath ||= inputFilepath
+    const jsonString = fs.readFileSync(inputFilepath, 'utf8')
+    const updatedString = jsonString
+        .replaceAll(/179\.99\d+,/g, '180,')
+        .replaceAll(/180\.00\d+,/g, '180,')
+
+    fs.writeFileSync(outputFilepath, updatedString);
+}
+
+function pruneProperties(topojson) {
+    for (const layer in topojson.objects) {
+        switch (layer) {
+            case 'countries':
+                topojson.objects[layer].geometries = topojson.objects[layer].geometries.map((geometry) => {
+                    const { properties } = geometry;
+                    if (properties) {
+                        geometry.id = properties.iso3cd;
+                        geometry.properties = {
+                            ct: properties.ct
+                        };
+                    }
+
+                    return geometry;
+                });
+                break;
+            case 'subunits':
+                topojson.objects[layer].geometries = topojson.objects[layer].geometries.map((geometry) => {
+                    const { properties } = geometry;
+                    if (properties) {
+                        geometry.id = properties.postal;
+                        geometry.properties = {
+                            ct: properties.ct,
+                            gu: properties.gu_a3
+                        };
+                    }
+
+                    return geometry;
+                });
+
+                break;
+            default:
+                topojson.objects[layer].geometries = topojson.objects[layer].geometries.map((geometry) => {
+                    delete geometry.id;
+                    delete geometry.properties;
+
+                    return geometry;
+                });
+
+                break;
+        }
+    }
+
+    return topojson;
+}
+
+function getCentroid(feature) {
+    const { type } = feature.geometry;
+    const projection = geoIdentity();
+    const path = geoPath(projection);
+
+    if (type === 'MultiPolygon') {
+        let maxArea = -Infinity;
+
+        for (const coordinates of feature.geometry.coordinates) {
+            const polygon = { type: 'Polygon', coordinates };
+            const area = path.area(polygon);
+            if (area > maxArea) {
+                maxArea = area;
+                feature = polygon;
+            }
+        }
+    }
+
+    return path.centroid(feature).map((coord) => +coord.toFixed(2));
+}
+
+async function createCountriesLayer({ bounds, filter, name, resolution, source }) {
+    const inputFilePath = `${outputDirGeojson}/${unFilename}_${resolution}m/${source}.geojson`;
+    const outputFilePath = `${outputDirGeojson}/${name}_${resolution}m/countries.geojson`;
+    const commands = [
+        inputFilePath,
+        bounds.length ? `-clip bbox=${bounds.join(',')}` : '',
+        filter ? `-filter '${filter}'` : '',
+        `-o ${outputFilePath}`
+    ].join(' ');
+    await mapshaper.runCommands(commands);
+    addCentroidsToGeojson(outputFilePath);
+}
+
 async function createLandLayer({ bounds, name, resolution, source }) {
     const inputFilePath = `${outputDirGeojson}/${name}_${resolution}m/countries.geojson`;
     const outputFilePath = `${outputDirGeojson}/${name}_${resolution}m/land.geojson`;
@@ -80,9 +159,12 @@ async function createCoastlinesLayer({ bounds, name, resolution, source }) {
         '-dissolve',
         '-lines',
         bounds.length ? `-clip bbox=${bounds.join(',')}` : '',
+        // Erase outer lines to avoid unpleasant lines through polygons crossing the antimeridian
+        ['antarctica', 'world'].includes(name) ? '-clip bbox=-179.999,-89.999,179.999,89.999' : '',
         `-o ${outputFilePath}`
     ].join(' ');
     await mapshaper.runCommands(commands);
+    if (['antarctica', 'world'].includes(name)) snapToAntimeridian(outputFilePath)
 }
 
 async function createOceanLayer({ bounds, name, resolution, source }) {
@@ -133,96 +215,19 @@ async function createSubunitsLayer({ name, resolution, source }) {
     addCentroidsToGeojson(outputFilePath);
 }
 
-function pruneProperties(topojson) {
-    for (const layer in topojson.objects) {
-        switch (layer) {
-            case 'countries':
-                topojson.objects[layer].geometries = topojson.objects[layer].geometries.map((geometry) => {
-                    const { properties } = geometry;
-                    if (properties) {
-                        geometry.id = properties.iso3cd;
-                        geometry.properties = {
-                            ct: properties.ct
-                        };
-                    }
-
-                    return geometry;
-                });
-                break;
-            case 'subunits':
-                topojson.objects[layer].geometries = topojson.objects[layer].geometries.map((geometry) => {
-                    const { properties } = geometry;
-                    if (properties) {
-                        geometry.id = properties.postal;
-                        geometry.properties = {
-                            ct: properties.ct,
-                            gu: properties.gu_a3
-                        };
-                    }
-
-                    return geometry;
-                });
-
-                break;
-            default:
-                topojson.objects[layer].geometries = topojson.objects[layer].geometries.map((geometry) => {
-                    delete geometry.id;
-                    delete geometry.properties;
-
-                    return geometry;
-                });
-
-                break;
-        }
-    }
-
-    return topojson;
-}
-
-function getCentroid(feature) {
-    const { type } = feature.geometry;
-    const projection = geoIdentity();
-    const path = geoPath(projection);
-
-    if (type === 'MultiPolygon') {
-        let maxArea = -Infinity;
-
-        for (const coordinates of feature.geometry.coordinates) {
-            const polygon = { type: 'Polygon', coordinates };
-            const area = path.area(polygon);
-            if (area > maxArea) {
-                maxArea = area;
-                feature = polygon;
-            }
-        }
-    }
-
-    return path.centroid(feature).map((coord) => +coord.toFixed(2));
-}
-
 async function convertLayersToTopojson({ name, resolution }) {
     const regionDir = path.join(outputDirGeojson, `${name}_${resolution}m`);
     if (!fs.existsSync(regionDir)) return;
 
     const outputFile = `${outputDirTopojson}/${name}_${resolution}m.json`;
+    // Scopes with polygons that cross the antimeridian need to be stitched (via the topology call)
     if (["antarctica", "world"].includes(name)) {
-    // if (false) {
-        const files = fs.readdirSync(regionDir)
         const geojsonObjects = {}
-        for (const file of files) {
-            const filePath = path.join(regionDir, file)
-            const layer = file.split(".")[0]
-            let stitchedGeojson = geoStitch(getJsonFile(filePath))
-            // stitchedGeojson = rewind(stitchedGeojson, true)
-            // fs.writeFileSync(filePath, JSON.stringify(stitchedGeojson));
-            geojsonObjects[layer] = stitchedGeojson
-            // geojsonObjects[layer] = getJsonFile(filePath)
+        for (const layer of Object.keys(config.layers)) {
+            const filePath = path.join(regionDir, `${layer}.geojson`)
+            geojsonObjects[layer] = rewindGeojson(getJsonFile(filePath))
         }
-        const topojsonTopology = topology(geojsonObjects)
-        // const topojsonTopology = topojsonLib.topology(geojsonObjects, {
-        //   verbose: true,
-        //   'property-transform': f => f.properties
-        // })
+        const topojsonTopology = topojsonLib.topology(geojsonObjects, { 'property-transform': f => f.properties })
         fs.writeFileSync(outputFile, JSON.stringify(topojsonTopology));
     } else {
         // Layer names default to file names
@@ -231,28 +236,17 @@ async function convertLayersToTopojson({ name, resolution }) {
     }
 
     // Remove extra information from features
-    // const topojson = getJsonFile(outputFile);
-    // const prunedTopojson = pruneProperties(topojson);
-    // fs.writeFileSync(outputFile, JSON.stringify(prunedTopojson));
+    const topojson = getJsonFile(outputFile);
+    const prunedTopojson = pruneProperties(topojson);
+    fs.writeFileSync(outputFile, JSON.stringify(prunedTopojson));
 }
 
 // Get polygon features from UN GeoJSON and patch Antarctica gap
 const inputFilePathUNGeojson = `${inputDir}/${unFilename}.geojson`;
 const inputFilePathUNGeojsonCleaned = `${inputDir}/${unFilename}_cleaned.geojson`;
-// TODO: Update all x-coords close to 180 to be exactly 180
-function snapToAntimeridian(inputFilepath, outputFilepath) {
-    const jsonString = fs.readFileSync(inputFilepath, 'utf8')
-    const updatedString = jsonString
-        .replaceAll(/179\.99\d+,/g, '180,')
-        .replaceAll(/180\.00\d+,/g, '180,')
-
-    fs.writeFileSync(outputFilepath, updatedString);
-}
 snapToAntimeridian(inputFilePathUNGeojson, inputFilePathUNGeojsonCleaned)
 const commandsAllFeaturesCommon = [
-    // TODO: Should I use the cleaned data or leave as is?
     inputFilePathUNGeojsonCleaned,
-    // inputFilePathUNGeojson,
     `-filter 'iso3cd === "ATA"' target=1 + name=antarctica`,
     // Use 'snap-interval' to patch gap in Antarctica
     '-clean snap-interval=0.015 target=antarctica',
@@ -269,7 +263,8 @@ const commandsAllFeaturesCommon = [
     '-erase source=caspian_sea target=all_features',
     // Update country codes for disputed territories at Egypt/Sudan border: https://en.wikipedia.org/wiki/Egypt%E2%80%93Sudan_border
     `-each 'if (globalid === "{CA12D116-7A19-41D1-9622-17C12CCC720D}") iso3cd = "XHT"'`, // Halaib Triangle
-    `-each 'if (globalid === "{9FD54A50-0BFB-4385-B342-1C3BDEE5ED9B}") iso3cd = "XBT"'` // Bir Tawil
+    `-each 'if (globalid === "{9FD54A50-0BFB-4385-B342-1C3BDEE5ED9B}") iso3cd = "XBT"'`, // Bir Tawil
+    `-each 'FID = iso3cd'`
 ]
 
 // Process 50m UN geodata
@@ -305,13 +300,11 @@ await mapshaper.runCommands(commandsLand50m);
 const inputFilePath110m = outputFilePath50m;
 const outputFilePath110m = `${outputDirGeojson}/${unFilename}_110m/all_features.geojson`;
 const commandsAllFeatures110m = [
-    // ...commandsAllFeaturesCommon,
     inputFilePath110m,
-    '-simplify 10% rdp',
+    '-simplify 20%',
     `-o target=1 ${outputFilePath110m}`
 ].join(" ")
 await mapshaper.runCommands(commandsAllFeatures110m);
-console.log(commandsAllFeatures110m)
 
 // Get countries from all polygon features
 const inputFilePathCountries110m = outputFilePath110m;
