ConvertgeospatialdataShapefileto3DdataVTKwithPythonGeopandasandPyvista.png

In our perspective, 3D visualization of geospatial data has been a long desired feature that has been covered in some features from SAGA GIS or in some plugins from QGIS. This time we developed a Python script that converts point / line / polygon ESRI shapefiles (or any vector file) to unstructured grid Vtk format type (Vtu) by the use of the Python libraries Geopandas and Pyvista. The tutorial has files, scripts, and videos that show the whole procedure with some remarks on the software and spatial files and a discussion about the nature of the spatial files that presents some challenges in the data conversion.

Notes for this tutorial:

All 3d geometries are related to the "Elev" attribute and the geometry field of the Pandas dataframe.
Lines shapefiles must be single parts
Polygons should not have holes or they wont be considered.
Polygons are exported as lines not as polygons due to the fact that Vtks were designed for convex polygons having poor performance with geospatial polygons.

As regular Shapefile and Vtk data format users, we know that the potential for the conversion is huge; however we consider this tutorial a good starting point for future improvement or the application in particular cases.

Know more about the VTK format on this link:

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

Tutorial

Scripts

This is the Python code related to the 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()