ConvertgeospatialdataShapefileto3DdataVTKwithPythonGeopandasandPyvista.png

En nuestra perspectiva, la visualización 3D de datos geoespaciales ha sido una característica deseada desde hace mucho tiempo en SIG y que se ha cubierto en algunas características de SAGA GIS o en algunos complementos de QGIS. Esta vez desarrollamos un script en Python que convierte punto / línea / polígono de shapefiles ESRI (o cualquier archivo vectorial) al formato Vtk de grilla no estructurada (Vtu) mediante el uso de las bibliotecas de Python: Geopandas y Pyvista. El tutorial tiene archivos, scripts y videos que muestran todo el procedimiento con algunos comentarios sobre el software y los archivos espaciales y una discusión sobre la naturaleza de los archivos espaciales que presenta algunos desafíos en la conversión de datos.

Notas para este tutorial:

Todas las geometrías 3D están relacionadas con el atributo "Elev" y el campo de geometría del dataframe de Geopandas.
Los shapefiles de líneas deben ser partes individuales.
Los polígonos no deben tener agujeros o no se considerarán.
Los polígonos se exportan como líneas y no como polígonos debido al hecho de que los Vtks se diseñaron para polígonos convexos que tienen un rendimiento deficiente con polígonos geoespaciales.

Como usuarios habituales de formatos de datos Shapefile y Vtk, sabemos que el potencial de conversión es enorme; sin embargo, consideramos que este tutorial es un buen punto de partida para futuras mejoras o la aplicación en casos particulares.

Conozca más del formato VTK en este enlace:

https://lorensen.github.io/VTKExamples/site/VTKFileFormats/

Tutorial

Código

Este es el código en Python relacionado a la conversion:

#import required packages
import itertools
import numpy as np
import pyvista as pv
import geopandas as gpd
#for windows users
from shapely import speedups
speedups.disable()

C:\Users\GIDA2\miniconda3\lib\site-packages\geopandas\_compat.py:84: UserWarning: The Shapely GEOS version (3.4.3-CAPI-1.8.3 r4285) is incompatible with the GEOS version PyGEOS was compiled with (3.8.1-CAPI-1.13.3). Conversions between both will be slow.
  warnings.warn(

#create geodataframes from all shapefiles
pointDf = gpd.read_file('../Shps/wellPoints.shp')
lineDf = gpd.read_file('../Shps/contoursLines.shp')
polyDf = gpd.read_file('../Shps/lakesPolygons.shp')

For point type geometry

#create emtpy lists to collect point information
cellSec = []
cellTypeSec = []
pointSec = []

#iterate over the points
i = 0
for index, valuesx in pointDf.iterrows():
    x, y, z = pointDf.loc[index].geometry.x, pointDf.loc[index].geometry.y, pointDf.loc[index].Elev
    pointSec.append([x,y,z])
    cellTypeSec.append([1])
    cellSec.append([1,i])
    i+=1

#convert list to numpy arrays
cellArray = np.array(cellSec)
cellTypeArray = np.array(cellTypeSec)
pointArray = np.array(pointSec)

#create the unstructured grid object
pointUgrid = pv.UnstructuredGrid(cellArray,cellTypeArray,pointArray)

#we can add some values to the point
pointUgrid.cell_arrays["Elev"] = pointDf.Elev.values

#plot and save as vtk
pointUgrid.plot()
pointUgrid.save('../Vtks/wellPoint.vtu',binary=False)

For line type geometry

#create emtpy dict to store the partial unstructure grids
lineTubes = {}

#iterate over the points
for index, values in lineDf.iterrows():
    cellSec = []
    linePointSec = []

    #iterate over the geometry coords
    zipObject = zip(values.geometry.xy[0],values.geometry.xy[1],itertools.repeat(values.Elev))
    for linePoint in zipObject:
        linePointSec.append([linePoint[0],linePoint[1],linePoint[2]])

    #get the number of vertex from the line and create the cell sequence
    nPoints = len(list(lineDf.loc[index].geometry.coords))
    cellSec = [nPoints] + [i for i in range(nPoints)]

    #convert list to numpy arrays
    cellSecArray = np.array(cellSec)
    cellTypeArray = np.array([4])
    linePointArray = np.array(linePointSec)

    partialLineUgrid = pv.UnstructuredGrid(cellSecArray,cellTypeArray,linePointArray)   
    #we can add some values to the point
    partialLineUgrid.cell_arrays["Elev"] = values.Elev
    lineTubes[str(index)] = partialLineUgrid

#merge all tubes and export resulting vtk
lineBlocks = pv.MultiBlock(lineTubes)
lineGrid = lineBlocks.combine()
lineGrid.save('../Vtks/contourLine.vtk',binary=False)
lineGrid.plot()

For poly type geometry

#create emtpy dict to store the partial unstructure grids
polyTubes = {}

#iterate over the points
for index, values in polyDf.iterrows():
    cellSec = []
    linePointSec = []

    #iterate over the geometry coords
    zipObject = zip(values.geometry.exterior.xy[0],
                    values.geometry.exterior.xy[1],
                    itertools.repeat(values.Elev))
    for linePoint in zipObject:
        linePointSec.append([linePoint[0],linePoint[1],linePoint[2]])

    #get the number of vertex from the line and create the cell sequence
    nPoints = len(list(polyDf.loc[index].geometry.exterior.coords))
    cellSec = [nPoints] + [i for i in range(nPoints)]

    #convert list to numpy arrays
    cellSecArray = np.array(cellSec)
    cellTypeArray = np.array([4])
    linePointArray = np.array(linePointSec)

    partialPolyUgrid = pv.UnstructuredGrid(cellSecArray,cellTypeArray,linePointArray)   
    #we can add some values to the point
    partialPolyUgrid.cell_arrays["Elev"] = values.Elev
    #    partialPolyUgrid.save('../vtk/partiallakePoly.vtk',binary=False)
    polyTubes[str(index)] = partialPolyUgrid

#merge all tubes and export resulting vtk
polyBlocks = pv.MultiBlock(polyTubes)
polyGrid = polyBlocks.combine()
polyGrid.save('../Vtks/lakePolyasLines.vtk',binary=False)
polyGrid.plot()