Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

I´m trying to make a 3D shapefile from a raster feature but not working. Tried already to run the command in OSGeo4W but nothing happen. I'm using the Qgis 1.8.0 Lisboa. The purpose is to get a final 3D dxf file! Can anyone help me?

share|improve this question
    
I'm unsure, but I don't think 1.8.0 handles 3D geometry. –  Ryan Garnett Apr 3 '13 at 20:12
    
is there any version which can handle 3D geometry? –  Jorge Marques Apr 3 '13 at 20:13
    
From what I have read, the develoent version of 1.9 or 2.0 is suppose to. I know the the newest version of PostGIS will handle 3D geometry –  Ryan Garnett Apr 3 '13 at 20:44

1 Answer 1

For the moment, the geometry system of QGIS does not handle z-values but z-coordinates can be read by parsing the binary format of the geometry (see Why is it not possible to extract the z value of a wkb geometry with PyQGIS while it is possible with ogr or Shapely ?) of which I take the example:

With a 3D point shapefile in the Python console of QGIS:

mylayer = qgis.utils.iface.activeLayer()
mylayer.geometryType() == QGis.Point
True
# select the first object
sel = mylayer.selectedFeatures()[0]
qgisgeom = sel.geometry()
qgisgeom.wkbType() == QGis.WKBPoint
False
qgisgeom.wkbType() == QGis.WKBPoint25D
True
# QGIS recognizes that the WKB geometry is 3D but is is impossible to extract the z value 
qgisgeom.asPoint()
(205553,89857.7)
# no z in the functions of geom.asPoint, only x and y in dir()
dir(qgisgeom.asPoint())
['__class__', ..., 'x', 'y']

But using other Python modules like osgeo.ogr or Shapely, you can extract the z value from the binary format of the QGIS geometry:

# with shapeky
from shapely.wkb import loads
wkb = qgisgeom.asWkb()
transfshapely = loads(wkb)
transfshapely.has_z
True
list(transfshapely.coords)
[(205552.628666262, 89857.66932836606, 222.1999969482422)]

# with osgeo.ogr
from osgeo import ogr
wkb = qgisgeom.asWkb()
geom_wkb = ogr.CreateGeometryFromWkb(wkb)
geom_wkb.GetGeometryType() == ogr.wkbPoint25D
True
geom_wkb.GetX(),geom_wkb.GetY(),geom_wkb.GetZ()
(205552.628666262, 89857.66932836606, 222.1999969482422)

So, for create a 3D shapefile in the Python console or in the Python shell you must use modules that are able to handle the z value like ogr, Fiona with Shapely or Pyshp (CreateElevationValues)

Example with Fiona:

import fiona
from shapely.geometry import mapping, Point
# example of a 3D point (Shapely)
point = Point(0, 0, 3)
# schema of the shapefile
schema = {'geometry': '3D Point','properties': {'id': 'int'},}
with fiona.collection('my_shp3D.shp', 'w', 'ESRI Shapefile', schema) as output:
   ## If there are multiple geometries, put the "for" loop here
   feat = {}
   feat['geometry'] = mapping(point)
   feat['properties'] = {'id': 1}
   output.write(feat)

Control:

d = fiona.open('my_shp3D.shp')
d.schema
{'geometry': '3D Point', 'properties': {u'id': 'int'}}
d.next()
{'geometry': {'type': 'Point', 'coordinates': (0.0, 0.0, 3.0)}, 'id': '0', 'properties': {'id': 1}}

As an example of extracting the z value from a DEM and the color value from a geological map, see "Python: Using vector and raster layer in a geological perspective without GIS software", Python: utilisation des couches vectorielles et matricielles dans une perspective géologique, sans logiciel SIG , in French but the Python scripts are universal...

share|improve this answer
    
My friends, thank you for your support! I've jus reached the result i wanted! Got the 3D dxf file from a shapefile! First we must reproject the raster image to projected coordinates, then run the raster tool extract contours, but in the expression shown below you must add "-3d" after the atribute used for elevation values. Run the tool!Then I opened it directly in Microstation and did a simple save as dxf file! –  Jorge Marques Apr 4 '13 at 16:25

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.