Specific process for geometries with Z or M coordinate with SpatiaLit…

…e dialect (#506)

* Temporary fix for geometries with Z in SpatiaLite.

Though most of the dialects support  geom_type ending with Z in SQL queries, it seems not working with SpatiaLite.
Hence replacing Z with a blank. Added a test to cover this behaviour.

* Add a more complete test and propose an extended fix

* Handle M coordinate

* Change postgis image in test_and_publish.yml to postgis/postgis:16-3.4

Add some formatting for postgresql dialect, MULTIPOINT geom_type.

* Enable postgis_raster extension for postgresql gis db.

* Fix docs, tests and docker helpers

* Fix postgis_raster extension in CI

* Update .github/workflows/test_and_publish.yml

Co-authored-by: Adrien Berchet <adrien.berchet@gmail.com>

* Update .github/workflows/test_and_publish.yml

Co-authored-by: Adrien Berchet <adrien.berchet@gmail.com>

* Update .github/workflows/test_and_publish.yml

Co-authored-by: Adrien Berchet <adrien.berchet@gmail.com>

---------

Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Co-authored-by: Adrien Berchet <adrien.berchet@gmail.com>

.github/workflows/test_and_publish.yml

@@ -41,7 +41,7 @@ jobs:
 
     services:
       postgres:
-        image: mdillon/postgis:11
+        image: postgis/postgis:16-3.4
         env:
           POSTGRES_DB: gis
           POSTGRES_PASSWORD: gis
@@ -84,11 +84,17 @@ jobs:
         psql -h localhost -p 5432 -U gis -d gis -c 'CREATE SCHEMA gis;'
 
         # Add PostGIS extension to "gis" database
-        psql -h localhost -p 5432 -U gis -d gis -c 'CREATE EXTENSION IF NOT EXISTS postgis;'
+        psql -h localhost -p 5432 -U gis -d gis -c 'CREATE EXTENSION IF NOT EXISTS postgis SCHEMA public;'
 
         # Drop PostGIS Tiger Geocoder extension to "gis" database
         psql -h localhost -p 5432 -U gis -d gis -c 'DROP EXTENSION IF EXISTS postgis_tiger_geocoder CASCADE;'
 
+        # Add PostGISRaster extension to "gis" database
+        psql -h localhost -p 5432 -U gis -d gis -c 'CREATE EXTENSION IF NOT EXISTS postgis_raster SCHEMA public;'
+
+        # Drop PostGIS Topology extension to "gis" database
+        psql -h localhost -p 5432 -U gis -d gis -c 'DROP EXTENSION IF EXISTS postgis_topology;'
+
     # Setup MySQL
     - name: Set up MySQL
       run: |

doc/conf.py

@@ -60,6 +60,9 @@
 #
 version = release = geoalchemy2.__version__
 
+# Remove some Sphinx warnings
+suppress_warnings = ["config.cache"]
+
 # The language for content autogenerated by Sphinx. Refer to documentation
 # for a list of supported languages.
 # language = None

doc/spatialite_tutorial.rst

@@ -54,7 +54,7 @@ Declare a Mapping
 Now that we have a working connection we can go ahead and create a mapping between
 a Python class and a database table::
 
-    >>> from sqlalchemy.ext.declarative import declarative_base
+    >>> from sqlalchemy.orm import declarative_base
     >>> from sqlalchemy import Column, Integer, String
     >>> from geoalchemy2 import Geometry
     >>>

geoalchemy2/elements.py

@@ -117,6 +117,7 @@ class WKTElement(_SpatialElement):
     """
 
     _REMOVE_SRID = re.compile("(SRID=([0-9]+); ?)?(.*)")
+    SPLIT_WKT_PATTERN = re.compile(r"((SRID=\d+) *; *)?([\w ]+) *(\([\d ,\(\)]+\))")
 
     geom_from: str = "ST_GeomFromText"
     geom_from_extended_version: str = "ST_GeomFromEWKT"

geoalchemy2/types/dialects/sqlite.py

@@ -1,22 +1,50 @@
 """This module defines specific functions for SQLite dialect."""
 
+import re
+
 from geoalchemy2.elements import RasterElement
 from geoalchemy2.elements import WKBElement
 from geoalchemy2.elements import WKTElement
 from geoalchemy2.shape import to_shape
 
 
+def format_geom_type(wkt, default_srid=None):
+    """Format the Geometry type for SQLite."""
+    match = re.match(WKTElement.SPLIT_WKT_PATTERN, wkt)
+    if match is None:
+        return wkt
+    _, srid, geom_type, coords = match.groups()
+    geom_type = geom_type.replace(" ", "")
+    if geom_type.endswith("ZM"):
+        geom_type = geom_type[:-2]
+    elif geom_type.endswith("Z"):
+        geom_type = geom_type[:-1]
+    if srid is None and default_srid is not None:
+        srid = f"SRID={default_srid}"
+    if srid is not None:
+        return "%s;%s%s" % (srid, geom_type, coords)
+    else:
+        return "%s%s" % (geom_type, coords)
+
+
 def bind_processor_process(spatial_type, bindvalue):
     if isinstance(bindvalue, WKTElement):
-        if bindvalue.extended:
-            return "%s" % (bindvalue.data)
-        else:
-            return "SRID=%d;%s" % (bindvalue.srid, bindvalue.data)
+        return format_geom_type(
+            bindvalue.data,
+            default_srid=bindvalue.srid if bindvalue.srid >= 0 else spatial_type.srid,
+        )
     elif isinstance(bindvalue, WKBElement):
         # With SpatiaLite we use Shapely to convert the WKBElement to an EWKT string
         shape = to_shape(bindvalue)
-        return "SRID=%d;%s" % (bindvalue.srid, shape.wkt)
+        # shapely.wkb.loads returns geom_type with a 'Z', for example, 'LINESTRING Z'
+        # which is a limitation with SpatiaLite. Hence, a temporary fix.
+        res = format_geom_type(
+            shape.wkt, default_srid=bindvalue.srid if bindvalue.srid >= 0 else spatial_type.srid
+        )
+        return res
     elif isinstance(bindvalue, RasterElement):
         return "%s" % (bindvalue.data)
+    elif isinstance(bindvalue, str):
+        return format_geom_type(bindvalue, default_srid=spatial_type.srid)
     else:
         return bindvalue

test_container/Dockerfile

@@ -5,7 +5,7 @@ RUN /install_requirements.sh
 
 COPY ./helpers/init_postgres.sh /
 env PGDATA="/var/lib/postgresql/data"
-env POSTGRES_PATH="/usr/lib/postgresql/14"
+env POSTGRES_PATH="/usr/lib/postgresql/16"
 RUN su postgres -c /init_postgres.sh
 
 ENV SPATIALITE_LIBRARY_PATH="/usr/lib/x86_64-linux-gnu/mod_spatialite.so"

test_container/helpers/install_requirements.sh

@@ -29,9 +29,9 @@ packages=(
     python3.12-venv
 
     # PostgreSQL and PostGIS
-    postgresql
-    postgresql-14-postgis-3
-    postgresql-14-postgis-3-scripts
+    postgresql-16
+    postgresql-16-postgis-3
+    postgresql-16-postgis-3-scripts
     libpq-dev
     libgeos-dev
 
@@ -56,8 +56,10 @@ packages=(
 export DEBIAN_FRONTEND=noninteractive
 
 apt-get update -y
-apt-get install --no-install-recommends -y software-properties-common gnupg2
+apt-get install --no-install-recommends -y software-properties-common gnupg2 wget
 add-apt-repository ppa:deadsnakes/ppa
+sh -c 'echo "deb https://apt.PostgreSQL.org/pub/repos/apt $(lsb_release -cs)-pgdg main" > /etc/apt/sources.list.d/pgdg.list'
+wget --quiet -O - https://www.PostgreSQL.org/media/keys/ACCC4CF8.asc | apt-key add -
 apt-get update -y
 
 apt-get install --no-install-recommends -y "${packages[@]}"

tests/test_functional.py

@@ -275,6 +275,140 @@ def test_insert(self, conn, Lake, setup_tables):
         srid = conn.execute(row[1].ST_SRID()).scalar()
         assert srid == 4326
 
+    @pytest.mark.parametrize(
+        "geom_type,wkt",
+        [
+            pytest.param("POINT", "(1 2)", id="Point"),
+            pytest.param("POINTZ", "(1 2 3)", id="Point Z"),
+            pytest.param("POINTM", "(1 2 3)", id="Point M"),
+            pytest.param("POINTZM", "(1 2 3 4)", id="Point ZM"),
+            pytest.param("LINESTRING", "(1 2, 3 4)", id="LineString"),
+            pytest.param("LINESTRINGZ", "(1 2 3, 4 5 6)", id="LineString Z"),
+            pytest.param("LINESTRINGM", "(1 2 3, 4 5 6)", id="LineString M"),
+            pytest.param("LINESTRINGZM", "(1 2 3 4, 5 6 7 8)", id="LineString ZM"),
+            pytest.param("POLYGON", "((1 2, 3 4, 5 6, 1 2))", id="Polygon"),
+            pytest.param("POLYGONZ", "((1 2 3, 4 5 6, 7 8 9, 1 2 3))", id="Polygon Z"),
+            pytest.param("POLYGONM", "((1 2 3, 4 5 6, 7 8 9, 1 2 3))", id="Polygon M"),
+            pytest.param(
+                "POLYGONZM", "((1 2 3 4,  5 6 7 8, 9 10 11 12, 1 2 3 4))", id="Polygon ZM"
+            ),
+            pytest.param("MULTIPOINT", "(1 2, 3 4)", id="Multi Point"),
+            pytest.param("MULTIPOINTZ", "(1 2 3, 4 5 6)", id="Multi Point Z"),
+            pytest.param("MULTIPOINTM", "(1 2 3, 4 5 6)", id="Multi Point M"),
+            pytest.param("MULTIPOINTZM", "(1 2 3 4, 5 6 7 8)", id="Multi Point ZM"),
+            pytest.param("MULTILINESTRING", "((1 2, 3 4), (10 20, 30 40))", id="Multi LineString"),
+            pytest.param(
+                "MULTILINESTRINGZ",
+                "((1 2 3, 4 5 6), (10 20 30, 40 50 60))",
+                id="Multi LineString Z",
+            ),
+            pytest.param(
+                "MULTILINESTRINGM",
+                "((1 2 3, 4 5 6), (10 20 30, 40 50 60))",
+                id="Multi LineString M",
+            ),
+            pytest.param(
+                "MULTILINESTRINGZM",
+                "((1 2 3 4, 5 6 7 8), (10 20 30 40, 50 60 70 80))",
+                id="Multi LineString ZM",
+            ),
+            pytest.param(
+                "MULTIPOLYGON",
+                "(((1 2, 3 4, 5 6, 1 2)), ((10 20, 30 40, 50 60, 10 20)))",
+                id="Multi Polygon",
+            ),
+            pytest.param(
+                "MULTIPOLYGONZ",
+                "(((1 2 3, 4 5 6, 7 8 9, 1 2 3)), ((10 20 30, 40 50 60, 70 80 90, 10 20 30)))",
+                id="Multi Polygon Z",
+            ),
+            pytest.param(
+                "MULTIPOLYGONM",
+                "(((1 2 3, 4 5 6, 7 8 9, 1 2 3)), ((10 20 30, 40 50 60, 70 80 90, 10 20 30)))",
+                id="Multi Polygon M",
+            ),
+            pytest.param(
+                "MULTIPOLYGONZM",
+                "(((1 2 3 4, 5 6 7 8, 9 10 11 12, 1 2 3 4)),"
+                " ((10 20 30 40, 50 60 70 80, 90 100 100 120, 10 20 30 40)))",
+                id="Multi Polygon ZM",
+            ),
+        ],
+    )
+    def test_insert_all_geom_types(self, dialect_name, base, conn, metadata, geom_type, wkt):
+        """Test insertion and selection of all geometry types."""
+        ndims = 2
+        if "Z" in geom_type[-2:]:
+            ndims += 1
+        if geom_type.endswith("M"):
+            ndims += 1
+            has_m = True
+        else:
+            has_m = False
+
+        if ndims > 2 and dialect_name == "mysql":
+            # Explicitly skip MySQL dialect to show that it can only work with 2D geometries
+            pytest.xfail(reason="MySQL only supports 2D geometry types")
+
+        class GeomTypeTable(base):
+            __tablename__ = "test_geom_types"
+            id = Column(Integer, primary_key=True)
+            geom = Column(Geometry(srid=4326, geometry_type=geom_type, dimension=ndims))
+
+        metadata.drop_all(bind=conn, checkfirst=True)
+        metadata.create_all(bind=conn)
+
+        inserted_wkt = f"{geom_type}{wkt}"
+
+        # Use the DB to generate the corresponding raw WKB
+        raw_wkb = conn.execute(
+            text("SELECT ST_AsBinary(ST_GeomFromText('{}', 4326))".format(inserted_wkt))
+        ).scalar()
+
+        wkb_elem = WKBElement(raw_wkb, srid=4326)
+        inserted_elements = [
+            {"geom": inserted_wkt},
+            {"geom": f"SRID=4326;{inserted_wkt}"},
+            {"geom": WKTElement(inserted_wkt, srid=4326)},
+            {"geom": WKTElement(f"SRID=4326;{inserted_wkt}")},
+        ]
+        if dialect_name not in ["postgresql", "sqlite"] or not has_m:
+            # Currently Shapely does not support geometry types with M dimension
+            inserted_elements.append({"geom": wkb_elem})
+            inserted_elements.append({"geom": wkb_elem.as_ewkb()})
+
+        # Insert the elements
+        conn.execute(
+            GeomTypeTable.__table__.insert(),
+            inserted_elements,
+        )
+
+        # Select the elements
+        query = select(
+            [
+                GeomTypeTable.__table__.c.id,
+                GeomTypeTable.__table__.c.geom.ST_AsText(),
+                GeomTypeTable.__table__.c.geom.ST_SRID(),
+            ],
+        )
+        results = conn.execute(query)
+        rows = results.all()
+
+        # Check that the selected elements are the same as the inputs
+        for row_id, row, srid in rows:
+            checked_wkt = row.upper().replace(" ", "")
+            expected_wkt = inserted_wkt.upper().replace(" ", "")
+            if "MULTIPOINT" in geom_type:
+                # Some dialects return MULTIPOINT geometries with nested parenthesis and others
+                # do not so we remove them before checking the results
+                checked_wkt = re.sub(r"\(([0-9]+)\)", "\\1", checked_wkt)
+            if row_id >= 5 and dialect_name in ["geopackage"] and has_m:
+                # Currently Shapely does not support geometry types with M dimension
+                assert checked_wkt != expected_wkt
+            else:
+                assert checked_wkt == expected_wkt
+            assert srid == 4326
+
     @test_only_with_dialects("postgresql", "sqlite")
     def test_insert_geom_poi(self, conn, Poi, setup_tables):
         conn.execute(
@@ -380,7 +514,7 @@ def test_WKT(self, session, Lake, setup_tables, dialect_name, postgis_version):
         lake = Lake("LINESTRING(0 0,1 1)")
         session.add(lake)
 
-        if (dialect_name == "postgresql" and postgis_version < 3) or dialect_name == "sqlite":
+        if dialect_name == "postgresql" and postgis_version < 3:
             with pytest.raises((DataError, IntegrityError)):
                 session.flush()
         else: