merge: improv GDAL import by cliping the source

Author: frodrigo

analysers/Analyser_Merge.py

@@ -36,6 +36,7 @@
 import re
 import fnmatch
 import shutil
+import subprocess
 from typing import Optional, Dict, Union, Callable
 from collections import defaultdict
 from .Analyser_Osmosis import Analyser_Osmosis
@@ -731,15 +732,25 @@ def import_(self, table, osmosis):
                 if info:
                     self.zip = info.filename
 
+            source_layer = [
+                ('/vsizip/' if self.zip else '' ) + tmp_file.name + (('/' + self.zip) if self.zip else ''),
+            ]
+            if self.layer:
+                source_layer.append(f"'{self.layer}'")
+
+            s_src = re.search('EPSG:([0-9]+)', subprocess.run(["gdalsrsinfo", "-e", *source_layer], stdout=subprocess.PIPE).stdout.decode('utf-8')).group(1)
+            wkt = PointInPolygon.PointInPolygon(self.polygon_id).polygon.as_wkt(s_src) if self.polygon_id else None
+
             select = "-select '{}'".format(','.join(self.fields)) if self.fields else ''
-            gdal = "ogr2ogr -f PostgreSQL 'PG:{}' -lco SCHEMA={} -nln '{}' -lco OVERWRITE=yes -lco GEOMETRY_NAME=geom -lco OVERWRITE=YES -lco LAUNDER=NO -skipfailures {} -t_srs EPSG:{} '{}' {}".format(
+            gdal = "ogr2ogr -f PostgreSQL 'PG:{}' -lco SCHEMA={} -nln '{}' {} -lco OVERWRITE=yes -lco GEOMETRY_NAME=geom -lco OVERWRITE=YES -lco LAUNDER=NO -skipfailures {} -t_srs EPSG:{} '{}' {}".format(
                 osmosis.config.osmosis_manager.db_string,
                 osmosis.config.osmosis_manager.db_user,
                 table,
+                f"-clipsrc '{wkt}'" if wkt else '',
                 select,
                 self.proj,
-                ('/vsizip/' if self.zip else '' ) + tmp_file.name + (('/' + self.zip) if self.zip else ''),
-                f"'{self.layer}'" if self.layer else '',
+                source_layer[0],
+                source_layer[1] if len(source_layer) >= 2 else '',
             )
             print(gdal)
             if os.system(gdal):
@@ -1168,6 +1179,7 @@ def init(self, url, name, parser, load = Load(), conflate = Conflate()):
             self.conflate.select.tags = [self.conflate.select.tags]
         self.conflate.mapping.eval_static(self)
 
+        self.parser.polygon_id = self.config.polygon_id
         self.load.osmosis = self
         self.load.polygon_id = self.config.polygon_id
         if "proj" in self.config.options:

modules/Polygon.py

@@ -22,6 +22,8 @@
 from shapely.wkt import loads
 from shapely.geometry import MultiPolygon
 from modules import downloader
+import pyproj
+from shapely.ops import transform
 
 
 class Polygon:
@@ -34,12 +36,19 @@ def __init__(self, polygon_id, cache_delay=60):
         polygon_url = u"http://polygons.openstreetmap.fr/"
         for id in polygon_id:
             url = polygon_url + "index.py?id="+str(id)
-            s = downloader.urlread(url, cache_delay)
+            downloader.urlread(url, cache_delay)
         url = polygon_url + "get_wkt.py?params=0&id=" + ",".join(map(str, polygon_id))
-        s = downloader.urlread(url, cache_delay)
-        if s.startswith("SRID="):
-            s = s.split(";", 1)[1]
-        self.polygon = loads(s)
+        self.wkt = wkt = downloader.urlread(url, cache_delay)
+        if wkt.startswith("SRID="):
+            wkt = wkt.split(";", 1)[1]
+        self.polygon = loads(wkt)
+
+    def as_wkt(self, srid) -> str:
+        wgs84 = pyproj.CRS('EPSG:4326')
+        t_src = pyproj.CRS(f'EPSG:{srid}')
+        project = pyproj.Transformer.from_crs(wgs84, t_src, always_xy=True).transform
+        t_poly = transform(project, self.polygon)
+        return t_poly
 
     def bboxes(self):
         bbox = self.polygon.bounds