You are here:GeoTux»Geo-Blogs»Desktop GIS»QGIS, visualización 3D de capas vectoriales con Python

Statistics

Invitados: 41
Usuarios registrados: 3153
Usuarios en línea:
-
Registrados hoy:
-

Register

RSS

Blogs and News:
Recibe las actualizaciones en Geo-Noticias y Geo-Blogs

Get them by e-mail
Recibir Geo-Noticias y Geo-Blogs por e-mail

¿What is this about?

Latest Geo-Tweets

Sunday, 21 October 2012 12:58

QGIS, visualización 3D de capas vectoriales con Python

Written by  martin Laloux
Rate this item
(8 votes)

En este artículo muestro cómo visualizar en 3D capas vectoriales (bien sean shapefiles de tipo 3D o shapefiles con atributo de altitud) con los módulos Python Shapely, Matplotlib y Visvis desde la consola Python de Quantum GIS.

Hace poco preguntaron (de nuevo) en la lista de desarrollo de QGIS ¿qué se puede hacer con el valor de altitud (z) de los shapefiles de tipo 3D (o 2.5D, PointZ, LineZ, PolygonZ, ...) o de los shapefiles con un atributo de altitud? ¿Es posible visualizarlos en 3 dimensiones?

 

Presenté una solución en la lista QGIS-developer, osgeo-org.1560.n6.nabble.com/Re-Qgis-use...alues-td5004272.html y osgeo-org.1560.n6.nabble.com/visualizing...ution-td5005360.html. También publiqué una versión más completa en francés en www.portailsig.org/content/qgis-represen...apefiles-avec-attrib. La presento aquí con el acuerdo del Portail SIG www.portailsig.org/.

 

El caso de la representación 3D en QGIS es una pregunta recurrente en los foros y listas de correo. Las soluciones propuestas hasta ahora son indirectas, pasando por el uso de GRASS GIS, PostGIS o SpatiaLite (www.faunalia.com/content/transfer-3d-shapefiles-z-values-table-attributes), o esperando por la extensión Globe (Globe), aún en desarollo.

  • El caso de las capas con un atributo de altura (z) es fácil de tratar. El plugin Interpolation de QGIS permite generar grids (formato ASCII ArcInfo grid, .asc), pero no ofrece ninguna forma  de visualizar los resultados en 3D.
  • El problema de los llamados shapefiles 3D, donde la valor de altitud z está integrada en la geometría, proviene de que QGIS o PyQGIS no saben tratarlos, los consideran como shapefiles 2D.

Sin embargo, se puede representar estas capas en 3D con los módulos Python que permiten dicha representación:

 

  1. Para el tratamiento de los shapefiles de tipo 3D, existen otros módulos de Python que son capaces de extraer el valor z, como ogr (de GDAL http://pypi.python.org/pypi/GDAL/1.9.1) y Shapely (http://pypi.python.org/pypi/Shapely/1.2.16), directamente desde una capa de QGIS, o Pyshp (http://code.google.com/p/pyshp/), desde un fichero shapefile. Estos módulos se pueden utilizar en la consola Python sin problema.

  2. Para la visualización en 3 dimensiones, vamos a utilizar 2 módulos Python que funcionan también muy bien en la consola:
    1. El clásico módulo Matplotlib http://matplotlib.org/).
    2. Un nuevo módulo, Visvis http://code.google.com/p/visvis/). Utiliza  NumPy y PyOpenGL.
    3. No hablaremos ni utilizaremos el potente módulo Mayavi http://code.enthought.com/projects/mayavi), demasiado complejo para lo que queremos hacer.

El problema de los shapefiles de tipo 3D en Quantum GIS

 

QGIS y PyQGIS reconocen el tipo 3D pero no ofrecen ninguna solución para utilizarlo:

Shapefile 3D (puntos, PointZ) que vamos a utilizar para el tratamiento:

Tratamiento con PyQGIS

 

Se puede observar que no hay ninguna manera de extraer la coordenada z de un punto3D (PointZ) aunque PyQGIS reconozca que es una geometria 3D (geom.wkbType()).

Tratamiento con otros módulos Python (ogr y Shapely)

Conversión de la geometría del punto en una geometría ogr (con el formato WKB):


 

Conversión de la geometría del punto en una geometría Shapely (con el formato WKB):

 

En los dos casos, vemos que es posible recuperar el valor de la coordenada z para utilizarla posteriormente.


Tratamientos en la consola Python

 

Los usuarios de QGIS que quieren trabajar en la consola Python empiezan con los tutoriales  http://www.qgis.org/pyqgis-cookbook/index.htmlhttp://www.qgisworkshop.org/html/workshop/index.html. Pero estos dos tutoriales, sin querer ofender a los autores, resultan bastante complejos para el principiante, especialmente en los procesos de iteración.

Por ejemplo, el proceso "iterating over a Vector Layer" se puede hacer en manera mucho mas sencilla en la consola Python:

  • Lo primero es seleccionar todos los elementos (geometría y atributos) de una capa. Se puede hacer con una simple función (sin QgsFeature = () = allAttrs provider.attributeIndexes (), provider.select (allAttrs) y otros ...):
def select_todo(capa):
    capa.select([])
    capa.setSelectedFeatures([obj.id() for obj in capa])
  • A continuación, se activa la capa de vector que queremos tratar y se utiliza la función anterior para seleccionar todos los elementos:
micapa = qgis.utils.iface.activeLayer()
select_todo(micapa)
  • Como resultado, todas las geometrías y todos los atributos de la capa están en la variable micapa.

El caso de las geometrías de tipo 3D

Las coordenadas x, y, z se extraen con la función loads de Shapely:

from shapely.wkb import loads
# lectura de todas la geometrías de la capa y extracción de las coordenadas  x,y,z
x=[]
y=[]
z=[]
for elem in micapa.selectedFeatures():
           # exportación de la geometrías al formato WKB
           geom= elem.geometry() 
           wkb = geom.asWkb()
           # extracción de las coordenadas x,y,z con Shapely 
           x.append(loads(wkb).x)
           y.append(loads(wkb).y)
           z.append(loads(wkb).z)

o una sugestión de Etienne Tourigny para mejorar el proceso cuando hay muchos elementos:

           data = loads(wkb)
           x.append(data.x)
           y.append(data.y)
z.append(data.z)

El caso de las geometrías 2D con un atributo de altitud

PyQGIS sabe extraer la coordinadas x,y (no hace falta otro módulo).

En la capa utilizada, la variable de altitud esta en segunda posición (elev).

  • El valor se adquiere con elem.attribute.map () [1] (se podría mejorar el proceso utilizando el nombre del atributo (elev), pero no es el propósito aquí).
  • .toFloat()[0] es la transformación de un objeto Qvariant de PyQt4 en una variable Python utilizable por los otros módulos.
# lectura de las coordenadas x,y  y extracción de los valores del atributo de elevación
x=[]
y=[]
z=[]
for elem in micapa.selectedFeatures():
           geom= elem.geometry() 
           x.append(geom.asPoint()[0])
           y.append(geom.asPoint()[1])
           z.append(elem.attributeMap()[1].toFloat()[0])

 

En los dos casos, los resultados son tres listas Python con la coordenadas x,y,z que se pueden utiliza para todos los tratamientos.


Visualización en 3D (puntos, grids, superficies , ASCII grid (.asc))

1) Los puntos

Es muy fácil utilizar estas 3 listas con Matplotlib

from mpl_toolkits.mplot3d.axes3d import *
import matplotlib.pyplot as plt
from matplotlib import cm
fig = plt.figure()
ax = Axes3D(fig)
ax.scatter3D(x,y,z,c=z,cmap=plt.cm.jet)
plt.show()

 

y el resultado es:

 

o con Visvis:

import visvis
# abrir un ventana y visualizar los puntos
f = visvis.gca()
m = visvis.plot(x,y,z, lc='k', ls='', mc='g', mw=2, lw=2, ms='.')
f.daspect = 1,1,10 # z x 10

 

resultado:

 

2) Grid, superficies a partir de los puntos

Como Matplotlib, NumPy o SciPy disponen de la mayoría de la funciones de interpolación, es muy sencillo calcular grids o superficies a partir de los puntos.

 

Con la función griddata de de Matplotlib

Es la más fácil de usar, utiliza un algoritmo de Delaunay o del "vecino más cercano" (según la instalación de Matplotlib y la repartición de los puntos) para interpolar una malla o rejilla de puntos (ver por ejemplo http://pybonacci.wordpress.com/2012/04/13/dibujando-lineas-de-nivel-en-python-con-matplotlib/)

import numpy as np
from matplotlib.mlab import griddata
# creación de una grilla 2D 
xi = np.linspace(min(x), max(x))
yi = np.linspace(min(y), max(y))
X, Y = np.meshgrid(xi, yi)
# interpolación
Z = griddata(x, y, z, xi, yi)

 

Resultado con Matplotlib:

fig = plt.figure()
ax = Axes3D(fig)
ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap=cm.jet,linewidth=1, antialiased=True)
plt.plot

 

Resultado con Visvis (usando también el griddata de Matplotlib):

f = visvis.gca()
m = visvis.grid(xi,yi,Z) # o m= =visvis.grid(yi,xi,Z) con una antigua versión 
f.daspect = 1,1,10 # z x 10

 

Grid con los puntos:

 

Relieve sombreado con una tabla de color (visvis.CM_JET)

m = visvis.surf(xi,yi,Z)
m.colormap = visvis.CM_JET

Con una de las funciónes spline de SciPy

Como ejemplo de las posibilidades, vamos a utilizar otro algoritmo de interpolación, el de "spline en placa delgada" (thin-plate spline) de SciPy.

import scipy as sp
import scipy.interpolate
# construcción de la grilla
spline = sp.interpolate.Rbf(x,y,z,function='thin-plate')
xi = np.linspace(min(x), max(x))
yi = np.linspace(min(y), max(y))
X, Y = np.meshgrid(xi, yi)
# interpolación
Z = spline(X,Y)

 

Resultado con Visvis

f = visvis.gca()
m = visvis.surf(xi,yi,Z)
m.colormap = visvis.CM_JET 
f.axis.visible = False

 

 

3) Con los grids en formato ESRI ASCII grid (.asc) resultantes del tratamiento del plugin Interpolation de QGIS.

En Python, hay dos maneras de leer estos ficheros:

  1. Con el módulo gdal.
  2. Con los resultados de los trabajos de David Finlayson que ofrecen varias soluciones para leer los ASCII grid o los grids de Golden Surfer. Las ultimas versiones se encuentran en bitbucket.org/david_finlayson/pyraster/src

Como son grillas regulares, hay varias funciones de Numpy o SciPy para tratarlas.

 

Resultado con el plugin Interpolation de QGIS:

 

Resultado del tratamiento en la consola Python con Visvis


 

Otros tratamientos posibles

 

Se pueden efectuar otros tratamientos como visualizar las curvas de nivel o de isolíneas en 3D (el plugin Contour de QGIS, también basado sobre Matplotlib, lo hace en 2D).

 

Todas las líneas y polígonos 3D pueden ser representados con Visvis:

 

Y como soy un geólogo, puedo representar los sondeos geológicos (y eventualmente generar las superficies de cada contacto).

Y ademas de la visualización, se pueden salvar los grids resultantes usando las funciones provistas por David Finlayson o utilizando NumPy / SciPy.


Comparación con los resultados de GRASS GIS

 

En ningún momento de estos tratamientos se habla de las proyecciones de las capas, excepto a través de las propias coordenadas. Por tanto, es interesante comparar los resultados obtenidos con los de GRASS GIS con los mismos datos, en un sistema de proyección:

resultado con Visvis (z x 1)                        resultado con GRASS GIS (z x 1)

resultado con Visvis (z x 10)                                      resultado con GRASS GIS (z x 10)

 

Las diferencias principales no resultan del hecho de que los datos sean proyectados o no (en un sistema de proyección cartesiano, x, y), sino del algoritmo de interpolación, del coeficiente de exageración de las alturas o del tipo de perspectiva.


Conclusiones

 

Esta claro que estos métodos no pueden sustituir totalmente el uso de aplicaciones como GRASS GIS o del plugin Globe aún en desarrollo.

Figura de www.faunalia.com/quantumgis

 

Pero son muy útiles como primer tratamiento, es decir, permiten visualizar directamente lo que se quiere con los puntos seleccionados y/ o combinar diferentes selecciones sin tener que lanzar  GRASS GIS.

¿Por qué no un plugin?

Hacer un plugin me parece muy difícil porque:

Todos los tratamientos han sido realizados con Mac OS X y las versiones  "master branch"  de QGIS de Larry Shaffer, disponible en qgis.dakotacarto.com/ y 1.8 de William Kyngesburye (www.kyngchaos.com/)


Last modified on Saturday, 03 November 2012 11:28

comments  

 
0 # no z coordinateenzo 2013-11-05 10:48
hi,

I tried this system but in win7 I have this problem:

raise DimensionError( "This point has no z coordinate.")

why?

If I use in linux It work

thanks!

enzo
Reply | Reply with quote | Quote
 
 
0 # no z coorddinategene 2013-11-13 17:12
I do not know because I don't use Windows.
But, my colleagues working on Windows have never had any problem.
Sorry.
Reply | Reply with quote | Quote
 
 
0 # visualizar en 3D desde *.tifxunilk 2013-11-13 11:05
Hola, Martin

Me interesó tanto tu artículo que lo puse en práctica considerando la lectura de los arrays x,y,z desde un GeoTIFF (*.tif). Mi artículo es el siguiente:

http://joseguerreroa.wordpress.com/2013/11/13/como-extraer-las-coordenadas-x-y-z-de-un-raster-tif-para-visualizacion-3d-con-python/

Saludos
Reply | Reply with quote | Quote
 
 
0 # this point has no z coordinateenzo 2014-03-20 00:35
hi martin,

I continue to have this error on every SO

"this point has no z coordinate"
line 67 in Shapley/geometry/point.py

do you have any idea?
Reply | Reply with quote | Quote
 

Add comment


Security code
Refresh

 

On-line users

Latest Geo-Forums

More Topics »

Latest Comments