ogr2ogr, PostGIS and Spatialite: create 3D vectors and transfer geometries Z values to attributes table

Sample data: dataset3dvectors.zip

Create a 3D shapefile from a 2D one

first check what ogrinfo says about the sample data:

$ ogrinfo -so shapefile_2d.shp shapefile_2d
INFO: Open of `shapefile_2d.shp'
      using driver `ESRI Shapefile' successful.

Layer name: shapefile_2d
Geometry: Point
Feature Count: 5
Extent: (484628.515471, 4227667.038043) - (617462.831823, 4605830.183414)
Layer SRS WKT:
PROJCS["WGS_1984_UTM_Zone_29N",
    GEOGCS["GCS_WGS_1984",
        DATUM["WGS_1984",
            SPHEROID["WGS_84",6378137,298.257223563]],
        PRIMEM["Greenwich",0],
        UNIT["Degree",0.017453292519943295]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",-9],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["Meter",1]]
z_value: Integer (10.0)

Now generate a real 3D shapefile by transferring the Z values stored in the “z_value” column in the attribute table to the output geometries:

$ ogr2ogr -f "ESRI Shapefile" shapefile_3d.shp shapefile_2d.shp -zfield z_value

Check the result with ogrinfo:

$ ogrinfo -so shapefile_3d.shp shapefile_3d
INFO: Open of `shapefile_3d.shp'
      using driver `ESRI Shapefile' successful.

Layer name: shapefile_3d
Geometry: 3D Point
Feature Count: 5
Extent: (484628.515471, 4227667.038043) - (617462.831823, 4605830.183414)
Layer SRS WKT:
PROJCS["WGS_1984_UTM_Zone_29N",
    GEOGCS["GCS_WGS_1984",
        DATUM["WGS_1984",
            SPHEROID["WGS_84",6378137,298.257223563]],
        PRIMEM["Greenwich",0],
        UNIT["Degree",0.017453292519943295]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",-9],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["Meter",1]]
z_value: Integer (10.0)

The same would apply also when the input layer is a line or a polygon one.

Transfer Z values from a 3D shapefile to the table of attributes of a 2D one

points:

$ ogr2ogr shapefile_with_z_from_3d.shp shapefile_3d.shp -dialect sqlite -sql "SELECT *, st_z(geometry) AS z_geom FROM shapefile_3d"

force output as 2D, add the “-dim 2” flag

$ ogr2ogr shapefile_with_z_from_3d.shp shapefile_3d.shp -dim 2 -dialect sqlite -sql "SELECT *, st_z(geometry) AS z_geom FROM shapefile_3d"

lines:

$ ogr2ogr shapefile_lines_with_z_from_3d.shp shapefile_3d_lines.shp -dim 2 -dialect sqlite -sql "SELECT *, st_z(st_pointn(geometry,1)) AS z_geom FROM shapefile_3d_lines"

polygons:

$ ogr2ogr shapefile_polygons_with_z_from_3d.shp shapefile_3d_polygons.shp -dim 2 -dialect sqlite -sql "SELECT *, st_z(st_pointn(ST_ExteriorRing(geometry),1)) AS z_geom FROM shapefile_3d_polygons"

Rationale: for polygons is first needed to get the exterior ring (as a linestring) of the geometries with the ST_ExteriorRing function, and then for both lines and polygons is needed to get Nth point (in this examples the first one) of the linestring with the ST_PointN function. After that the Z value is extracted with ST_Z.

Import a 2D Shapefile into PostGIS as 3D vector

$ ogr2ogr -progress --config PG_USE_COPY YES -f PostgreSQL PG:"host=*** port=5432 dbname=*** password=*** user=curso" /path/to/shapefile_2d.shp shapefile_2d -lco SCHEMA=public -lco GEOMETRY_NAME=geom -lco FID=gid -nln 3d_vector -zfield z_value

check the result with ogrinfo:

$ ogrinfo -so PG:"host=** user=** dbname=** password=***" 3d_vector                                      INFO: Open of `PG:host=localhost user=curso dbname=curso password=curso'
      using driver `PostgreSQL' successful.

Layer name: 3d_vector
Geometry: 3D Point
Feature Count: 5
Extent: (484628.515471, 4227667.038043) - (617462.831823, 4605830.183414)
Layer SRS WKT:
PROJCS["WGS 84 / UTM zone 29N",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4326"]],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",-9],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    AUTHORITY["EPSG","32629"],
    AXIS["Easting",EAST],
    AXIS["Northing",NORTH]]
FID Column = gid
Geometry Column = geom
z_value: Integer (10.0)

Then from the PostgreSQL console add a new empty column to the PostGIS 3D vector:

ALTER TABLE 3d_vector ADD COLUMN z_geom integer;

and then populate it with the Z values from the geometries:

UPDATE 3d_vector SET z_geom = st_z(geom);

For Spatialite the notes are the same as for PostGIS.

ogr2ogr, PostGIS and Spatialite: create 3D vectors and transfer geometries Z values to attributes table