fix(datacenters): preserve missing geometry and compare ring topology

This commit is contained in:
German Khokhlov committed 2026-11-03 15:02:42 +04:00
1 parent 4318f51de4
commit 04d626a4b1
3 files changed
+148 -13

No files matched your search

+52 -12
View File
@@ -4,7 +4,8 @@
Downloads a dated Geofabrik Norway extract, verifies its MD5, filters
telecom=data_center / building=data_center with osmium, and merges the result
into datacenters.geojsonl: features inside the Geofabrik Norway boundary are
added, updated or removed; every other line is kept byte-for-byte.
added, updated or removed once no longer tagged; every other line is kept
byte-for-byte.
Geofabrik keeps only its newest dated extract, so --extract defaults to the
current one and a recorded extract can be re-verified only while it is listed.
@@ -16,6 +17,7 @@ import math
from pathlib import Path
import shutil
import subprocess
import sys
import urllib.request
ROOT = Path(__file__).resolve().parent.parent
@@ -95,10 +97,15 @@ def osm_key(feature):
def export(osmium, pbf, cache, extract):
"""Exported features by key, and the keys of every tagged object in the extract."""
filtered = cache / (extract + '-datacenters.osm.pbf')
listed = cache / (extract + '-datacenters.opl')
exported = cache / (extract + '-datacenters.geojsonseq')
subprocess.run([osmium, 'tags-filter', '-O', '-o', str(filtered), str(pbf),
'nwr/telecom=data_center', 'nwr/building=data_center'], check=True)
subprocess.run([osmium, 'tags-filter', '-O', '-R', '-f', 'opl', '-o', str(listed), str(filtered),
'nwr/telecom=data_center', 'nwr/building=data_center'], check=True)
tagged = {line.split(' ', 1)[0] for line in listed.read_text(encoding='utf-8').splitlines() if line}
subprocess.run([osmium, 'export', '-O', '-f', 'geojsonseq', '--add-unique-id=type_id',
'-o', str(exported), str(filtered)], check=True)
fresh = {}
@@ -122,22 +129,48 @@ def export(osmium, pbf, cache, extract):
elif key in fresh or geometry['type'] == 'LineString':
continue # closed ways are taken from their area export
fresh[key] = {'geometry': geometry, 'tags': tags}
return fresh
return fresh, tagged
def same_geometry(old, new):
"""Same type and vertex set within tolerance; ring start and winding may differ."""
"""Same type and rings within tolerance; ring start, winding and part or hole order may differ."""
if old['type'] != new['type']:
return False
a, b = points(old['coordinates']), points(new['coordinates'])
if len(a) != len(b):
return False
def meters(p, q):
dx = (p[0] - q[0]) * 111320 * math.cos(math.radians(p[1]))
return math.hypot(dx, (p[1] - q[1]) * 110540)
return all(min(meters(p, q) for q in b) <= SAME_GEOMETRY_METERS for p in a) and \
all(min(meters(q, p) for p in a) <= SAME_GEOMETRY_METERS for q in b)
def same_line(a, b):
return len(a) == len(b) and all(meters(p, q) <= SAME_GEOMETRY_METERS for p, q in zip(a, b))
def same_ring(a, b):
a, b = a[:-1], b[:-1] # closed rings repeat their first vertex
return len(a) == len(b) and \
any(same_line(a[i:] + a[:i], c) for c in (b, b[::-1]) for i in range(len(a)))
def same_parts(a, b, same):
rest = list(b)
for part in a:
match = next((i for i, other in enumerate(rest) if same(part, other)), None)
if match is None:
return False
rest.pop(match)
return not rest
def same_polygon(a, b):
return same_ring(a[0], b[0]) and same_parts(a[1:], b[1:], same_ring)
a, b = old['coordinates'], new['coordinates']
if old['type'] == 'Point':
return same_line([a], [b])
if old['type'] == 'LineString':
return same_line(a, b)
if old['type'] == 'Polygon':
return same_polygon(a, b)
if old['type'] == 'MultiPolygon':
return same_parts(a, b, same_polygon)
return a == b
def feature_line(key, sequence, fresh):
@@ -148,9 +181,9 @@ def feature_line(key, sequence, fresh):
return json.dumps(feature, separators=(',', ':'), ensure_ascii=False)
def merge(lines, fresh, rings):
def merge(lines, fresh, rings, tagged):
out, seen, sequence = [], set(), 0
report = {'kept': 0, 'updated': [], 'removed': [], 'added': []}
report = {'kept': 0, 'updated': [], 'removed': [], 'added': [], 'unexported': []}
for line in lines:
feature = json.loads(line)
sequence = max(sequence, feature['properties']['id'])
@@ -165,6 +198,10 @@ def merge(lines, fresh, rings):
out.append(feature_line(key, feature['properties']['id'], fresh[key]))
report['updated'].append(key)
continue
if key in tagged: # still tagged, but osmium built no usable geometry: keep the old line
report['unexported'].append(key)
out.append(line)
continue
lon, lat = points(feature['geometry']['coordinates'])[0][:2]
if inside(rings, lon, lat):
report['removed'].append(key)
@@ -196,9 +233,12 @@ def main():
if actual != expected:
raise RuntimeError('MD5 mismatch for %s: %s != %s' % (pbf.name, actual, expected))
rings = read_poly(fetch(GEOFABRIK + 'norway.poly', args.cache / 'norway.poly'))
fresh = export(args.osmium, pbf, args.cache, extract)
fresh, tagged = export(args.osmium, pbf, args.cache, extract)
current = [line for line in DATASET.read_text(encoding='utf-8').splitlines() if line.strip()]
merged, report = merge(current, fresh, rings)
merged, report = merge(current, fresh, rings, tagged)
if report['unexported']:
print('warning: kept the old line of tagged features osmium did not export: %s'
% ', '.join(report['unexported']), file=sys.stderr)
payload = '\n'.join(merged) + '\n'
if args.verify:
if payload != DATASET.read_text(encoding='utf-8'):
@@ -0,0 +1,93 @@
import test from 'node:test';
import assert from 'node:assert/strict';
import { spawnSync } from 'node:child_process';
import { fileURLToPath } from 'node:url';
const SCRIPT = fileURLToPath(
new URL('../../scripts/build-datacenters-norway.py', import.meta.url),
);
function python(t, source) {
const result = spawnSync('python3', ['-c', source, SCRIPT], {
encoding: 'utf8',
});
if (result.error?.code === 'ENOENT') return t.skip('python3 is not installed');
assert.equal(result.status, 0, result.stderr);
}
test('Norway refresh compares rings, not vertex sets', (t) => {
python(
t,
String.raw`
import runpy, sys
same = runpy.run_path(sys.argv[1])['same_geometry']
poly = lambda *rings: {'type': 'Polygon', 'coordinates': list(rings)}
multi = lambda *polys: {'type': 'MultiPolygon', 'coordinates': [p['coordinates'] for p in polys]}
a, b, c, d, e = [0, 0], [4, 0], [4, 4], [0, 4], [2, 1]
assert not same(poly([a, b, c, d, e, a]), poly([a, b, e, c, d, a])), 'reordered edges'
assert same(poly([a, b, c, d, e, a]), poly([c, b, a, e, d, c])), 'ring start and winding'
assert same(poly([a, b, c, d, a]), poly([[0, 1e-7], b, c, d, [0, 1e-7]])), 'tolerance'
outer = [[0, 0], [9, 0], [9, 9], [0, 9], [0, 0]]
one, two = [[1, 1], [2, 1], [2, 2], [1, 1]], [[5, 5], [6, 5], [6, 6], [5, 5]]
assert same(poly(outer, one, two), poly(outer, two, one)), 'hole order'
assert not same(poly(outer, one), poly(one, outer)), 'outer and inner roles swapped'
assert same(multi(poly(one), poly(two)), multi(poly(two), poly(one))), 'part order'
p, q, r, s, u, v = [0, 0], [1, 0], [0, 1], [5, 5], [6, 5], [5, 6]
assert not same(multi(poly([p, q, r, p]), poly([s, u, v, s])),
multi(poly([p, q, s, p]), poly([r, u, v, r]))), 'vertices moved between parts'
`,
);
});
test('Norway refresh keeps a still-tagged feature osmium did not export', (t) => {
python(
t,
String.raw`
import json, runpy, sys
merge = runpy.run_path(sys.argv[1])['merge']
square = [[[10, 60], [10.01, 60], [10.01, 60.01], [10, 60.01], [10, 60]]]
def line(osm_id, sequence):
return json.dumps({'id': osm_id, 'type': 'Feature', 'geometry': {'type': 'Polygon', 'coordinates': square},
'properties': {'id': sequence, 'tags': {'building': 'data_center'},
'type': 'data_center', 'osm_id': osm_id}})
norway = [([(0, 0), (20, 0), (20, 80), (0, 80)], False)]
lines = [line(7, 1), line(8, 2)]
out, report = merge(lines, {}, norway, {'w7'})
assert out == [lines[0]], out
assert report['unexported'] == ['w7'] and report['removed'] == ['w8'], report
`,
);
});
test('Norway refresh lists tagged objects apart from the exported geometry', (t) => {
python(
t,
String.raw`
import json, runpy, subprocess, sys, tempfile
from pathlib import Path
export = runpy.run_path(sys.argv[1])['export']
tags = {'building': 'data_center'}
ring = [[10, 60], [10.01, 60], [10.01, 60.01], [10, 60]]
features = [
{'id': 'n1', 'properties': {'telecom': 'data_center'}, 'geometry': {'type': 'Point', 'coordinates': [10, 60]}},
{'id': 'w5', 'properties': tags, 'geometry': {'type': 'LineString', 'coordinates': ring}},
{'id': 'a10', 'properties': tags, 'geometry': {'type': 'MultiPolygon', 'coordinates': [[ring]]}},
{'id': 'w7', 'properties': tags, 'geometry': {'type': 'LineString', 'coordinates': ring[:3]}},
]
outputs = {
'.opl': 'n1 v1 Ttelecom=data_center x10 y60\nw5 v1 Tbuilding=data_center Nn2,n3,n4,n2\n'
'w7 v1 Tbuilding=data_center Nn2,n3,n4\nr9 v1 Ttelecom=data_center Mw8@outer\n',
'.geojsonseq': ''.join(chr(30) + json.dumps({'type': 'Feature', **f}) + '\n' for f in features),
}
def run(args, check=False):
target = Path(args[args.index('-o') + 1])
target.write_text(outputs.get(target.suffix, ''), encoding='utf-8')
subprocess.run = run
with tempfile.TemporaryDirectory() as cache:
fresh, tagged = export('osmium', Path(cache) / 'norway-test.osm.pbf', Path(cache), 'norway-test')
assert sorted(fresh) == ['n1', 'w5'], fresh
assert fresh['w5']['geometry'] == {'type': 'Polygon', 'coordinates': [ring]}, fresh['w5']
assert tagged == {'n1', 'w5', 'w7', 'r9'}, tagged
`,
);
});
+3 -1
View File
@@ -27,7 +27,9 @@ reflects OSM as of the extract below, everything else the original snapshot.
2026-10-05T20:21:35Z, MD5 `1d2cfdeb98980fce11b8891d0f9146c0`
- Filter: `osmium tags-filter nwr/telecom=data_center nwr/building=data_center`
- Contact tags are dropped as above; unchanged features keep their original
line, so geometry is only rewritten when tags or vertices actually changed.
line, so geometry is only rewritten when tags or rings actually changed.
A feature is removed only once it is no longer tagged; one still tagged but
without an exportable geometry keeps its line and is reported as unexported.
- Refresh: `python3 scripts/build-datacenters-norway.py --cache <dir>` uses
Geofabrik's current Norway extract and prints the dated name it used; record
that name, its OSM timestamp and MD5 here. Requires `osmium-tool`.