1"""EUDR Parcel Polygons Austria – österreichische Grundstücksnummer -> EUDR-taugliches GeoJSON-Polygon.
2
3Datenquelle: BEV – Bundesamt für Eich- und Vermessungswesen, CC BY 4.0 (kataster.bev.gv.at).
4Alle Module (Kachel-Decoder, Geometrie, Kataster-Client, Kernlogik, Apify-Einstieg) in einer Datei.
5"""
6
7from __future__ import annotations
8
9import asyncio
10import gzip
11import json
12import math
13import re
14import struct
15import sys
16import threading
17import time
18import urllib.error
19import urllib.parse
20import urllib.request
21from collections import Counter, defaultdict
22from dataclasses import dataclass, field
23from datetime import datetime, timezone
24
25from apify import Actor
26
27
28geo = mvt = sys.modules[__name__]
29
30
31
32
33
34POINT, LINESTRING, POLYGON = 1, 2, 3
35
36
37@dataclass
38class Feature:
39 id: int | None
40 type: int
41 props: dict
42 rings: list[list[tuple[int, int]]]
43
44
45@dataclass
46class Layer:
47 name: str
48 extent: int = 4096
49 features: list[Feature] = field(default_factory=list)
50
51
52class _Reader:
53 __slots__ = ("buf", "pos", "end")
54
55 def __init__(self, buf: bytes, start: int = 0, end: int | None = None) -> None:
56 self.buf = buf
57 self.pos = start
58 self.end = len(buf) if end is None else end
59
60 def varint(self) -> int:
61 result, shift = 0, 0
62 while True:
63 b = self.buf[self.pos]
64 self.pos += 1
65 result |= (b & 0x7F) << shift
66 if not b & 0x80:
67 return result
68 shift += 7
69
70 def fields(self):
71 """Liefert (Feldnummer, Wiretyp, Wert) – bei Längen-Feldern (start, end)."""
72 while self.pos < self.end:
73 key = self.varint()
74 num, wtype = key >> 3, key & 7
75 if wtype == 0:
76 yield num, wtype, self.varint()
77 elif wtype == 2:
78 length = self.varint()
79 start = self.pos
80 self.pos += length
81 yield num, wtype, (start, start + length)
82 elif wtype == 1:
83 start = self.pos
84 self.pos += 8
85 yield num, wtype, (start, start + 8)
86 elif wtype == 5:
87 start = self.pos
88 self.pos += 4
89 yield num, wtype, (start, start + 4)
90 else:
91 raise ValueError(f"Unbekannter Protobuf-Wiretyp {wtype}")
92
93
94def _zigzag(n: int) -> int:
95 return (n >> 1) ^ -(n & 1)
96
97
98def _signed64(n: int) -> int:
99 return n - (1 << 64) if n >= (1 << 63) else n
100
101
102def _packed(buf: bytes, span: tuple[int, int]) -> list[int]:
103 r = _Reader(buf, *span)
104 out = []
105 while r.pos < r.end:
106 out.append(r.varint())
107 return out
108
109
110def _value(buf: bytes, span: tuple[int, int]):
111 val = None
112 for num, wtype, v in _Reader(buf, *span).fields():
113 if num == 1:
114 val = buf[v[0]:v[1]].decode("utf-8")
115 elif num == 2:
116 val = struct.unpack("<f", buf[v[0]:v[1]])[0]
117 elif num == 3:
118 val = struct.unpack("<d", buf[v[0]:v[1]])[0]
119 elif num == 4:
120 val = _signed64(v)
121 elif num == 5:
122 val = v
123 elif num == 6:
124 val = _zigzag(v)
125 elif num == 7:
126 val = bool(v)
127 return val
128
129
130def _geometry(cmds: list[int], gtype: int) -> list[list[tuple[int, int]]]:
131 rings: list[list[tuple[int, int]]] = []
132 cur: list[tuple[int, int]] | None = None
133 x = y = 0
134 i = 0
135 while i < len(cmds):
136 c = cmds[i]
137 i += 1
138 cid, count = c & 7, c >> 3
139 if cid == 7:
140 if cur:
141 rings.append(cur)
142 cur = None
143 continue
144 for _ in range(count):
145 x += _zigzag(cmds[i])
146 y += _zigzag(cmds[i + 1])
147 i += 2
148 if cid == 1:
149 if cur:
150 rings.append(cur)
151 cur = [(x, y)]
152 else:
153 if cur is None:
154 cur = [(x, y)]
155 else:
156 cur.append((x, y))
157 if cur and gtype != POLYGON:
158 rings.append(cur)
159 return rings
160
161
162def decode(buf: bytes, layers: set[str] | None = None) -> dict[str, Layer]:
163 """Dekodiert eine Kachel. Mit `layers` werden nur diese Ebenen gelesen."""
164 result: dict[str, Layer] = {}
165 for num, _, span in _Reader(buf).fields():
166 if num != 3:
167 continue
168 name = None
169 extent = 4096
170 keys: list[str] = []
171 value_spans = []
172 feature_spans = []
173 for n2, _, v in _Reader(buf, *span).fields():
174 if n2 == 1:
175 name = buf[v[0]:v[1]].decode("utf-8")
176 elif n2 == 2:
177 feature_spans.append(v)
178 elif n2 == 3:
179 keys.append(buf[v[0]:v[1]].decode("utf-8"))
180 elif n2 == 4:
181 value_spans.append(v)
182 elif n2 == 5:
183 extent = v
184 if name is None or (layers is not None and name not in layers):
185 continue
186 values = [_value(buf, s) for s in value_spans]
187 layer = Layer(name=name, extent=extent)
188 for fspan in feature_spans:
189 fid, gtype, tags, cmds = None, 0, [], []
190 for n3, _, v in _Reader(buf, *fspan).fields():
191 if n3 == 1:
192 fid = v
193 elif n3 == 2:
194 tags = _packed(buf, v)
195 elif n3 == 3:
196 gtype = v
197 elif n3 == 4:
198 cmds = _packed(buf, v)
199 props = {keys[tags[k]]: values[tags[k + 1]] for k in range(0, len(tags) - 1, 2)}
200 layer.features.append(Feature(fid, gtype, props, _geometry(cmds, gtype)))
201 result[name] = layer
202 return result
203
204
205
206
207Point = tuple[float, float]
208Ring = list[Point]
209
210
211
212
213def ring_area(ring: Ring) -> float:
214 """Vorzeichenbehaftete Fläche (Gaußsche Trapezformel, ohne Schlusspunkt)."""
215 s = 0.0
216 n = len(ring)
217 for i in range(n):
218 x1, y1 = ring[i]
219 x2, y2 = ring[(i + 1) % n]
220 s += x1 * y2 - x2 * y1
221 return s / 2.0
222
223
224def _dedupe(ring: Ring) -> Ring:
225 out: Ring = []
226 for p in ring:
227 if not out or out[-1] != p:
228 out.append(p)
229 while len(out) > 1 and out[0] == out[-1]:
230 out.pop()
231 return out
232
233
234def point_in_ring(pt: Point, ring: Ring) -> bool:
235 x, y = pt
236 inside = False
237 n = len(ring)
238 j = n - 1
239 for i in range(n):
240 xi, yi = ring[i]
241 xj, yj = ring[j]
242 if (yi > y) != (yj > y):
243 xcross = xi + (y - yi) * (xj - xi) / (yj - yi)
244 if x < xcross:
245 inside = not inside
246 j = i
247 return inside
248
249
250
251
252def clip_ring(ring: Ring, xmin: float, ymin: float, xmax: float, ymax: float) -> Ring:
253 """Sutherland-Hodgman-Zuschnitt eines Rings auf ein Rechteck (ganzzahlig gerundet)."""
254
255 def clip(points: Ring, inside, intersect) -> Ring:
256 if not points:
257 return []
258 out: Ring = []
259 prev = points[-1]
260 prev_in = inside(prev)
261 for cur in points:
262 cur_in = inside(cur)
263 if cur_in:
264 if not prev_in:
265 out.append(intersect(prev, cur))
266 out.append(cur)
267 elif prev_in:
268 out.append(intersect(prev, cur))
269 prev, prev_in = cur, cur_in
270 return out
271
272 def ix(x0):
273 def f(a, b):
274 t = (x0 - a[0]) / (b[0] - a[0])
275 return (x0, round(a[1] + t * (b[1] - a[1])))
276 return f
277
278 def iy(y0):
279 def f(a, b):
280 t = (y0 - a[1]) / (b[1] - a[1])
281 return (round(a[0] + t * (b[0] - a[0])), y0)
282 return f
283
284 pts = list(ring)
285 pts = clip(pts, lambda p: p[0] >= xmin, ix(xmin))
286 pts = clip(pts, lambda p: p[0] <= xmax, ix(xmax))
287 pts = clip(pts, lambda p: p[1] >= ymin, iy(ymin))
288 pts = clip(pts, lambda p: p[1] <= ymax, iy(ymax))
289 return _dedupe(pts)
290
291
292
293
294def _cluster(values: list[float], tol: float, grid: int) -> dict[float, float]:
295 """Werte, die näher als `tol` beieinander liegen, auf einen Vertreter einrasten.
296 Exakte Rastervielfache (Kachelecken) haben Vorrang."""
297 mapping: dict[float, float] = {}
298 vals = sorted(set(values))
299 group: list[float] = []
300
301 def flush():
302 if not group:
303 return
304 rep = next((v for v in group if v % grid == 0), group[len(group) // 2])
305 for v in group:
306 mapping[v] = rep
307
308 for v in vals:
309 if group and v - group[-1] > tol:
310 flush()
311 group = []
312 group.append(v)
313 flush()
314 return mapping
315
316
317def _simplify(ring: Ring, tol: float) -> Ring:
318 """Entfernt (fast) kollineare Zwischenpunkte, z. B. die Schnittpunkte an Kachelrändern."""
319 changed = True
320 pts = list(ring)
321 while changed and len(pts) > 3:
322 changed = False
323 out: Ring = []
324 n = len(pts)
325 for i in range(n):
326 a = out[-1] if out else pts[i - 1]
327 p = pts[i]
328 b = pts[(i + 1) % n]
329 dx, dy = b[0] - a[0], b[1] - a[1]
330 length = math.hypot(dx, dy)
331 if length == 0:
332 changed = True
333 continue
334 dist = abs((p[0] - a[0]) * dy - (p[1] - a[1]) * dx) / length
335
336 proj = ((p[0] - a[0]) * dx + (p[1] - a[1]) * dy) / (length * length)
337 if dist <= tol and 0 < proj < 1:
338 changed = True
339 continue
340 out.append(p)
341 pts = _dedupe(out)
342 return pts
343
344
345def assemble(rings: list[Ring], grid: int, snap_tol: float = 8, simplify_tol: float = 1.5) -> list[tuple[Ring, list[Ring]]]:
346 """Setzt zugeschnittene Ringe (globale Koordinaten) zu Polygonen zusammen.
347
348 Rückgabe: Liste von (Außenring, [Löcher]); Außenringe haben positive Fläche.
349 """
350
351 on_v: dict[float, list[float]] = defaultdict(list)
352 on_h: dict[float, list[float]] = defaultdict(list)
353 for ring in rings:
354 for x, y in ring:
355 if x % grid == 0:
356 on_v[x].append(y)
357 if y % grid == 0:
358 on_h[y].append(x)
359 vmap = {x: _cluster(ys, snap_tol, grid) for x, ys in on_v.items()}
360 hmap = {y: _cluster(xs, snap_tol, grid) for y, xs in on_h.items()}
361
362 def snap(p: Point) -> Point:
363 x, y = p
364 if x % grid == 0 and y % grid == 0:
365 return p
366 if x % grid == 0:
367 return (x, vmap[x][y])
368 if y % grid == 0:
369 return (hmap[y][x], y)
370 return p
371
372 snapped = []
373 for ring in rings:
374 r = _dedupe([snap(p) for p in ring])
375 if len(r) >= 3:
376 snapped.append(r)
377
378
379 v_pts: dict[float, set[float]] = defaultdict(set)
380 h_pts: dict[float, set[float]] = defaultdict(set)
381 for ring in snapped:
382 for x, y in ring:
383 if x % grid == 0:
384 v_pts[x].add(y)
385 if y % grid == 0:
386 h_pts[y].add(x)
387
388 edges: list[tuple[Point, Point]] = []
389 for ring in snapped:
390 n = len(ring)
391 for i in range(n):
392 a, b = ring[i], ring[(i + 1) % n]
393 if a == b:
394 continue
395 if a[0] == b[0] and a[0] % grid == 0:
396 lo, hi = sorted((a[1], b[1]))
397 cuts = sorted(y for y in v_pts[a[0]] if lo < y < hi)
398 if a[1] > b[1]:
399 cuts.reverse()
400 seq = [a] + [(a[0], y) for y in cuts] + [b]
401 elif a[1] == b[1] and a[1] % grid == 0:
402 lo, hi = sorted((a[0], b[0]))
403 cuts = sorted(x for x in h_pts[a[1]] if lo < x < hi)
404 if a[0] > b[0]:
405 cuts.reverse()
406 seq = [a] + [(x, a[1]) for x in cuts] + [b]
407 else:
408 seq = [a, b]
409 for k in range(len(seq) - 1):
410 edges.append((seq[k], seq[k + 1]))
411
412
413 count = Counter(edges)
414 remaining: list[tuple[Point, Point]] = []
415 seen = set()
416 for (a, b) in count:
417 key = (a, b) if a <= b else (b, a)
418 if key in seen:
419 continue
420 seen.add(key)
421 net = count.get((a, b), 0) - count.get((b, a), 0)
422 if net > 0:
423 remaining.append((a, b))
424 elif net < 0:
425 remaining.append((b, a))
426
427
428 outgoing: dict[Point, list[Point]] = defaultdict(list)
429 for a, b in remaining:
430 outgoing[a].append(b)
431
432 result_rings: list[Ring] = []
433 while True:
434 start = next((p for p, lst in outgoing.items() if lst), None)
435 if start is None:
436 break
437 ring = [start]
438 prev, cur = start, outgoing[start].pop()
439 guard = 0
440 while cur != start and guard < 1_000_000:
441 guard += 1
442 ring.append(cur)
443 cands = outgoing.get(cur)
444 if not cands:
445 break
446 if len(cands) == 1:
447 nxt = cands.pop()
448 else:
449 din = (cur[0] - prev[0], cur[1] - prev[1])
450
451 def turn(c):
452 dout = (c[0] - cur[0], c[1] - cur[1])
453 cross = din[0] * dout[1] - din[1] * dout[0]
454 dot = din[0] * dout[0] + din[1] * dout[1]
455 return math.atan2(cross, dot)
456
457 nxt = max(cands, key=turn)
458 cands.remove(nxt)
459 prev, cur = cur, nxt
460 ring = _simplify(_dedupe(ring), simplify_tol)
461 if len(ring) >= 3 and abs(ring_area(ring)) > 0:
462 result_rings.append(ring)
463
464
465 shells = [r for r in result_rings if ring_area(r) > 0]
466 holes = [r for r in result_rings if ring_area(r) < 0]
467 polygons: list[tuple[Ring, list[Ring]]] = [(s, []) for s in shells]
468 for h in holes:
469 a, b = h[0], h[1]
470 probe = ((a[0] + b[0]) / 2, (a[1] + b[1]) / 2)
471 owners = [i for i, (s, _) in enumerate(polygons) if point_in_ring(probe, s)]
472 if owners:
473 i = min(owners, key=lambda k: ring_area(polygons[k][0]))
474 polygons[i][1].append(h)
475 polygons.sort(key=lambda p: -ring_area(p[0]))
476 return polygons
477
478
479
480
481def lonlat_to_tile(lon: float, lat: float, z: int) -> tuple[float, float]:
482 n = 2 ** z
483 x = (lon + 180.0) / 360.0 * n
484 y = (1.0 - math.asinh(math.tan(math.radians(lat))) / math.pi) / 2.0 * n
485 return x, y
486
487
488def global_to_lonlat(X: float, Y: float, z: int, extent: int) -> tuple[float, float]:
489 n = extent * (2 ** z)
490 lon = X / n * 360.0 - 180.0
491 lat = math.degrees(math.atan(math.sinh(math.pi * (1.0 - 2.0 * Y / n))))
492 return lon, lat
493
494
495_A = 6378137.0
496_E2 = 0.00669437999014
497
498
499def area_m2(polygons_lonlat: list[list[list[tuple[float, float]]]]) -> float:
500 """Fläche auf dem WGS84-Ellipsoid (lokale flächentreue Näherung, Fehler << 0,1 % bei Grundstücken)."""
501 pts = [p for poly in polygons_lonlat for ring in poly for p in ring]
502 if not pts:
503 return 0.0
504 lat0 = math.radians(sum(p[1] for p in pts) / len(pts))
505 lon0 = sum(p[0] for p in pts) / len(pts)
506 m0 = _A * (1 - _E2) / (1 - _E2 * math.sin(lat0) ** 2) ** 1.5
507
508 def proj(p):
509 lat = math.radians(p[1])
510 n = _A / math.sqrt(1 - _E2 * math.sin(lat) ** 2)
511 return (math.radians(p[0] - lon0) * n * math.cos(lat), (lat - lat0) * m0)
512
513 total = 0.0
514 for poly in polygons_lonlat:
515 for k, ring in enumerate(poly):
516 a = abs(ring_area([proj(p) for p in ring]))
517 total += a if k == 0 else -a
518 return total
519
520
521
522
523BASE_URL = "https://kataster.bev.gv.at"
524USER_AGENT = "EUDR-Parcel-Polygons-Austria/0.1 (Apify Actor; Datenquelle BEV, CC BY 4.0)"
525
526
527class BevClient:
528 def __init__(self, base_url: str = BASE_URL, timeout: float = 30.0, retries: int = 3, pause: float = 0.05) -> None:
529 self.base_url = base_url.rstrip("/")
530 self.timeout = timeout
531 self.retries = retries
532 self.pause = pause
533
534 def _get(self, path: str, accept: str) -> bytes | None:
535 url = self.base_url + path
536 last_err: Exception | None = None
537 for attempt in range(self.retries + 1):
538 req = urllib.request.Request(url, headers={
539 "User-Agent": USER_AGENT,
540 "Accept": accept,
541 "Accept-Encoding": "gzip",
542 })
543 try:
544 with urllib.request.urlopen(req, timeout=self.timeout) as resp:
545 data = resp.read()
546 if resp.headers.get("Content-Encoding", "").lower() == "gzip" or data[:2] == b"\x1f\x8b":
547 data = gzip.decompress(data)
548 if self.pause:
549 time.sleep(self.pause)
550 return data
551 except urllib.error.HTTPError as e:
552 if e.code in (204, 404):
553 return None
554 last_err = e
555 if e.code not in (429, 500, 502, 503, 504):
556 raise
557 except (urllib.error.URLError, TimeoutError, ConnectionError) as e:
558 last_err = e
559 time.sleep(min(2 ** attempt, 10))
560 raise RuntimeError(f"BEV-Dienst nicht erreichbar ({url}): {last_err}")
561
562 def get_json(self, path: str) -> dict | None:
563 data = self._get(path, "application/json")
564 return None if data is None else json.loads(data.decode("utf-8"))
565
566 def search_parcel(self, kg: str, gnr: str) -> dict | None:
567 term = urllib.parse.quote(f"{kg}-{gnr}", safe="")
568 return self.get_json(f"/api/search?term={term}&layers=gst")
569
570 def parcel_details(self, parcel_id: int) -> dict | None:
571 return self.get_json(f"/api/gst/{int(parcel_id)}")
572
573 def tile(self, z: int, x: int, y: int) -> bytes | None:
574 return self._get(f"/tiles/kataster/{z}/{x}/{y}.pbf", "application/x-protobuf")
575
576
577
578
579ZOOM = 16
580ATTRIBUTION = "Datenquelle: BEV – Bundesamt für Eich- und Vermessungswesen, CC BY 4.0 (kataster.bev.gv.at)"
581EUDR_POLYGON_THRESHOLD_HA = 4.0
582
583
584
585NUTZUNG_NAMES = {
586 "40": "Dauerkulturanlagen oder Erwerbsgärten", "41": "Gebäude", "42": "Parkplätze",
587 "48": "Äcker, Wiesen oder Weiden", "52": "Gärten", "53": "Weingärten", "54": "Alpen",
588 "55": "Krummholzflächen", "56": "Wälder", "57": "Verbuschte Flächen", "58": "Forststraßen",
589 "59": "Fließende Gewässer", "60": "Stehende Gewässer", "61": "Feuchtgebiete",
590 "62": "Vegetationsarme Flächen", "63": "Betriebsflächen", "64": "Gewässerrandflächen",
591 "65": "Verkehrsrandflächen", "72": "Friedhöfe", "73": "Rechtlich nicht Wald", "74": "Rechtlich Wald",
592 "77": "Rechtlich Weingarten", "78": "Rechtlich kein Weingarten", "83": "Gebäudenebenflächen",
593 "84": "Abbauflächen, Halden und Deponien", "87": "Fels und Geröllflächen", "88": "Gletscher",
594 "92": "Schienenverkehrsanlagen", "95": "Straßenverkehrsanlagen", "96": "Freizeitflächen",
595}
596FOREST_CODES = ("55", "56", "74")
597
598
599def _nutzung_name(n: dict) -> str | None:
600 name = n.get("nutzung")
601 if name is None or str(name).strip() in ("", "0"):
602 return NUTZUNG_NAMES.get(str(n.get("ns")), None if name is None else str(name))
603 return name
604
605_PARCEL_RE = re.compile(r"^\s*(?:KG\.?\s*)?(\d{1,5})\s*[-;,:\s]\s*(\.?\d+(?:/\d+)?)\s*$", re.IGNORECASE)
606
607
608class InputError(ValueError):
609 pass
610
611
612@dataclass
613class ParcelRef:
614 kg: str
615 gnr: str
616 reference: str | None = None
617 raw: str = ""
618
619 @property
620 def label(self) -> str:
621 return f"{self.kg}-{self.gnr}"
622
623
624def parse_parcel(item) -> ParcelRef:
625 """Akzeptiert '72106-116/81', '72106 116/81', '72106;.116' oder {'kg':..,'gnr':..,'reference':..}."""
626 if isinstance(item, dict):
627 kg = str(item.get("kg") or item.get("kgnr") or item.get("katastralgemeinde") or "").strip()
628 gnr = str(item.get("gnr") or item.get("grundstuecksnummer") or "").strip().replace(" ", "")
629 ref = item.get("reference") or item.get("referenz")
630 raw = f"{kg}-{gnr}"
631 if not re.fullmatch(r"\d{1,5}", kg) or not re.fullmatch(r"\.?\d+(?:/\d+)?", gnr):
632 raise InputError(f"Ungültige Angabe: {item!r}")
633 return ParcelRef(kg.zfill(5), gnr, str(ref) if ref else None, raw)
634 raw = str(item)
635 m = _PARCEL_RE.match(raw)
636 if not m:
637 raise InputError(
638 f"'{raw}' nicht erkannt. Erwartet: Katastralgemeindenummer und Grundstücksnummer, z. B. 72106-116/81"
639 )
640 return ParcelRef(m.group(1).zfill(5), m.group(2), None, raw)
641
642
643def _tile_range(bbox: tuple[float, float, float, float], z: int) -> list[tuple[int, int]]:
644 minlon, minlat, maxlon, maxlat = bbox
645 x0, y0 = geo.lonlat_to_tile(minlon, maxlat, z)
646 x1, y1 = geo.lonlat_to_tile(maxlon, minlat, z)
647 return [(x, y) for x in range(int(x0), int(x1) + 1) for y in range(int(y0), int(y1) + 1)]
648
649
650def _bbox_of(geometry: dict) -> tuple[float, float, float, float] | None:
651 coords = []
652
653 def walk(c):
654 if c and isinstance(c[0], (int, float)):
655 coords.append(c)
656 else:
657 for cc in c:
658 walk(cc)
659
660 if not geometry:
661 return None
662 walk(geometry.get("coordinates", []))
663 if not coords:
664 return None
665 lons = [c[0] for c in coords]
666 lats = [c[1] for c in coords]
667 return min(lons), min(lats), max(lons), max(lats)
668
669
670def _interior_point(polygons: list[list[list[tuple[float, float]]]]) -> tuple[float, float]:
671 """Ein Punkt, der sicher im (größten) Polygon liegt – für die EUDR-Punktangabe."""
672 shell = polygons[0][0]
673 holes = polygons[0][1:]
674 ys = sorted(p[1] for p in shell)
675 best = None
676 for frac in (0.5, 0.35, 0.65, 0.2, 0.8):
677 y = ys[0] + (ys[-1] - ys[0]) * frac
678 xs = []
679 for ring in [shell] + holes:
680 n = len(ring)
681 for i in range(n):
682 (x1, y1), (x2, y2) = ring[i], ring[(i + 1) % n]
683 if (y1 > y) != (y2 > y):
684 xs.append(x1 + (y - y1) * (x2 - x1) / (y2 - y1))
685 xs.sort()
686 for i in range(0, len(xs) - 1, 2):
687 w = xs[i + 1] - xs[i]
688 if best is None or w > best[0]:
689 best = (w, ((xs[i] + xs[i + 1]) / 2, y))
690 if best:
691 break
692 if best is None:
693 return shell[0]
694 return best[1]
695
696
697class ParcelResolver:
698 def __init__(self, client, zoom: int = ZOOM, max_tiles: int = 400, precision: int = 7) -> None:
699 self.client = client
700 self.zoom = zoom
701 self.max_tiles = max_tiles
702 self.precision = precision
703 self._tiles: dict[tuple[int, int], dict] = {}
704 self._lock = threading.Lock()
705
706
707 def _tile(self, x: int, y: int) -> dict:
708 key = (x, y)
709 with self._lock:
710 if key in self._tiles:
711 return self._tiles[key]
712 data = self.client.tile(self.zoom, x, y)
713 layers = mvt.decode(data, {"gst", "kg", "pg"}) if data else {}
714 with self._lock:
715 self._tiles[key] = layers
716 if len(self._tiles) > 3000:
717 self._tiles.pop(next(iter(self._tiles)))
718 return layers
719
720
721 def _find(self, ref: ParcelRef) -> tuple[dict | None, list[dict]]:
722 res = self.client.search_parcel(ref.kg, ref.gnr) or {}
723 feats = ((res.get("data") or {}).get("features")) or []
724 feats = [f for f in feats if str((f.get("properties") or {}).get("kg", "")).zfill(5) == ref.kg]
725 exact = [f for f in feats if (f.get("properties") or {}).get("gnr") == ref.gnr]
726 return (exact[0] if exact else None), feats
727
728 def _find_in_tiles(self, ref: ParcelRef, siblings: list[dict]) -> tuple[int | None, tuple | None]:
729 """Plan B, wenn die Suche zu viele Treffer liefert: rund um die Nachbar-Grundstücke
730 mit derselben Stammnummer in den Kacheln nachsehen."""
731 boxes = [b for b in (_bbox_of(f.get("geometry")) for f in siblings) if b]
732 if not boxes:
733 return None, None
734 bbox = (min(b[0] for b in boxes), min(b[1] for b in boxes), max(b[2] for b in boxes), max(b[3] for b in boxes))
735 tiles = _tile_range(bbox, self.zoom)
736 if len(tiles) > self.max_tiles:
737 return None, None
738 for tx, ty in tiles:
739 layer = self._tile(tx, ty).get("gst")
740 if not layer:
741 continue
742 for f in layer.features:
743 if f.props.get("gnr") == ref.gnr and str(f.props.get("kg", "")).zfill(5) == ref.kg:
744 xs = [p[0] for r in f.rings for p in r]
745 ys = [p[1] for r in f.rings for p in r]
746 e = layer.extent
747 lon0, lat0 = geo.global_to_lonlat(tx * e + min(xs), ty * e + max(ys), self.zoom, e)
748 lon1, lat1 = geo.global_to_lonlat(tx * e + max(xs), ty * e + min(ys), self.zoom, e)
749
750 return f.id, (lon0, lat0, lon1, lat1)
751 return None, None
752
753 def _municipality_at(self, lonlat: tuple[float, float]) -> str | None:
754 fx, fy = geo.lonlat_to_tile(lonlat[0], lonlat[1], self.zoom)
755 tx, ty = int(fx), int(fy)
756 layer = self._tile(tx, ty).get("pg")
757 if not layer:
758 return None
759 pt = ((fx - tx) * layer.extent, (fy - ty) * layer.extent)
760 for f in layer.features:
761 if f.type == mvt.POLYGON and sum(geo.point_in_ring(pt, r) for r in f.rings) % 2 == 1:
762 return f.props.get("pg")
763 return None
764
765
766 def resolve(self, ref: ParcelRef) -> dict:
767 now = datetime.now(timezone.utc).isoformat(timespec="seconds")
768 base = {
769 "input": ref.raw or ref.label,
770 "reference": ref.reference,
771 "kg": ref.kg,
772 "gnr": ref.gnr,
773 "retrievedAt": now,
774 "attribution": ATTRIBUTION,
775 }
776
777 hit, candidates = self._find(ref)
778 parcel_id, bbox = None, None
779 if hit:
780 parcel_id = hit.get("id")
781 bbox = _bbox_of(hit.get("geometry"))
782 else:
783 stem = ref.gnr.split("/")[0]
784 siblings = [f for f in candidates if str(f["properties"].get("gnr", "")).split("/")[0] == stem]
785 parcel_id, bbox = self._find_in_tiles(ref, siblings)
786 if parcel_id is None or bbox is None:
787 hint = ""
788 alt = "." + ref.gnr if not ref.gnr.startswith(".") else ref.gnr[1:]
789 if any(f["properties"].get("gnr") == alt for f in candidates):
790 hint = f" Meinten Sie {ref.kg}-{alt}? (Punkt = Baufläche)"
791 return {**base, "status": "NICHT_GEFUNDEN",
792 "message": f"Grundstück {ref.label} wurde im Kataster nicht gefunden.{hint}"}
793
794 details = (self.client.parcel_details(parcel_id) or {}).get("properties") or {}
795
796 land_use = [
797 {"nutzung": _nutzung_name(n), "code": str(n.get("ns")), "flaecheM2": n.get("fl")}
798 for n in details.get("nutzungen") or []
799 if n.get("gnr") in (None, ref.gnr)
800 ]
801 official_m2 = sum(int(n["flaecheM2"] or 0) for n in land_use) or None
802
803
804 pad = 2e-5
805 tiles = _tile_range((bbox[0] - pad, bbox[1] - pad, bbox[2] + pad, bbox[3] + pad), self.zoom)
806 if len(tiles) > self.max_tiles:
807 return {**base, "status": "FEHLER", "parcelId": parcel_id,
808 "message": f"Grundstück zu groß ({len(tiles)} Kacheln > Limit {self.max_tiles}). Limit im Input erhöhen."}
809
810 pieces: list[list[tuple[int, int]]] = []
811 extent = None
812 kg_name = None
813 pg_names: set[str] = set()
814 rstatus = details.get("rstatus")
815
816
817 queue = list(tiles)
818 visited: set[tuple[int, int]] = set()
819 while queue:
820 tx, ty = queue.pop(0)
821 if (tx, ty) in visited:
822 continue
823 visited.add((tx, ty))
824 if len(visited) > self.max_tiles:
825 return {**base, "status": "FEHLER", "parcelId": parcel_id,
826 "message": f"Grundstück zu groß (> {self.max_tiles} Kacheln). Limit im Input erhöhen."}
827 layers = self._tile(tx, ty)
828 if kg_name is None and layers.get("kg") and layers["kg"].features:
829 for f in layers["kg"].features:
830 if str(f.props.get("kg_nr", "")).zfill(5) == ref.kg:
831 kg_name = f.props.get("kg")
832 if layers.get("pg"):
833 pg_names.update(f.props.get("pg") for f in layers["pg"].features if f.props.get("pg"))
834 gst = layers.get("gst")
835 if not gst:
836 continue
837 e = gst.extent
838 extent = extent or e
839 for f in gst.features:
840 if f.type != mvt.POLYGON:
841 continue
842 if f.id != parcel_id and not (f.props.get("gnr") == ref.gnr and str(f.props.get("kg", "")).zfill(5) == ref.kg):
843 continue
844 rstatus = rstatus or f.props.get("rstatus")
845 x0, y0, x1, y1 = tx * e, ty * e, (tx + 1) * e, (ty + 1) * e
846 for ring in f.rings:
847 g = [(x0 + x, y0 + y) for x, y in ring]
848 c = geo.clip_ring(g, x0, y0, x1, y1)
849 if len(c) < 3:
850 continue
851 pieces.append(c)
852 for (nx, ny), touches in (((tx - 1, ty), any(p[0] == x0 for p in c)),
853 ((tx + 1, ty), any(p[0] == x1 for p in c)),
854 ((tx, ty - 1), any(p[1] == y0 for p in c)),
855 ((tx, ty + 1), any(p[1] == y1 for p in c))):
856 if touches and (nx, ny) not in visited:
857 queue.append((nx, ny))
858
859 if not pieces or extent is None:
860 return {**base, "status": "FEHLER", "parcelId": parcel_id,
861 "message": "Grundstück gefunden, aber keine Umrissdaten in den Kartenkacheln."}
862
863 polys = geo.assemble(pieces, extent)
864 if not polys:
865 return {**base, "status": "FEHLER", "parcelId": parcel_id,
866 "message": "Umriss konnte nicht zusammengesetzt werden."}
867
868
869 out_polys: list[list[list[tuple[float, float]]]] = []
870 for shell, holes in polys:
871 rings = []
872 for k, ring in enumerate([shell] + holes):
873 ll = [geo.global_to_lonlat(x, y, self.zoom, extent) for x, y in ring]
874 a = geo.ring_area(ll)
875 if (k == 0 and a < 0) or (k > 0 and a > 0):
876 ll.reverse()
877 rings.append(ll)
878 out_polys.append(rings)
879
880 area = geo.area_m2(out_polys)
881 p = self.precision
882
883 def fmt(ring):
884 pts = [[round(x, p), round(y, p)] for x, y in ring]
885 return pts + [pts[0]]
886
887 if len(out_polys) == 1:
888 geometry = {"type": "Polygon", "coordinates": [fmt(r) for r in out_polys[0]]}
889 else:
890 geometry = {"type": "MultiPolygon", "coordinates": [[fmt(r) for r in poly] for poly in out_polys]}
891
892 ipt = _interior_point(out_polys)
893
894
895 municipality = self._municipality_at(ipt) or (next(iter(pg_names)) if len(pg_names) == 1 else None)
896 deviation = None
897 status, message = "OK", "Polygon erstellt."
898 if official_m2:
899 deviation = (area - official_m2) / official_m2 * 100
900 tolerance = max(3.0, 100 * 25 / official_m2)
901 if abs(deviation) > tolerance:
902 status = "PRUEFEN"
903 dev_txt = f"{deviation:+.1f}".replace(".", ",")
904 message = (f"Berechnete Fläche weicht um {dev_txt} % von der Katasterfläche ab – "
905 "bitte Umriss kontrollieren.")
906 forest_m2 = sum(int(n["flaecheM2"] or 0) for n in land_use
907 if n["code"] in FOREST_CODES or str(n["nutzung"] or "").lower().startswith(("wald", "wäld")))
908
909 forest_share = round(100 * forest_m2 / official_m2, 1) if official_m2 else None
910 if forest_share is not None and forest_share < 90:
911 if forest_m2 == 0:
912 note = "Laut Kataster kein Wald auf diesem Grundstück."
913 else:
914 forest_ha = f"{forest_m2 / 10_000:.2f}".replace(".", ",")
915 note = (f"Gemischtes Grundstück: {forest_share:.0f} % Wald ({forest_ha} ha). "
916 "Ob ein Polygon nötig ist, hängt davon ab, welche Teilfläche als Erzeugungsfläche gilt; "
917 "das Polygon umfasst das ganze Grundstück.")
918 message = f"{message} {note}"
919 area_ha = (official_m2 or area) / 10_000
920
921 feature = {
922 "type": "Feature",
923 "geometry": geometry,
924 "properties": {
925 "ProductionPlace": f"KG {ref.kg} {kg_name or ''} Gst. {ref.gnr}".replace(" ", " ").strip(),
926 "ProducerCountry": "AT",
927 "Area": round(area_ha, 4),
928 "kg": ref.kg,
929 "gnr": ref.gnr,
930 "reference": ref.reference,
931 },
932 }
933
934 return {
935 **base,
936 "status": status,
937 "message": message,
938 "parcelId": parcel_id,
939 "kgName": kg_name,
940 "municipality": municipality,
941 "cadastreType": {"G": "Grenzkataster", "E": "Grundsteuerkataster"}.get(rstatus or "", rstatus),
942 "landUse": land_use,
943 "forestAreaM2": forest_m2,
944 "forestSharePercent": forest_share,
945 "officialAreaM2": official_m2,
946 "computedAreaM2": round(area, 1),
947 "areaDeviationPercent": None if deviation is None else round(deviation, 2),
948 "areaHa": round(area_ha, 4),
949 "eudrPolygonRequired": area_ha > EUDR_POLYGON_THRESHOLD_HA,
950 "interiorPoint": [round(ipt[0], p), round(ipt[1], p)],
951 "vertexCount": sum(len(r) for poly in out_polys for r in poly),
952 "tilesUsed": len(visited),
953 "geometry": geometry,
954 "geojsonFeature": feature,
955 }
956
957
958def feature_collection(items: list[dict]) -> dict:
959 return {
960 "type": "FeatureCollection",
961 "features": [i["geojsonFeature"] for i in items if i.get("geojsonFeature")],
962 }
963
964
965
966
967
968MAX_PARALLEL = 3
969
970
971async def main() -> None:
972 async with Actor:
973 inp = await Actor.get_input() or {}
974 raw_items = inp.get("parcels") or []
975 if isinstance(raw_items, str):
976 raw_items = [line for line in raw_items.splitlines() if line.strip()]
977 if not raw_items:
978 await Actor.fail(status_message="Keine Grundstücke angegeben (Feld 'parcels').")
979 return
980
981 resolver = ParcelResolver(
982 BevClient(),
983 max_tiles=int(inp.get("maxTilesPerParcel") or 400),
984 precision=int(inp.get("coordinatePrecision") or 7),
985 )
986 include_geojson = inp.get("includeGeoJsonFeature", True)
987
988 refs, results = [], []
989 for item in raw_items:
990 try:
991 refs.append(parse_parcel(item))
992 except InputError as e:
993 row = {"input": str(item), "status": "UNGUELTIGE_EINGABE", "message": str(e)}
994 results.append(row)
995 await Actor.push_data(row)
996
997 Actor.log.info(f"{len(refs)} Grundstücke werden verarbeitet …")
998 sem = asyncio.Semaphore(MAX_PARALLEL)
999 done = 0
1000
1001 async def work(ref):
1002 nonlocal done
1003 async with sem:
1004 try:
1005 row = await asyncio.to_thread(resolver.resolve, ref)
1006 except Exception as e:
1007 row = {"input": ref.raw, "kg": ref.kg, "gnr": ref.gnr, "status": "FEHLER", "message": str(e)}
1008 stored = dict(row)
1009 if not include_geojson:
1010 stored.pop("geojsonFeature", None)
1011 await Actor.push_data(stored)
1012 if row.get("status") in ("OK", "PRUEFEN"):
1013 try:
1014 await Actor.charge(event_name="parcel-polygon")
1015 except Exception:
1016 pass
1017 results.append(row)
1018 done += 1
1019 await Actor.set_status_message(f"{done}/{len(refs)} Grundstücke verarbeitet")
1020
1021 await asyncio.gather(*(work(r) for r in refs))
1022
1023 fc = feature_collection([r for r in results if r.get("status") in ("OK", "PRUEFEN")])
1024 await Actor.set_value(
1025 "EUDR_POLYGONE", json.dumps(fc, ensure_ascii=False), content_type="application/geo+json"
1026 )
1027 ok = sum(1 for r in results if r.get("status") == "OK")
1028 check = sum(1 for r in results if r.get("status") == "PRUEFEN")
1029 await Actor.set_status_message(
1030 f"Fertig: {ok} OK, {check} bitte prüfen, {len(results) - ok - check} ohne Ergebnis. "
1031 "Alle Polygone als GeoJSON im Key-Value-Store unter EUDR_POLYGONE."
1032 )
