Geoalchemy: insert 3D elements in PostGIS database

Question

I have shp files that I would like to insert into a PostGIS database. I'm planning to use Geoalchemy2 and geopandas to do so. However, my data have the geometry types LINESTRING Z and POINT Z. As I don't know the geometry type of object, I used Geometry("GEOMETRY") dtype as parameter of my to_sql call.

I get the following error as a result:

(psycopg2.DataError) Geometry has Z dimension but column does not

Here is a functioning code:

from geoalchemy2 import Geometry, WKTElement
import geopandas as gpd
import sqlalchemy

engine = sqlalchemy.create_engine('postgresql://username:password@host:port/username')

gdf = gpd.read_file("file.shp")
gdf['geom'] = gdf['geometry'].apply(lambda x: WKTElement(x.wkt))
gdf.drop('geometry', 1, inplace=True)
gdf.to_sql(
    "tablename",
    engine,
    if_exists='append',
    index=False,
    dtype={'geom': Geometry('Geometry', dimension=3)}
)

The dimension=3 parameter doesn't change anything.

I saw PostGIS Column has Z dimension but geometry does not, Geojson in postgis Geometry : Z dimension error and PostGIS error reads "Geometry has Z dimension but column does not" propositions but as I create the tables with the to_sql command, I cannot alter the fields. I tried to alter the column afterwards and to insert then the 3D elements (e.g. LINESTRING Z (...)) and that did work.

Do you think I should:

• Create the tables manually first by reading the columns?
• Stop looking for Python shp import and use something else (shp2pgsql)?
• Any other idea?

Answer 1

I would like to give my own experience about it :

I wanted to import shapefile from Qgis in Postgres, but i needed to have 3D data in postgres regardless of the type in Qgis.

I mean, the input shapefile could be a 2D geometric dataset with an attribute Z or a 3D geometric dataset. And in the end we should have 3D data in our postgres database.

So I used shp2pgsql to import the shapefile in temporary schema called "work" with :

1) import 3D or import the 2D shapefile. i used a bit of python script before the call to shp2pgsl to check if there is a z attribute and raise a warning if no possibility is found to have an elevation. then, you could ask your user on Qgis to create the attribute, depends on what you want. if the user go for a full 2D dataset, i raise a warning and try :

1 - use a DEM to find an elevation

2 - go for St_Force_3D of the 2D data (will give default 0 value for Z).

2) then with some SQL request re-pull the data in the database, with request inside the database.

to show some code, here is the python function which allow you to build the part of the SQL statement which handle the geometry.

        if __json_keys__[table]['geom']['typeimport'] == 'ZPOINT' :
            if self.geomtype == 'ZPOINT' or self.geomtype == 'ZMPOINT' :
                geomstatement = "ST_SetSRID(ST_MakePoint(ST_X({point}),ST_Y({point}),ST_Z({point})),ST_SRID({point}))".format(point = self.connection_table.cellWidget(0,1).layout().itemAt(0).widget().currentText())
            elif self.geomtype == 'POINT' :
                row_z_ground = [i+1 for i,x in enumerate(__json_keys__[table]["columns"]) if x == "z_ground"]
                if len(row_z_ground) == 1:
                    row_z_ground = row_z_ground[0]
                    if self.connection_table.cellWidget(row_z_ground,1).layout().itemAt(0).widget().currentText() :
                        geomstatement = "ST_SetSRID(ST_MakePoint(ST_X({point}),ST_Y({point}),{z_data}),ST_SRID({point}))".format(point = self.connection_table.cellWidget(0,1).layout().itemAt(0).widget().currentText(),z_data=self.connection_table.cellWidget(row_z_ground,1).layout().itemAt(0).widget().currentText())
                    else :
                        geomstatement = "ST_SetSRID(ST_MakePoint(ST_X({point}),ST_Y({point}),project.altitude(ST_MakePoint(ST_X({point}),ST_Y({point})))),ST_SRID({point}))".format(point = self.connection_table.cellWidget(0,1).layout().itemAt(0).widget().currentText())
                else :
                    geomstatement = "ST_SetSRID(ST_Force_3DZ(ST_MakePoint(ST_X({point}),ST_Y({point}))),ST_SRID({point}))".format(point = self.connection_table.cellWidget(0,1).layout().itemAt(0).widget().currentText())
        elif __json_keys__[table]['geom']['typeimport'] == 'MULTILINESTRINGZ' :
            geomstatement = "ST_Force_3DZ((ST_Dump("+ self.connection_table.cellWidget(0,1).layout().itemAt(0).widget().currentText() +")).geom)"
        return geomstatement

IMO shp2pgsql is a good option if you want to import .shp in postgres