ogr2ogr doesn't convert 3D and 4D geometries correct

[horst.duester@…]

The conversion of 3D and 4D geometries leads to bad results. You can download test data from:

 http://www.sogis1.so.ch/sogis/dl/aww_gewiso_dump.zip

This file contains a PostGIS dump of a 4D Linestring Layer and a shape created with pgsql2shp.

ogr2ogr is not able to create a correct 4D or 3D shape from the attached PostGIS layer.

ogrinfo PG:dbname=.... aww_gewiso

reports:

Layer name: aww_gewiso
Geometry: 3D Multi Line String
Feature Count: 6609
Extent: (589671.375000, 212041.000000) - (645546.750000, 261743.703125)
Layer SRS WKT:
(unknown)
FID Column = ogc_fid
Geometry Column = wkb_geometry
....

The conversion with ogr2ogr from PostGIS to shape produces a shape without geometry.

[mloskot]

The URL given in the report does not work anymore, file is unavailable.

[mloskot]

Perhaps, this issue is directly related or even duplicate of the Ticket #1323

[mloskot]

I contacted directly with Horst asking him for new URL to the data package and if he could verify if the issue still occurs after recent fixes related to PostGIS.

[mloskot]

I got reply from Horst. Here it is with my inlined response:

Düster Horst wrote:

Hi Mateusz Please excuse the dead link. I just reactivated it.

Horst,

Thanks, I downloaded it.

I made a short check and the bug seems to persist ogr 1.4.0.

Yes, this version is still affected.

Since which version your mentioned fixes will be available?

1.4.1 and higher.

Here you can see related reports already fixed:

• 1036 XYMZ data in WKT leads to an error

• 1195 Naked ewkb from postigs (queried through python) causes OGR to segfault

• 1323 XYZM data cannot be read from wkb

[mloskot]

Here is Horst's reply:

Düster Horst wrote:

I tried it with gdal-1.4.1 and it works with restriction. ogr2ogr produces from 4d table a 3d shape. Try it with the download files.

---

I have the test data and I'm going to test it soon.

[mloskot]

I finally tested current GDAL SVN and export of PostGIS 3D/4D data to Shapefile works almost correctly. Here are steps I took while testing:

1. I imported Horst's data from his Shapefile (unfortunately, the aww_gewiso.dmp is useless due to missing data/schema dump):

shp2pgsql -W LATIN1 aww_gewiso_pgsql2shp.shp gewiso | psql -d gewiso

2. Export from PostGIS back to Shapefile using ogr2ogr

mloskot@dog:~$ ogr2ogr -f "ESRI Shapefile" out.shp PG:dbname=gewiso
OGR_PG: DBName="gewiso"
OGR_PG: PostSIS version string: '1.3 USE_GEOS=1 USE_PROJ=1 USE_STATS=1' -> '1.3'
OGR_PG: POSTGIS_VERSION=1.3 USE_GEOS=1 USE_PROJ=1 USE_STATS=1
OGR_PG: Using column 'ogc_fid' as FID for table 'gewiso'
OGR_PG: Layer 'gewiso' geometry type: MULTILINESTRING:3D Multi Line String, Dim=3
OGR: OGROpen(PG:dbname=gewiso/0x804f870) succeeded as PostgreSQL.
OGR_PG: PQexec(DECLARE OGRPGLayerReader CURSOR for SELECT "ogc_fid", AsEWKT("the_geom"), "abschnitt", "anno", "archive", "archive_da", "arcid25040", "bearbeiter", "bearbeitun", "bende", "dh_gew", "dh_km", "erfasser", "erfassung", "flip", "gem2", "gemeinde", "gewissnr", "gid", "gnrso", "groesse", "hauptgnr", "meas_oben", "meas_unten", "name", "nebengnr", "neg_imp", "new_date", "privat", "qualitaet", "se_connect", "se_select", "strahler", "symbol", "typ", "vorf_gewis", "vorf_neben" FROM "gewiso" )
PG: 12031 features read on layer 'gewiso'.

Note, the geometry is detected as 3 DIM, not 4.

3. Compare out.shp with original file aww_gewiso_pgsql2shp.shp using shpdump from Shapelib

• original data

mloskot@dog:~/dev/gdal/bugs/1229$ shpdump aww_gewiso_pgsql2shp.shp | head -n 10
Shapefile Type: ArcZ   # of Shapes: 12031

File Bounds: (  589671.375,  212041.000,0,0)
         to  (  645546.733,  261743.703,116424,116424)

Shape:0 (ArcZ)  nVertices=2, nParts=1
  Bounds:(  621683.642,  237491.198, 192.528, 192.528)
      to (  621712.874,  237557.632, 265.109, 265.109)
     (  621712.874,  237491.198, 192.528, 192.528) Ring 
     (  621683.642,  237557.632, 265.109, 265.109)  

Note, what coordinates we have of the 4th DIM.

• exported data - out.shp

mloskot@dog:~/dev/gdal/bugs/1229$ shpdump out.shp | head -n 10
Shapefile Type: ArcZ   # of Shapes: 12031

File Bounds: (  589671.375,  212041.000,0,0)
         to  (  645546.733,  261743.703,116424,0)

Shape:0 (ArcZ)  nVertices=2, nParts=1
  Bounds:(  621683.642,  237491.198, 192.528, 0)
      to (  621712.874,  237557.632, 265.109, 0)
     (  621712.874,  237491.198, 192.528, 0) Ring 
     (  621683.642,  237557.632, 265.109, 0)  

As you see, the original shapefile is not exactly the same as the exported shapefile (out.shp). In exported file, coordinates of M dimension are Zeros.

Original shapefile is correctly detected and imported by shp2pgsql:

gewiso=> select * from geometry_columns ;
 f_table_catalog | f_table_schema | f_table_name | f_geometry_column | coord_dimension | srid |      type       
-----------------+----------------+--------------+-------------------+-----------------+------+-----------------
                 | public         | gewiso       | the_geom          |               4 |   -1 | MULTILINESTRING

I'm checking OGR PG driver what's going wrong.

[mloskot]

According to my understanding of the PG driver code, 4 DIM is intentionally not supported and if a 4 DIM is detected for a layer, it's processed as 3 DIM. Here is the code: ogrpgtablelayer.cpp:382

Frank,

AFAIR, we've decided (see comments of #1323) to not to touch the PG driver code now and wait for official specification of WKB with M dimension well defined. If my understanding is correct, then I believe we can close this ticked or move to future after 1.5.

Could you confirm it?

[mloskot]

I'm closing this bug for now.