[GIS] How to convert 3D DXF to PostGIS 2.0

batchconvertdxfpostgis

I'm not very skilled on using postgis but I face a problem that I must solve, i.e. create tables to load geographical data stored as dxf.

I have a thousand dxf files (3d) with polylines, points and polygons and I want to convert them to a postgis/postgres db.

I'm in windows xp using postgres 9 postgis 2.0.

Gdal/OGR are installed using Osgeo4Win installer.
I can write some basic python script to handle data in batch.

Can you please direct me to some useful tutorial/hint on how to face this problem.
Thanks in advance

Best Answer

I'm a little bit late, but hope it helps.

Tuts and books

If you want something learn about ogr, gdal, python and postGIS you can look at first on these pages and book.

Geoprocessing with Python using Open Source GIS

http://www.gis.usu.edu/~chrisg/python/

ogr

http://www.gdal.org/ogr/

ogr2ogr - cheatsheet

http://www.bostongis.com/?content_name=ogr_cheatsheet#41

for an programming example there is an an apitutorial for C++, C and Python

http://www.gdal.org/ogr/ogr_apitut.html

you can also begin here:

http://wiki.osgeo.org/wiki/OSGeo_Python_Library

http://trac.osgeo.org/gdal/wiki/GdalOgrInPython

PostGIS

http://www.bostongis.com/ is allways a good starting point. And also the book "PostGIS in Action" by them is an great guid to PostGIS.

An other nice overview you can find in the following book:

Erik Westra, Python Geospatial Development.

http://www.packtpub.com/python-geospatial-development/book

Now to your specific problem.

First Problem - Geometry

GEOMETRY in dxf vs GEOMETRY in GIS

One problem is the definition of geometries in dxf and in GIS

Geometries in dxf:

Points are allways points, there is no problem here.

Linestrings are linestrings, theres also no problem so far.

But when it comes to polygons there is an problem with in dxf. Polygons will allways represented as LINESTRING!!! Whenever you choose polygon or 3dpoly as option, to create an polygon, it will be an linestring.

If your read the features from an dxf-file with ogr the representation of polygons is ever an representation of an linestring.

Example:

WKT-string: LINESTRING(1 1 1,2 1 1,2 2 1,1 2 1,1 1 1)

But should be:

WKT: POLYGON(1 1 1,2 1 1,2 2 1,1 2 1,1 1 1)

See Wikipedia for WKT-notation

ogr2ogr doesn't change this. You have to do it yourself. We cover this later.

Second Problem - ogr2ogr

The great problem with ogr2ogr is to know what all this options are for and where to learn about. There isn't much help with ogr2ogr --help or ogr2ogr --long-usage on the commandline. So by trial and error I have figured out the following.

ogr2ogr -dsco INIT_WITH_EPSG=yes -lco LAUNDER=yes -lco SPATIAL_INDEX=yes \
    -append -skipfailures -f PostgreSQL \
    PG:"dbname={dbname} user={user} password={password} schemas={schema}" \
    -nln {yourtable} \
    -nlt {geometry}25D \
    -a_srs EPSG:{EPSG-Code}
    {srs-file}.dxf

Here are the options (for Version 1.9):

  • dsco INIT_WITH_EPSG -- don't sure we need this for PostgreSQL, works for Spatialite
  • lco -- layer creation option; LAUNDER=yes all strings to lower and some special charackters will be presented as underbar
  • lco SPATIAL_INDEX=yes; creates an spatial index
  • append (or/and update) -- append to existing table
  • skipfailures -- use this option, otherwise ogr2ogr will stop if it find an geometry that isn't defined in the table
  • f PostgreSQL -- write to PostgreSQL
  • nln -- yourtablename
    • create a table with this name in your db, without! columns, this will ogr2ogr do, it takes the dxf-structure (layer, subclasses, extendedentity, linetype, entityhandle, text)
    • if you you want your own colums, you have to figure out it by yourself
    • without nln it creates an table "entities"
  • nlt -- new layer type; here you can choose between an 2D or an 3D geometry
    • from ogr2ogr --long-usage
    • nlt type: Force a geometry type for new layer. One of NONE, GEOMETRY, POINT, LINESTRING, POLYGON, GEOMETRYCOLLECTION, MULTIPOINT, MULTIPOLYGON, or MULTILINESTRING. Add "25D" for 3D layers Default is type of source layer.
    • default representation for dxf-features is 2D!, but with three coordinates, so ogr2ogr sets Z-values allways to zero
  • a_srs EPSG: -- assign an output EPSG-code

You have to do this for points, linestrings and geometrycollections.

Then you can write an shell script like this one (created by underdark)

for f in `ls *dxf` do
ogr2ogr -dsco INIT_WITH_EPSG=yes -lco LAUNDER=yes -lco SPATIAL_INDEX=yes \
  -append -skipfailures -f PostgreSQL \
  PG:"dbname={dbname} user={user} password={password} schemas={schema}" \
  -nln $f \
  -nlt {geometry}25D \
  -a_srs EPSG:{EPSG-Code}
  $f
done

So then you have an table with points, an table with linestrings and an table with geometrycollections. Now you have to figure out where your polygons are. You can find them with the following SQL-statement. It returns all linestrings where the first point and the last point are equal. You can put this in an SELECT INTO-statement and create an new table. After that you can also clean your linestring table with this statement. It is also possible to create an trigger so it will be a little more comfortable ;-).

SELECT ogc_fid, layer,
       ST_Polygon(ST_GeomFromEWKT(ST_asText(wkb_geometry)), {EPSG-Code}) as the_geom
INTO {your-polygon-table} 
FROM {your-linestring-table}
WHERE ST_StartPoint(wkb_geometry) = ST_EndPoint(wkb_geometry);
Related Question