Re: [Qgis-user] Using/visualizing 3D data (Z values)

classic Classic list List threaded Threaded
12 messages Options
Reply | Threaded
Open this post in threaded view
|

Re: [Qgis-user] Using/visualizing 3D data (Z values)

Etienne Tourigny-3
Hi all

reviving this old thread because I'm interested in using the z-value (elevation) from kml files (type wkbMultiLineString25D) in a qgis plugin.

Are there any plans to introduce z-value (25D) support in 2.0?

Experimenting has shown that doing geom.exportToWkt() strips the z-value.

How can I currently handle z-values? I understand that I would have to parse the wkb directly, is this possible in python? If so can someone give me some pointers?

Other option - external python library like shapely? It seems shapely does support z levels.

Or would it be simpler to parse the kml files by hand using an xml parser, as syntax is relatively simple?


regards,
Etienne

On Mon, Dec 5, 2011 at 8:02 AM, Marco Hugentobler <[hidden email]> wrote:












Hi Paul

Unfortunately, the current geometry system does not handle z- and 
m-values in a sound way.
Some tools consider z and others don't. M-values are even more difficult 
as they cannot be accessed by the QGIS API (z-coordinates can be read by 
parsing the binary format of the geometry).

An overhaul of the geometry system is planned for version 2.

Regards,
Marco

On 05.12.2011 10:42, Paul Lens wrote:
> Hi all,
>
> 1. I made a TIN surface, using the interpolation plugin, from a point 
> shapefile (with "ALTITUDE" field) and  3D line shapefiles. I lost time 
> trying to extract the relevant lines from the line shapefiles without 
> loosing the 3D information (Z values of the nodes).
> If I remember well, I first used "Merge Layers" from the mmqgis 
> plugin, and it kept the Z information. Than, I used the Intersection 
> command in the Vector menu, and the Z information was lost. Finally, 
> after manually selecting the relevant lines and using "save as", I got 
> the expected 3D line shapefile keeping the Z values.
> Is this loss of Z information a known limitation of the gdal vector 
> commands in QGIS (vector menu) and is this limitation affecting all 
> these commands? Or did I do something wrong?
>
> 2. In a second step, I wanted to be able to easily see (through 
> labeling a Z field) or ask (with the identify command) those Z values 
> from the 3D lines, so as to fine tune the interpretation of the 
> topography.
> Again the gdal tool (in Vector menu) "Node Extraction" lost all Z 
> information.
> Is there a workaround to create a point shapefile with a Z field? Or, 
> at least, to keep the Z values in a point shapefile, rasterize it so 
> as to be able to visualize the Z value using  the "Identify" command 
> on the raster?
>
> Thanking all the developers for this wonderful program,
>


-- 
Dr. Marco Hugentobler
Sourcepole -  Linux&  Open Source Solutions
Churerstrasse 22, CH-8808 Pfäffikon SZ, Switzerland
[hidden email] http://www.sourcepole.ch
Technical Advisor QGIS Project Steering Committee


_______________________________________________
Qgis-user mailing list
[hidden email]
http://lists.osgeo.org/mailman/listinfo/qgis-user



_______________________________________________
Qgis-developer mailing list
[hidden email]
http://lists.osgeo.org/mailman/listinfo/qgis-developer
Reply | Threaded
Open this post in threaded view
|

Re: [Qgis-user] Using/visualizing 3D data (Z values)

gene
Hello, In Python:
- For shapefiles, , I use the pyshp (pyshp: Python Shapefile Library/) module of Joel Lawhead to read or write 3D shapefiles:
(see CreateElevationValues  : How to work with elevation (a.k.a. z) values.
>>>import shapefile
>>>r = shapefile.Reader("MyPolyZ")
>>>r.shapes()[0].points
[[-89.0, 33.0], [-90.0, 31.0], [-91.0, 30.0]]
>>>r.shapes()[0].z
[12, 11, 12]

I use it without problem in the Python console  to process the results with Shapely, Numpy or Grass.script

- For kml files, it is relatively easy to parse files with standard modules like Elementtree,  or with more specialized modules (http://pypi.python.org/pypi?%3Aaction=search&term=kml&submit=search)

I also use this approach in the Python console
Reply | Threaded
Open this post in threaded view
|

Re: [Qgis-user] Using/visualizing 3D data (Z values)

Etienne Tourigny-3
Thanks for the info.

I was able to use shapely pretty easily, by passing the wkb from qgis
to shapely via shapely.wkb.loads(wkb)

This removed the need for actual file parsing, which qgis seem to do
well (apart from the z-level values, which it reads but does not
report),

here is the relevant python code:

from shapely.wkb import loads
feat = QgsFeature()
provider = layer.dataProvider()
provider.select(provider.attributeIndexes())
if provider.nextFeature(feat):
  lines = loads(feat.geometry().asWkb())
  if len(lines) == 1:
    coords = lines[0].coords
    (x,y,z[i]) = coords[i]

cheers
Etienne

On Wed, Sep 26, 2012 at 3:26 AM, gene <[hidden email]> wrote:

> Hello, In Python:
> - For shapefiles, , I use the pyshp ( pyshp: Python Shapefile Library/
> <http://code.google.com/p/pyshp/>  ) module of Joel Lawhead to read or write
> 3D shapefiles:
> (see  CreateElevationValues  : How to work with elevation (a.k.a. z) values.
> <http://code.google.com/p/pyshp/wiki/CreateElevationValues>
>>>>import shapefile
>>>>r = shapefile.Reader("MyPolyZ")
>>>>r.shapes()[0].points
> [[-89.0, 33.0], [-90.0, 31.0], [-91.0, 30.0]]
>>>>r.shapes()[0].z
> [12, 11, 12]
>
> I use it without problem in the Python console  to process the results with
> Shapely, Numpy or Grass.script
>
> - For kml files, it is relatively easy to parse files with standard modules
> like Elementtree,  or with more specialized modules (
> http://pypi.python.org/pypi?%3Aaction=search&term=kml&submit=search
> <http://pypi.python.org/pypi?%3Aaction=search&term=kml&submit=search>  )
>
> I also use this approach in the Python console
>
>
>
> --
> View this message in context: http://osgeo-org.1560.n6.nabble.com/Re-Qgis-user-Using-visualizing-3D-data-Z-values-tp5004272p5004487.html
> Sent from the Quantum GIS - Developer mailing list archive at Nabble.com.
> _______________________________________________
> Qgis-developer mailing list
> [hidden email]
> http://lists.osgeo.org/mailman/listinfo/qgis-developer
_______________________________________________
Qgis-developer mailing list
[hidden email]
http://lists.osgeo.org/mailman/listinfo/qgis-developer
Reply | Threaded
Open this post in threaded view
|

Re: [Qgis-user] Using/visualizing 3D data (Z values)

Nathan Woodrow
Very handy.  Is there any issues with exposing the Z values from QGIS?
 If it's there why not let people get to it.

- Nathan

On Wed, Sep 26, 2012 at 9:42 PM, Etienne Tourigny
<[hidden email]> wrote:

> Thanks for the info.
>
> I was able to use shapely pretty easily, by passing the wkb from qgis
> to shapely via shapely.wkb.loads(wkb)
>
> This removed the need for actual file parsing, which qgis seem to do
> well (apart from the z-level values, which it reads but does not
> report),
>
> here is the relevant python code:
>
> from shapely.wkb import loads
> feat = QgsFeature()
> provider = layer.dataProvider()
> provider.select(provider.attributeIndexes())
> if provider.nextFeature(feat):
>   lines = loads(feat.geometry().asWkb())
>   if len(lines) == 1:
>     coords = lines[0].coords
>     (x,y,z[i]) = coords[i]
>
> cheers
> Etienne
>
> On Wed, Sep 26, 2012 at 3:26 AM, gene <[hidden email]> wrote:
>> Hello, In Python:
>> - For shapefiles, , I use the pyshp ( pyshp: Python Shapefile Library/
>> <http://code.google.com/p/pyshp/>  ) module of Joel Lawhead to read or write
>> 3D shapefiles:
>> (see  CreateElevationValues  : How to work with elevation (a.k.a. z) values.
>> <http://code.google.com/p/pyshp/wiki/CreateElevationValues>
>>>>>import shapefile
>>>>>r = shapefile.Reader("MyPolyZ")
>>>>>r.shapes()[0].points
>> [[-89.0, 33.0], [-90.0, 31.0], [-91.0, 30.0]]
>>>>>r.shapes()[0].z
>> [12, 11, 12]
>>
>> I use it without problem in the Python console  to process the results with
>> Shapely, Numpy or Grass.script
>>
>> - For kml files, it is relatively easy to parse files with standard modules
>> like Elementtree,  or with more specialized modules (
>> http://pypi.python.org/pypi?%3Aaction=search&term=kml&submit=search
>> <http://pypi.python.org/pypi?%3Aaction=search&term=kml&submit=search>  )
>>
>> I also use this approach in the Python console
>>
>>
>>
>> --
>> View this message in context: http://osgeo-org.1560.n6.nabble.com/Re-Qgis-user-Using-visualizing-3D-data-Z-values-tp5004272p5004487.html
>> Sent from the Quantum GIS - Developer mailing list archive at Nabble.com.
>> _______________________________________________
>> Qgis-developer mailing list
>> [hidden email]
>> http://lists.osgeo.org/mailman/listinfo/qgis-developer
> _______________________________________________
> Qgis-developer mailing list
> [hidden email]
> http://lists.osgeo.org/mailman/listinfo/qgis-developer
_______________________________________________
Qgis-developer mailing list
[hidden email]
http://lists.osgeo.org/mailman/listinfo/qgis-developer
Reply | Threaded
Open this post in threaded view
|

Re: [Qgis-user] Using/visualizing 3D data (Z values)

gene
In reply to this post by Etienne Tourigny-3
The problem is that Shapely is based on GEOS and GEOS isn't a 3D library at all. The z of the geometries are not considered during any spatial operations but you can use matplotlib to represent the 3D geometry

>>> from shapely.geometry import Point
>>> p = geometry.Point(1,2,3)
>>> list(p.coords)
[(1.0, 2.0, 3.0)]
>>> p.has_z
True
>>> p.z
3.

and wkt et wkb formats do not recognize the z value
>>> p.wkt
'POINT (1.0000000000000000 2.0000000000000000)'

The geo_interface format (GeoJson), yes
>>> from shapely.geometry import mapping, shape
>>> p.__geo_interface__
{'type': 'Point', 'coordinates': (1.0, 2.0, 3.0)}


Reply | Threaded
Open this post in threaded view
|

Re: [Qgis-user] Using/visualizing 3D data (Z values)

Noli Sicad
On 9/26/12, gene <[hidden email]> wrote:
> The problem is that Shapely is based on GEOS and GEOS isn't a 3D library at
> all. *The z of the geometries are not considered during any spatial
> operations but you can use matplotlib to represent the 3D geometry*

Spatialite 4.0 has 3D and it is using GEOS library for geometry. I
wonder where Spatialite gets its 3D?

~~~~

POINT 2D, XY
AddGeometryColumn('tbl', 'geom', 4326, 'POINT', 2);
AddGeometryColumn('tbl', 'geom', 4326, 'POINT', 'XY')
AddGeometryColumn('tbl', 'geom', 4326, 'POINT');

3D, XYZ
AddGeometryColumn('tbl', 'geom', 4326, 'POINT', 3);
AddGeometryColumn('tbl', 'geom', 4326, 'POINT', 'XYZ')
AddGeometryColumn('tbl', 'geom', 4326, 'POINTZ')
AddGeometryColumn('tbl', 'geom', 4326, 'POINTZ', 3);
AddGeometryColumn('tbl', 'geom', 4326, 'POINTZ', 'XYZ');

2D + measure, XYM
AddGeometryColumn('tbl', 'geom', 4326, 'POINT', 'XYM');
AddGeometryColumn('tbl', 'geom', 4326, 'POINTM');
AddGeometryColumn('tbl', 'geom', 4326, 'POINTM', 'XYM');

3D + measure, XYZM
AddGeometryColumn('tbl', 'geom', 4326, 'POINT', 4);
AddGeometryColumn('tbl', 'geom', 4326, 'POINT', 'XYZM')
AddGeometryColumn('tbl', 'geom', 4326, 'POINTZM');
AddGeometryColumn('tbl', 'geom', 4326, 'POINTZM', 4);
AddGeometryColumn('tbl', 'geom', 4326, 'POINTZM', 'XYZM');

https://www.gaia-gis.it/fossil/libspatialite/wiki?name=switching-to-4.0

Noli
_______________________________________________
Qgis-developer mailing list
[hidden email]
http://lists.osgeo.org/mailman/listinfo/qgis-developer
Reply | Threaded
Open this post in threaded view
|

Re: [Qgis-user] Using/visualizing 3D data (Z values)

gene
It is a problem of Shapely. According to the Shapely manual (http://toblerity.github.com/shapely/manual.html):
    "A third z coordinate value may be used when constructing instances, but has no effect on geometric analysis. All operations are performed in the x-y plane"

Nothing prevents you, like me, to use something other with the resulting numpy.arrays or the GeoJson representation of Shapely.
Reply | Threaded
Open this post in threaded view
|

Re: [Qgis-user] Using/visualizing 3D data (Z values)

Etienne Tourigny-3
In reply to this post by gene
On Wed, Sep 26, 2012 at 9:24 AM, gene <[hidden email]> wrote:
> The problem is that Shapely is based on GEOS and GEOS isn't a 3D library at
> all. *The z of the geometries are not considered during any spatial
> operations but you can use matplotlib to represent the 3D geometry*
>

QGis is also based on GEOS.

I'm not talking about doing 3d operations, just having z level
available.  Currently it is read and copied internally (as I
experienced by exporting wkb to shapely), but not output to wkt and
others.  QgsPoint does not contain z value, but we could add a z() and
has_z() members like in shapely.  Enabling full 3d support is another
ball game entirely.

>>>> from shapely.geometry import Point
>>>> p = geometry.Point(1,2,3)
>>>> list(p.coords)
> [(1.0, 2.0, 3.0)]
>>>> p.has_z
> True
>>>> p.z
> 3.
>
> *and wkt et wkb formats do not recognize the z value*
>>>> p.wkt
> 'POINT (1.0000000000000000 2.0000000000000000)'
>

You are right, strictly speaking 3d points should be represented with
wkt POINTZ(x,y,z).

But perhaps we can be practical and allow for a simple extension?  We
could expose the z-level in the api, and *perhaps* expose it to the
wkt (since it's already exposed in the wkb).

This seems to be what is done in ogr (as it reports geometries with z values).

$ ogrinfo  2012-09-23\ 1359__20120923_1359.kml  Tracks
Had to open data source read-only.
INFO: Open of `2012-09-23 1359__20120923_1359.kml'
      using driver `KML' successful.

Layer name: Tracks
Geometry: 3D Multi Line String
Feature Count: 1
Extent: (-45.900415, -23.194376) - (-45.741626, -23.085176)
Layer SRS WKT:
[...]
Name: String (0.0)
Description: String (0.0)
OGRFeature(Tracks):0
  Name (String) = Segment 1
  Description (String) =
MULTILINESTRING ((-45.7416262 -23.08517574 561.3199,-45.74180873
-23.08526774 560.202, [...] ))


> *The geo_interface format (GeoJson), yes*
>>>> from shapely.geometry import mapping, shape
>>>> p.__geo_interface__
> {'type': 'Point', 'coordinates': (1.0, 2.0, 3.0)}
>
>
>
>
>
>
> --
> View this message in context: http://osgeo-org.1560.n6.nabble.com/Re-Qgis-user-Using-visualizing-3D-data-Z-values-tp5004272p5004608.html
> Sent from the Quantum GIS - Developer mailing list archive at Nabble.com.
> _______________________________________________
> Qgis-developer mailing list
> [hidden email]
> http://lists.osgeo.org/mailman/listinfo/qgis-developer
_______________________________________________
Qgis-developer mailing list
[hidden email]
http://lists.osgeo.org/mailman/listinfo/qgis-developer
Reply | Threaded
Open this post in threaded view
|

Re: [Qgis-user] Using/visualizing 3D data (Z values)

gene
In reply to this post by gene
Using a script like:

def select_all(couche):
    couche.select([])
    couche.setSelectedFeatures([obj.id() for obj in couche])

macouche = qgis.utils.iface.activeLayer()
select_all(macouche)

from shapely.wkb import loads
from shapely.geometry import shape
lignes = []
for elem in macouche.selectedFeatures():
    geom= elem.geometry()
    wkb = geom.asWkb()
    lignes.append(loads(wkb))
    list(lignes.coords)
   
You have all the elements (x,y and z) to plot your geometry in 3D with matplotlib, visvis or others or to use the Python grass.script module to create a 3D layer in GRASS


Reply | Threaded
Open this post in threaded view
|

Re: [Qgis-user] Using/visualizing 3D data (Z values)

Nathan Woodrow
That's cool.  It would be cool if there was a QGIS plugin that took
the selected features and rendered them in 3D using matplotlib.... :)

- Nathan

On Wed, Sep 26, 2012 at 11:26 PM, gene <[hidden email]> wrote:

> Using a script like:
>
> def select_all(couche):
>     couche.select([])
>     couche.setSelectedFeatures([obj.id() for obj in couche])
>
> macouche = qgis.utils.iface.activeLayer()
> select_all(macouche)
>
> from shapely.wkb import loads
> from shapely.geometry import shape
> lignes = []
> for elem in macouche.selectedFeatures():
>     geom= elem.geometry()
>     wkb = geom.asWkb()
>     lignes.append(loads(wkb))
>     list(lignes.coords)
>
> You have all the elements (x,y and z) to plot your geometry in 3D with
> matplotlib, visvis or others or to use the Python grass.script module to
> create a 3D layer in GRASS
>
>
>
>
>
>
> --
> View this message in context: http://osgeo-org.1560.n6.nabble.com/Re-Qgis-user-Using-visualizing-3D-data-Z-values-tp5004272p5004652.html
> Sent from the Quantum GIS - Developer mailing list archive at Nabble.com.
> _______________________________________________
> Qgis-developer mailing list
> [hidden email]
> http://lists.osgeo.org/mailman/listinfo/qgis-developer
_______________________________________________
Qgis-developer mailing list
[hidden email]
http://lists.osgeo.org/mailman/listinfo/qgis-developer
Reply | Threaded
Open this post in threaded view
|

Re: [Qgis-user] Using/visualizing 3D data (Z values)

gene
This solution works in the Python console

x=[]
y=[]
z=[]
for elem in macouche.selectedFeatures():
           geom= elem.geometry()
           wkb = geom.asWkb()
           x.append(loads(wkb).x)
           y.append(loads(wkb).y)
           z.append(loads(wkb).z)

For 3d scatter points
from mpl_toolkits.mplot3d.axes3d import *
import pylab as plt
fig = plt.figure()
ax = fig.gca(projection='3d')
ax.plot(x,y,z,'o')
plt.show()
or
fig = plt.figure()
ax = Axes3D(fig)
ax.scatter3D(x,y,z,c=z,cmap=plt.cm.jet)
plt.show()



and with a simple grid from the points
from matplotlib.mlab import griddata
fig = plt.figure()
ax = fig.gca(projection='3d')
xi = np.linspace(min(x), max(x))
yi = np.linspace(min(y), max(y))
X, Y = np.meshgrid(xi, yi)
Z = griddata(x, y, z, xi, yi)
surf = ax.plot_surface(X, Y, Z, rstride=2, cstride=2, cmap=cm.jet,linewidth=1, antialiased=True)
ax.scatter3D(x,y,z,c=z,cmap=plt.cm.jet)
plt.show()


       
No color here because it is a  "discontinuous" surface (griddata).
Reply | Threaded
Open this post in threaded view
|

Re: [Qgis-user] Using/visualizing 3D data (Z values)

Vincent Picavet (ml)
In reply to this post by Etienne Tourigny-3
Hi all,

Le mardi 25 septembre 2012 15:47:14, Etienne Tourigny a écrit :
> reviving this old thread because I'm interested in using the z-value
> (elevation) from kml files (type wkbMultiLineString25D) in a qgis plugin.
>
> Are there any plans to introduce z-value (25D) support in 2.0?

We at Oslandia are working on 3D processing in databases and visualization.

We have an early piece of code opening a 3D vector layer in QGIS and feeding
it to QGIS Globe. This works with real 3D objects loaded from database, and
should work with 2.5D (x, y, z) shapefiles too theoritically.

We have to get it in a better shape before starting to provide pull requests
on this, but you should see this coming in a couple of months.

Then this will open a whole new range of needs with 3D data, but one step
after the other...

Vincent


>
> Experimenting has shown that doing geom.exportToWkt() strips the z-value.
>
> How can I currently handle z-values? I understand that I would have to
> parse the wkb directly, is this possible in python? If so can someone give
> me some pointers?
>
> Other option - external python library like shapely? It seems shapely does
> support z levels.
>
> Or would it be simpler to parse the kml files by hand using an xml parser,
> as syntax is relatively simple?
>
>
> regards,
> Etienne
>
> On Mon, Dec 5, 2011 at 8:02 AM, Marco Hugentobler <
>
> [hidden email]> wrote:
> > Hi Paul
> >
> > Unfortunately, the current geometry system does not handle z- and
> > m-values in a sound way.
> > Some tools consider z and others don't. M-values are even more difficult
> > as they cannot be accessed by the QGIS API (z-coordinates can be read by
> > parsing the binary format of the geometry).
> >
> > An overhaul of the geometry system is planned for version 2.
> >
> > Regards,
> > Marco
> >
> > On 05.12.2011 10:42, Paul Lens wrote:
> > > Hi all,
> > >
> > > 1. I made a TIN surface, using the interpolation plugin, from a point
> > > shapefile (with "ALTITUDE" field) and  3D line shapefiles. I lost time
> > > trying to extract the relevant lines from the line shapefiles without
> > > loosing the 3D information (Z values of the nodes).
> > > If I remember well, I first used "Merge Layers" from the mmqgis
> > > plugin, and it kept the Z information. Than, I used the Intersection
> > > command in the Vector menu, and the Z information was lost. Finally,
> > > after manually selecting the relevant lines and using "save as", I got
> > > the expected 3D line shapefile keeping the Z values.
> > > Is this loss of Z information a known limitation of the gdal vector
> > > commands in QGIS (vector menu) and is this limitation affecting all
> > > these commands? Or did I do something wrong?
> > >
> > > 2. In a second step, I wanted to be able to easily see (through
> > > labeling a Z field) or ask (with the identify command) those Z values
> > > from the 3D lines, so as to fine tune the interpretation of the
> > > topography.
> > > Again the gdal tool (in Vector menu) "Node Extraction" lost all Z
> > > information.
> > > Is there a workaround to create a point shapefile with a Z field? Or,
> > > at least, to keep the Z values in a point shapefile, rasterize it so
> > > as to be able to visualize the Z value using  the "Identify" command
> > > on the raster?
> > >
> > > Thanking all the developers for this wonderful program,
> >
> > --
> > Dr. Marco Hugentobler
> > Sourcepole -  Linux&  Open Source Solutions
> > Churerstrasse 22, CH-8808 Pfäffikon SZ,
> > [hidden email] http://www.sourcepole.ch
> > Technical Advisor QGIS Project Steering Committee
> >
> >
> >
> > _______________________________________________
> > Qgis-user mailing list
> > [hidden email]
> > http://lists.osgeo.org/mailman/listinfo/qgis-user
_______________________________________________
Qgis-developer mailing list
[hidden email]
http://lists.osgeo.org/mailman/listinfo/qgis-developer