diff options
| -rw-r--r-- | webmap-download-mrr.py | 189 |
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: |