aboutsummaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
-rw-r--r--webmap-download-mrr.py189
1 files changed, 98 insertions, 91 deletions
diff --git a/webmap-download-mrr.py b/webmap-download-mrr.py
index 6a25e98..e273ed2 100644
--- a/webmap-download-mrr.py
+++ b/webmap-download-mrr.py
@@ -303,6 +303,103 @@ class WMS:
logging.debug('is_uptodate(%s, %s, …) == True', layername, path)
return True
+ def getFeatures(self, layername, path_strings, gt, progress=None):
+ crs = None
+ max_probe_count = 3
+ probes = []
+ for d in path_strings:
+ d0 = d
+ m = POINT0_RE.match(d)
+ if m is None:
+ raise Exception(f"Can't parse <path>'s \"d\" attribute value: \"{d0}\"")
+ d = d.removeprefix(m.group())
+ x = float(m.group(1))
+ y = float(m.group(2))
+ path = [ [x, y] ]
+
+ # get first part of the path (until the pen is moved)
+ while True:
+ m = POINT_RE.match(d)
+ if m is None:
+ break
+ d = d.removeprefix(m.group())
+ x = float(m.group(1))
+ y = float(m.group(2))
+ if path[-1] != [x, y]:
+ path.append([x, y])
+ if d.startswith('Z') and len(path) > 0:
+ if path[0] != path[-1]:
+ path.append(path[0])
+ d = d.removeprefix("Z")
+ logging.debug('d="%s" -> path=%s, rest="%s"', d0, path, d)
+
+ # convert to an outer ring
+ ring0 = ogr.Geometry(ogr.wkbLinearRing)
+ for px, py in path:
+ x, y = gdal.ApplyGeoTransform(gt, px, py)
+ ring0.AddPoint_2D(x, y)
+
+ poly0 = ogr.Geometry(ogr.wkbPolygon)
+ poly0.AddGeometry(ring0)
+
+ # randomly probe once per 10km² (1000ha), but not more than max_probe_count
+ probe_count = poly0.GetArea() / 1e7
+ if probe_count < 1:
+ probe_count = 1
+ elif probe_count > max_probe_count:
+ probe_count = max_probe_count
+ else:
+ probe_count = round(probe_count)
+
+ probes.append( (poly0.MakeValid(), probe_count) )
+
+ if len(probes) == 0:
+ return {}, crs
+
+ features = {}
+ probe_count_sum = sum([n for _, n in probes])
+ pbar = None
+ try:
+ if progress is not None:
+ pbar = progress(total=probe_count_sum, leave=False, unit='')
+ for geom, probe_count in probes:
+ triangulation = triangulate(geom)
+
+ # always try the centroid first (later fallback to a random
+ # point if the centroid is not in geom)
+ x, y = nearCentroidWithin(geom)
+
+ # crude probing and add feature(s) found along the way
+ for i in range(probe_count):
+ if i > 0 or x is None or y is None:
+ x, y = getRandomPoint(triangulation)
+
+ data = self.getFeatureInfo(x, y, layername, res=1024)
+ crs2 = data.get('crs', None)
+ data = data['features']
+ if crs2 is not None:
+ if crs is None:
+ crs = crs2
+ elif crs != crs2 and len(data) > 0:
+ # the server should never return different CRSs
+ raise Exception(f'{json.dumps(crs)} != {json.dumps(crs2)}')
+
+ for feature in data:
+ fpr, feature = mkFeature(feature, layername)
+ features[fpr] = feature
+ data = None
+
+ if pbar is not None:
+ pbar.update(1)
+ geom = None
+ finally:
+ del probes[:]
+ if pbar is not None:
+ pbar.close()
+
+ logging.info('Found %d unique feature(s) after crude probing', len(features))
+ return features, crs
+
# Triangulate the input geometry
def triangulate(geom):
geomType = geom.GetGeometryType()
@@ -517,97 +614,7 @@ def download(dl, dest, dir_fd=None, headers={}, session=requests, progress=None)
path_strings.add(d)
logging.info('Found %d unique SVG path(s)', len(path_strings))
- crs = None
- max_probe_count = 3
- probes = []
- for d in path_strings:
- d0 = d
- m = POINT0_RE.match(d)
- if m is None:
- raise Exception(f"Can't parse <path>'s \"d\" attribute value: \"{d0}\"")
- d = d.removeprefix(m.group())
- x = float(m.group(1))
- y = float(m.group(2))
- path = [ [x, y] ]
-
- # get first part of the path (until the pen is moved)
- while True:
- m = POINT_RE.match(d)
- if m is None:
- break
- d = d.removeprefix(m.group())
- x = float(m.group(1))
- y = float(m.group(2))
- if path[-1] != [x, y]:
- path.append([x, y])
- if d.startswith('Z') and len(path) > 0:
- if path[0] != path[-1]:
- path.append(path[0])
- d = d.removeprefix("Z")
- logging.debug('d="%s" -> path=%s, rest="%s"', d0, path, d)
-
- # convert to an outer ring
- ring0 = ogr.Geometry(ogr.wkbLinearRing)
- for px, py in path:
- x, y = gdal.ApplyGeoTransform(gt, px, py)
- ring0.AddPoint_2D(x, y)
-
- poly0 = ogr.Geometry(ogr.wkbPolygon)
- poly0.AddGeometry(ring0)
-
- # randomly probe once per 10km² (1000ha), but not more than max_probe_count
- probe_count = poly0.GetArea() / 1e7
- if probe_count < 1:
- probe_count = 1
- elif probe_count > max_probe_count:
- probe_count = max_probe_count
- else:
- probe_count = round(probe_count)
-
- probes.append( (poly0.MakeValid(), probe_count) )
-
- features = {}
- probe_count_sum = sum([n for _, n in probes])
- pbar = None
- try:
- if progress is not None:
- pbar = progress(total=probe_count_sum, leave=False, unit='')
- for geom, probe_count in probes:
- triangulation = triangulate(geom)
-
- # always try the centroid first (later fallback to a random
- # point if the centroid is not in geom)
- x, y = nearCentroidWithin(geom)
-
- # crude probing and add feature(s) found along the way
- for i in range(probe_count):
- if i > 0 or x is None or y is None:
- x, y = getRandomPoint(triangulation)
-
- data = wms.getFeatureInfo(x, y, layername, res=1024)
- crs2 = data.get('crs', None)
- data = data['features']
- if crs2 is not None:
- if crs is None:
- crs = crs2
- elif crs != crs2 and len(data) > 0:
- # the server should never return different CRSs
- raise Exception(f'{json.dumps(crs)} != {json.dumps(crs2)}')
-
- for feature in data:
- fpr, feature = mkFeature(feature, layername)
- features[fpr] = feature
- data = None
-
- if pbar is not None:
- pbar.update(1)
- geom = None
- finally:
- del probes[:]
- if pbar is not None:
- pbar.close()
-
- logging.info('Found %d unique feature(s) after crude probing', len(features))
+ features, crs = wms.getFeatures(layername, path_strings, gt, progress=progress)
# refine geometries
if len(features) > 0: