Repository: geofileops/geofileops
Pull request #584 (2024-12-08 14:04:06)

ENH: support a list of input files for create_table_as_sql

codecov[bot]: ## [Codecov](https://app.codecov.io/gh/geofileops/geofileops/pull/584?dropdown=coverage&src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=geofileops) Report
Attention: Patch coverage is `98.24561%` with `1 line` in your changes missing coverage. Please review.
> Project coverage is 90.25%. Comparing base [(`5f08842`)](https://app.codecov.io/gh/geofileops/geofileops/commit/5f08842e475a5aaff2f4ace75e7bc77c9c2fe023?dropdown=coverage&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=geofileops) to head [(`f39a1cb`)](https://app.codecov.io/gh/geofileops/geofileops/commit/f39a1cb77207c7ea347df39c28d3760cdec90ed5?dropdown=coverage&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=geofileops).
> Report is 72 commits behind head on main.

| [Files with missing lines](https://app.codecov.io/gh/geofileops/geofileops/pull/584?dropdown=coverage&src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=geofileops) | Patch % | Lines |
|---|---|---|
| [geofileops/util/\_geoops\_sql.py](https://app.codecov.io/gh/geofileops/geofileops/pull/584?src=pr&el=tree&filepath=geofileops%2Futil%2F_geoops_sql.py&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=geofileops#diff-Z2VvZmlsZW9wcy91dGlsL19nZW9vcHNfc3FsLnB5) | 94.73% | [1 Missing :warning: ](https://app.codecov.io/gh/geofileops/geofileops/pull/584?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=geofileops) |

<details><summary>Additional details and impacted files</summary>


```diff
@@            Coverage Diff             @@
##             main     #584      +/-   ##
==========================================
- Coverage   93.11%   90.25%   -2.87%     
==========================================
  Files          21       21              
  Lines        4113     4268     +155     
==========================================
+ Hits         3830     3852      +22     
- Misses        283      416     +133     
```

</details>

[:umbrella: View full report in Codecov by Sentry](https://app.codecov.io/gh/geofileops/geofileops/pull/584?dropdown=coverage&src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=geofileops).   
:loudspeaker: Have feedback on the report? [Share it here](https://about.codecov.io/codecov-pr-comment-feedback/?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=geofileops).


Patch:
diff --git a/geofileops/util/_geoops_sql.py b/geofileops/util/_geoops_sql.py
index 442f182..f50d6fc 100644
--- a/geofileops/util/_geoops_sql.py
+++ b/geofileops/util/_geoops_sql.py
@@ -628,7 +628,7 @@ def _single_layer_vector_operation(
                 sql_tmp = sql_template.format(batch_filter="")
                 cols = _sqlite_util.get_columns(
                     sql_stmt=sql_tmp,
-                    input1_path=processing_params.input1_path,
+                    input_path=processing_params.input1_path,
                 )
                 geom_selected = input_layerinfo.geometrycolumn in cols
 
@@ -647,7 +647,7 @@ def _single_layer_vector_operation(
             # Get all columns of the sql_template
             sql_tmp = sql_template.format(batch_filter="")
             cols = _sqlite_util.get_columns(
-                sql_stmt=sql_tmp, input1_path=processing_params.input1_path
+                sql_stmt=sql_tmp, input_path=processing_params.input1_path
             )
             attributes = [
                 col for col in cols if col.lower() != input_layerinfo.geometrycolumn
@@ -2658,6 +2658,9 @@ def _two_layer_vector_operation(
     if tmp_dir is None:
         tmp_dir = _io_util.create_tempdir(f"geofileops/{operation_name}")
 
+    # Check if crs are the same in the input layers + use it (if there is one)
+    output_crs = _check_crs(input1_path, input1_layer, input2_path, input2_layer)
+
     # Prepare output filename
     tmp_output_path = tmp_dir / output_path.name
     tmp_output_path.parent.mkdir(exist_ok=True, parents=True)
@@ -2776,8 +2779,7 @@ def _two_layer_vector_operation(
         )
         column_datatypes = _sqlite_util.get_columns(
             sql_stmt=sql_stmt,
-            input1_path=processing_params.input1_path,
-            input2_path=processing_params.input2_path,
+            input_path=[processing_params.input1_path, processing_params.input2_path],
         )
 
         # Apply gridsize if it is specified
@@ -2894,16 +2896,15 @@ def _two_layer_vector_operation(
 
                 # Remark: this temp file doesn't need spatial index
                 future = calculate_pool.submit(
-                    calculate_two_layers,
+                    _calculate_two_layers,
                     input1_path=processing_params.batches[batch_id]["input1_path"],
-                    input1_layer=processing_params.batches[batch_id]["input1_layer"],
                     input2_path=processing_params.batches[batch_id]["input2_path"],
-                    input2_layer=processing_params.batches[batch_id]["input2_layer"],
                     output_path=tmp_partial_output_path,
                     sql_stmt=sql_stmt,
                     output_layer=output_layer,
                     explodecollections=explodecollections_now,
                     force_output_geometrytype=output_geometrytype_now,
+                    output_crs=output_crs,
                     use_ogr=use_ogr,
                     create_spatial_index=False,
                     column_datatypes=column_datatypes,
@@ -3008,16 +3009,48 @@ def _two_layer_vector_operation(
             shutil.rmtree(tmp_dir, ignore_errors=True)
 
 
-def calculate_two_layers(
+def _check_crs(
+    input1_path: Path, input1_layer: str, input2_path: Path, input2_layer: str
+) -> int:
+    input1_info = gfo.get_layerinfo(input1_path, input1_layer, raise_on_nogeom=False)
+    input2_info = None
+    if input2_path is not None:
+        input2_info = gfo.get_layerinfo(
+            input2_path, input2_layer, raise_on_nogeom=False
+        )
+    crs_epsg = -1
+    if input1_info.crs is not None:
+        crs_epsg1 = input1_info.crs.to_epsg()
+        if crs_epsg1 is not None:
+            crs_epsg = crs_epsg1
+        # If input 2 also has a crs, check if it is the same.
+        if (
+            input2_info is not None
+            and input2_info.crs is not None
+            and crs_epsg1 != input2_info.crs.to_epsg()
+        ):
+            warnings.warn(
+                "input1 layer doesn't have the same crs as input2 layer: "
+                f"{input1_info.crs} vs {input2_info.crs}",
+                stacklevel=5,
+            )
+    elif input2_info is not None and input2_info.crs is not None:
+        crs_epsg2 = input2_info.crs.to_epsg()
+        if crs_epsg2 is not None:
+            crs_epsg = crs_epsg2
+
+    return crs_epsg
+
+
+def _calculate_two_layers(
     input1_path: Path,
-    input1_layer: str,
     input2_path: Path,
-    input2_layer: str,
     output_path: Path,
     sql_stmt: str,
     output_layer: str,
     explodecollections: bool,
     force_output_geometrytype: Optional[GeometryType],
+    output_crs: int,
     create_spatial_index: bool,
     column_datatypes: dict,
     use_ogr: bool,
@@ -3030,14 +3063,12 @@ def calculate_two_layers(
             output_name = f"{output_path.stem}_tmp{output_path.suffix}"
             output_tmp_path = output_path.parent / output_name
         _sqlite_util.create_table_as_sql(
-            input1_path=input1_path,
-            input1_layer=input1_layer,
-            input2_path=input2_path,
-            input2_layer=input2_layer,
+            input_path=[input1_path, input2_path],
             output_path=output_tmp_path,
             sql_stmt=sql_stmt,
             output_layer=output_layer,
             output_geometrytype=force_output_geometrytype,
+            output_crs=output_crs,
             create_spatial_index=create_spatial_index,
             profile=_sqlite_util.SqliteProfile.SPEED,
             column_datatypes=column_datatypes,
diff --git a/geofileops/util/_sqlite_util.py b/geofileops/util/_sqlite_util.py
index f508dc7..6763998 100644
--- a/geofileops/util/_sqlite_util.py
+++ b/geofileops/util/_sqlite_util.py
@@ -7,6 +7,7 @@ import pprint
 import shutil
 import sqlite3
 import tempfile
+from collections.abc import Iterable
 from pathlib import Path
 from typing import Optional, Union
 
@@ -125,16 +126,19 @@ def create_new_spatialdb(path: Path, crs_epsg: Optional[int] = None):
 
 def get_columns(
     sql_stmt: str,
-    input1_path: Path,
-    input2_path: Optional[Path] = None,
+    input_path: Union[Path, list[Path]],
     empty_output_ok: bool = True,
     use_spatialite: bool = True,
     output_geometrytype: Optional[GeometryType] = None,
 ) -> dict[str, str]:
     # Create temp output db to be sure the output DB is writable, even though we only
     # create a temporary table.
+    if isinstance(input_path, (str, Path)):
+        input_path = [input_path]
+    assert isinstance(input_path, Iterable)
+
     tmp_dir = Path(tempfile.mkdtemp(prefix="geofileops/get_columns_"))
-    tmp_path = tmp_dir / f"temp{input1_path.suffix}"
+    tmp_path = tmp_dir / f"temp{input_path[0].suffix}"
     create_new_spatialdb(path=tmp_path)
 
     sql = None
@@ -147,27 +151,27 @@ def get_columns(
                 sql = "SELECT EnableGpkgMode();"
                 conn.execute(sql)
 
-        # Attach to input1
-        input1_databasename = "input1"
-        sql = f"ATTACH DATABASE ? AS {input1_databasename}"
-        dbSpec = (str(input1_path),)
-        conn.execute(sql, dbSpec)
-
-        # If input2 isn't the same database input1, attach to it
-        input2_databasename = None
-        if input2_path is not None:
-            if input2_path == input1_path:
-                input2_databasename = input1_databasename
-            else:
-                input2_databasename = "input2"
-                sql = f"ATTACH DATABASE ? AS {input2_databasename}"
-                dbSpec = (str(input2_path),)
+            # Attach to all input databases
+            db_path_to_placeholder: dict[Path, str] = {}
+            db_placeholder_to_value: dict[str, str] = {}
+            for index, path in enumerate(input_path):
+                db_placeholder = f"input{index + 1}_databasename"
+                if path in db_path_to_placeholder:
+                    db_placeholder_to_value[db_placeholder] = db_placeholder_to_value[
+                        db_path_to_placeholder[path]
+                    ]
+                    continue
+
+                input_databasename = f"input{index + 1}"
+                sql = f"ATTACH DATABASE ? AS {input_databasename}"
+                dbSpec = (str(path),)
                 conn.execute(sql, dbSpec)
+                db_path_to_placeholder[path] = db_placeholder
+                db_placeholder_to_value[db_placeholder] = input_databasename
 
         # Prepare sql statement for execute
         sql_stmt_prepared = sql_stmt.format(
-            input1_databasename=input1_databasename,
-            input2_databasename=input2_databasename,
+            **db_placeholder_to_value,
             batch_filter="",
         )
 
@@ -271,14 +275,12 @@ def get_columns(
 
 
 def create_table_as_sql(
-    input1_path: Path,
-    input1_layer: str,
-    input2_path: Path,
-    input2_layer: str,
+    input_path: Union[Path, list[Path]],
     output_path: Path,
     sql_stmt: str,
     output_layer: str,
     output_geometrytype: Optional[GeometryType],
+    output_crs: int,
     append: bool = False,
     update: bool = False,
     create_spatial_index: bool = True,
@@ -289,14 +291,12 @@ def create_table_as_sql(
     """Execute sql statement and save the result in the output file.
 
     Args:
-        input1_path (Path): the path to the 1st input file.
-        input1_layer (str): the layer/table to select from in het 1st input file
-        input2_path (Path): the path to the 2nd input file.
-        input2_layer (str): the layer/table to select from in het 2nd input file
+        input_path (Path or list): the path(s) to the input file(s).
         output_path (Path): the path where the output file needs to be created/appended.
         sql_stmt (str): SELECT statement to run on the input files.
         output_layer (str): layer/table name to use.
         output_geometrytype (Optional[GeometryType]): geometry type of the output.
+        output_crs (int): epsg code of crs of the output file.
         append (bool, optional): True to append to an existing file. Defaults to False.
         update (bool, optional): True to append to an existing layer. Defaults to False.
         create_spatial_index (bool, optional): True to create a spatial index on the
@@ -321,42 +321,22 @@ def create_table_as_sql(
     # Check input parameters
     if append is True or update is True:
         raise ValueError("append=True nor update=True are implemented.")
+
+    # All input files and the output file must have the same suffix.
+    if isinstance(input_path, (str, Path)):
+        input_path = [input_path]
+    assert isinstance(input_path, Iterable)
+
     output_suffix_lower = output_path.suffix.lower()
-    if output_suffix_lower != input1_path.suffix.lower():
-        raise ValueError("output_path and both input paths must have the same suffix!")
-    if input2_path is not None and output_suffix_lower != input2_path.suffix.lower():
-        raise ValueError("output_path and both input paths must have the same suffix!")
-
-    # Check if crs are the same in the input layers + use it (if there is one)
-    input1_info = gfo.get_layerinfo(input1_path, input1_layer, raise_on_nogeom=False)
-    input2_info = None
-    if input2_path is not None:
-        input2_info = gfo.get_layerinfo(
-            input2_path, input2_layer, raise_on_nogeom=False
-        )
-    crs_epsg = -1
-    if input1_info.crs is not None:
-        crs_epsg1 = input1_info.crs.to_epsg()
-        if crs_epsg1 is not None:
-            crs_epsg = crs_epsg1
-        # If input 2 also has a crs, check if it is the same.
-        if (
-            input2_info is not None
-            and input2_info.crs is not None
-            and crs_epsg1 != input2_info.crs.to_epsg()
-        ):
-            logger.warning(
-                "input1 layer doesn't have the same crs as input2 layer: "
-                f"{input1_info.crs} vs {input2_info.crs}"
+    for path in input_path:
+        if output_suffix_lower != path.suffix.lower():
+            raise ValueError(
+                "output_path and all input paths must have the same suffix!"
             )
-    elif input2_info is not None and input2_info.crs is not None:
-        crs_epsg2 = input2_info.crs.to_epsg()
-        if crs_epsg2 is not None:
-            crs_epsg = crs_epsg2
 
     # If output file doesn't exist yet, create and init it
     if not output_path.exists():
-        create_new_spatialdb(path=output_path, crs_epsg=crs_epsg)
+        create_new_spatialdb(path=output_path, crs_epsg=output_crs)
 
     sql = None
     conn = sqlite3.connect(output_path, detect_types=sqlite3.PARSE_DECLTYPES, uri=True)
@@ -390,26 +370,23 @@ def create_table_as_sql(
             sql = f"PRAGMA soft_heap_limit={1024*1024*1024};"
             conn.execute(sql)
 
-            # If attach to input1
-            input1_databasename = "input1"
-            sql = f"ATTACH DATABASE ? AS {input1_databasename}"
-            dbSpec = (str(input1_path),)
-            # input1_uri = f"file:{input1_path}?immutable=1"
-            # dbSpec = (str(input1_uri),)
-            conn.execute(sql, dbSpec)
-
-            # If input2 isn't the same database input1, attach to it
-            input2_databasename = None
-            if input2_path is not None:
-                if input2_path == input1_path:
-                    input2_databasename = input1_databasename
-                else:
-                    input2_databasename = "input2"
-                    sql = f"ATTACH DATABASE ? AS {input2_databasename}"
-                    dbSpec = (str(input2_path),)
-                    # input2_uri = f"file:{input1_path}?immutable=1"
-                    # dbSpec = (str(input2_uri),)
-                    conn.execute(sql, dbSpec)
+            # Attach to all input databases
+            db_path_to_placeholder: dict[Path, str] = {}
+            db_placeholder_to_value: dict[str, str] = {}
+            for index, path in enumerate(input_path):
+                db_placeholder = f"input{index+1}_databasename"
+                if path in db_path_to_placeholder:
+                    db_placeholder_to_value[db_placeholder] = db_placeholder_to_value[
+                        db_path_to_placeholder[path]
+                    ]
+                    continue
+
+                input_databasename = f"input{index + 1}"
+                sql = f"ATTACH DATABASE ? AS {input_databasename}"
+                dbSpec = (str(path),)
+                conn.execute(sql, dbSpec)
+                db_path_to_placeholder[path] = db_placeholder
+                db_placeholder_to_value[db_placeholder] = input_databasename
 
             # Use the sqlite profile specified
             if profile is SqliteProfile.SPEED:
@@ -421,8 +398,7 @@ def create_table_as_sql(
                 # it doesn't hurt and maybe on other platforms...
                 for databasename in [
                     output_databasename,
-                    input1_databasename,
-                    input2_databasename,
+                    *db_placeholder_to_value.values(),
                 ]:
                     conn.execute(f"PRAGMA {databasename}.journal_mode=OFF;")
 
@@ -435,8 +411,7 @@ def create_table_as_sql(
             if column_types is None:
                 column_types = get_columns(
                     sql_stmt=sql_stmt,
-                    input1_path=input1_path,
-                    input2_path=input2_path,
+                    input_path=input_path,
                     empty_output_ok=empty_output_ok,
                     use_spatialite=True,
                     output_geometrytype=output_geometrytype,
@@ -447,10 +422,7 @@ def create_table_as_sql(
                 output_geometrytype = GeometryType(column_types["geom"])
 
             # Prepare sql statement
-            sql_stmt = sql_stmt.format(
-                input1_databasename=input1_databasename,
-                input2_databasename=input2_databasename,
-            )
+            sql_stmt = sql_stmt.format(**db_placeholder_to_value)
 
             # Now we can create the table
             # Create output table using the gpkgAddGeometryColumn() function
@@ -500,7 +472,7 @@ def create_table_as_sql(
                         table_name, data_type, identifier, description, last_change,
                         min_x, min_y, max_x, max_y, srs_id)
                     VALUES ('{output_layer}', '{data_type}', NULL, '', DATETIME(),
-                        NULL, NULL, NULL, NULL, {to_string_for_sql(crs_epsg)});
+                        NULL, NULL, NULL, NULL, {to_string_for_sql(output_crs)});
                 """
                 conn.execute(sql)
 
@@ -511,7 +483,7 @@ def create_table_as_sql(
                         INSERT INTO {output_databasename}.gpkg_geometry_columns (
                             table_name, column_name, geometry_type_name, srs_id, z, m)
                         VALUES ('{output_layer}', 'geom', '{output_geometrytype.name}',
-                            {to_string_for_sql(crs_epsg)}, 0, 0);
+                            {to_string_for_sql(output_crs)}, 0, 0);
                     """
                     conn.execute(sql)
 
@@ -526,8 +498,8 @@ def create_table_as_sql(
                     sql = f"""
                         SELECT RecoverGeometryColumn(
                           '{output_layer}', 'geom',
-                          {to_string_for_sql(crs_epsg)}, '{output_geometrytype.name}');
-                    """
+                          {to_string_for_sql(output_crs)}, '{output_geometrytype.name}');
+                    """  # noqa: E501
                     conn.execute(sql)
 
             # Insert data using the sql statement specified


Test patch:
diff --git a/tests/test_geofileops_twolayers.py b/tests/test_geofileops_twolayers.py
index cee2052..7b73088 100644
--- a/tests/test_geofileops_twolayers.py
+++ b/tests/test_geofileops_twolayers.py
@@ -1042,6 +1042,20 @@ def test_intersection_columns_fid(
             assert sorted(output_gdf.l2_FiD.unique().tolist()) == [1, 2, 3, 4, 5]
 
 
+def test_intersection_different_crs(tmp_path):
+    # Prepare test data
+    input1_path = test_helper.get_testfile("polygon-parcel", suffix=".gpkg")
+    input2_path = test_helper.get_testfile("polygon-parcel", suffix=".gpkg", epsg=4326)
+    output_path = tmp_path / "output.gpkg"
+
+    with pytest.warns(match="input1 layer doesn't have the same crs as input2 layer"):
+        gfo.intersection(
+            input1_path=input1_path,
+            input2_path=input2_path,
+            output_path=output_path,
+        )
+
+
 @pytest.mark.parametrize(
     "suffix, explodecollections, where_post, exp_featurecount",
     [
diff --git a/tests/test_sqlite_util.py b/tests/test_sqlite_util.py
index 7cf6a8c..7e31434 100644
--- a/tests/test_sqlite_util.py
+++ b/tests/test_sqlite_util.py
@@ -44,13 +44,11 @@ def test_create_table_as_sql(tmp_path, create_spatial_index):
             """
 
     sqlite_util.create_table_as_sql(
-        input1_path=input1_path,
-        input1_layer="parcels",
-        input2_path=input2_path,
-        input2_layer="zones",
+        input_path=[input1_path, input2_path],
         output_path=output_path,
         output_layer=output_path.stem,
         output_geometrytype=gfo.GeometryType.MULTIPOLYGON,
+        output_crs=31370,
         sql_stmt=sql_stmt,
         profile=sqlite_util.SqliteProfile.SPEED,
         create_spatial_index=create_spatial_index,
@@ -74,29 +72,26 @@ def test_create_table_as_sql(tmp_path, create_spatial_index):
         ({"append": True}, "append=True nor update=True are implemented."),
         ({"update": True}, "append=True nor update=True are implemented."),
         (
-            {"input1_path": Path("input1.sqlite")},
-            "output_path and both input paths must have the same suffix!",
+            {"input_path": [Path("input1.sqlite"), Path("input2.gpkg")]},
+            "output_path and all input paths must have the same suffix!",
         ),
         (
-            {"input2_path": Path("input2.sqlite")},
-            "output_path and both input paths must have the same suffix!",
+            {"input_path": [Path("input1.gpkg"), Path("input2.sqlite")]},
+            "output_path and all input paths must have the same suffix!",
         ),
         (
-            {"output_path": Path("output.sqlite")},
-            "output_path and both input paths must have the same suffix!",
+            {
+                "input_path": [Path("input1.gpkg"), Path("input2.gpkg")],
+                "output_path": Path("output.sqlite"),
+            },
+            "output_path and all input paths must have the same suffix!",
         ),
     ],
 )
 def test_create_table_as_sql_invalidparams(kwargs, expected_error):
     # Set default values for kwargs that are not specified:
-    if "input1_path" not in kwargs:
-        kwargs["input1_path"] = Path("input1.gpkg")
-    if "input1_layer" not in kwargs:
-        kwargs["input1_layer"] = "input1_layer"
-    if "input2_path" not in kwargs:
-        kwargs["input2_path"] = Path("input2.gpkg")
-    if "input2_layer" not in kwargs:
-        kwargs["input2_layer"] = "input2_layer"
+    if "input_path" not in kwargs:
+        kwargs["input_path"] = [Path("input1.gpkg"), Path("input2.gpkg")]
     if "output_path" not in kwargs:
         kwargs["output_path"] = Path("output.gpkg")
     if "output_layer" not in kwargs:
@@ -105,6 +100,8 @@ def test_create_table_as_sql_invalidparams(kwargs, expected_error):
         kwargs["sql_stmt"] = "SELECTje"
     if "output_geometrytype" not in kwargs:
         kwargs["output_geometrytype"] = None
+    if "output_crs" not in kwargs:
+        kwargs["output_crs"] = -1
 
     with pytest.raises(ValueError, match=expected_error):
         sqlite_util.create_table_as_sql(**kwargs)
@@ -153,7 +150,7 @@ def test_get_columns():
 
     # Run test
     columns = sqlite_util.get_columns(
-        sql_stmt=sql_stmt, input1_path=input1_path, input2_path=input2_path
+        sql_stmt=sql_stmt, input_path=[input1_path, input2_path]
     )
 
     assert len(columns) == 5
@@ -177,13 +174,11 @@ def test_create_table_as_sql_single_input(tmp_path):
           FROM "parcels" layer
     """
     sqlite_util.create_table_as_sql(
-        input1_path=input_path,
-        input1_layer="parcels",
-        input2_path=None,
-        input2_layer=None,
+        input_path=input_path,
         output_path=output_path,
         output_layer=None,
         output_geometrytype=None,
+        output_crs=31370,
         sql_stmt=sql_stmt,
     )