Previously we discussed the basics of the LandXML Surface format and how to read and edit the format. It's a relatively simple format to manipulate using Python. If we look at the primary use of the LandXML Surface format, it's one of the few ways to transfer a surface between AutoCAD's Civil 3D, MicroStation InRoads and ESRI ArcGIS without constructing the surface data from scratch. The three big CAD/GIS vendors all have proprietary TIN formats that can only be read using their software. Also when using LandXML to transfer surfaces, make sure you understand the limitations of InRoads with LandXML as detailed in the previous post.
There is one major issue with the transfer of surface data in the LandXML format with ArcGIS. Both Civil 3D and InRoads have the option to both import and export LandXML data, while ArcGIS is able to import the LandXML format but not export surface data from its TIN format to LandXML. This issue was highlighted on StackExchange and by ESRI users as an ArcGIS Idea.
The good news that we laid out in the previous post is that LandXML is a relatively simple format to read and write using Python and open source tools. Let's use those easily available open source tools and Python to come up with a solution.
When it comes to reading the polygon shapefile, there are several commonly used python libraries that fit the bill. The GDAL library is a great swiss army library that can be downloaded easily by navigating to Chris Gohlke's website and scrolling down to the GDAL library. For windows users this is one of the simplest ways to get a library installed that is compiled with c by using "pip install {insert/your/path/to/wheel/here.whl}".
For this project we are going to keep things simple, which means you can reach for the Python Shapefile library which can be installed with "pip install pyshp". The only two libraries we will use outside of the Python standard library are pyshp and lxml. We covered in the previous LandXML post how to install lxml. The script will be run using Python 3.6, but should also work in 3.5 and possibly 2.7.
The first part of the problem is that the ESRI TIN format is a proprietary format that we cannot directly read. But there is a workaround for this dilemma... See the TIN Triangle tool under the 3D Analyst Tools. We won't go through the code for this as you can either use ArcPy or just the standard ArcGIS dialog to access the tool under ArcToolbox.
When it comes to overall approach, it will be fairly simple and take fewer than 100 lines to write a complete solution. The LandXML format is a collection of a faces/triangles which is represented by the ids of the three points that make up the faces/triangles and another collection of the individual points with their assigned id.
Below is the original TIN file I will be working with.
Below is the polygon created based on the TIN file using the TIN Triangle tool.
Let's get started by adding the necessary imports and reading the input TIN shapefile...
import sys
import os
import shapefile
import lxml.etree as ET
# Input polygon shapefile converted from TIN surface
tin_shp = r'C:\KnickKnackCivil\scripts\ExampleData\TinTriangles.shp'
# Outputs
out_xml = os.path.splitext(tin_shp)[0]+'_Surface.xml'
surf_name = os.path.splitext(os.path.basename(tin_shp))[0]
# Reading input TIN shapefile using PyShp
in_shp = shapefile.Reader(tin_shp)
shapeRecs = in_shp.shapeRecords()
Now we want to define our LandXML surface elements...
# Initializing landxml surface items
landxml = ET.Element('LandXML')
units = ET.SubElement(landxml, 'Units')
surfaces = ET.SubElement(landxml, 'Surfaces')
surface = ET.SubElement(surfaces, 'Surface', name=surf_name)
definition = ET.SubElement(surface, 'Definition',
surfType="TIN")
pnts = ET.SubElement(definition, 'Pnts')
faces = ET.SubElement(definition, 'Faces')
# Change units here if using meters or international feet
unit_imperial = ET.SubElement(units,
'Imperial',
areaUnit="squareFoot",
linearUnit="USSurveyFoot",
volumeUnit="cubicFeet",
temperatureUnit="fahrenheit",
pressureUnit="inHG")
We are going to define a point dictionary and a list of faces with each list item being a collection of the three point ids for each face/triangle.
# Initializing output variables
pnt_dict = {}
face_list = []
cnt = 0
Now we come to the actual creation of the points and triangles. We will iterate through the list of input triangle polygons, and build our previously defined dictionary of point ids for each new unique coordinate. We will also build our list of triangles based on point ids, and then iterate and write those after all of our points have been created.
# Creating reference point dictionary/id for each coordinate
# As well as LandXML points, and list of faces
for sr in shapeRecs:
shape_pnt_ids = [] # id of each shape point
# Each shape should only have 3 points
for pnt in range(3):
# Coordinate with y, x, z format
coord = (sr.shape.points[pnt][1],
sr.shape.points[pnt][0],
sr.shape.z[pnt])
# If element is new, add to dictionary and
# write xml point element
if coord not in pnt_dict:
cnt+=1
pnt_dict[coord] = cnt
shape_pnt_ids.append(cnt) # Add point id to list
# Individual point landxml features
pnt_text = f'{coord[0]:.5f} {coord[1]:.5f} {coord[2]:.3f}'
pnt = ET.SubElement(pnts, 'P', id=str(cnt)).text = pnt_text
# If point is already in the point dictionary, append existing point id
else:
shape_pnt_ids.append(pnt_dict[coord])
# Reference face list for each shape
face_list.append(shape_pnt_ids)
# Writing faces to landxml
for face in face_list:
ET.SubElement(faces, 'F').text = f'{face[0]} {face[1]} {face[2]}'
For the finale, write the defined xml to the output file...
# Writing output
tree = ET.ElementTree(landxml)
tree.write(out_xml, pretty_print=True, xml_declaration=True, encoding="iso-8859-1")
Here is the output LandXML surface visualized in Civil 3D...
And there you have it, a relatively painless landxml TIN surface converted from an ESRI TIN in less than 100 lines using simple, free off the shelf tools... It wasn't that difficult was it ESRI?
If you want the complete script that can be run in the command line see below. No warranty or guarantees come with this script.
'''MIT License
Copyright (c) 2018 David Hostetler
Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:
The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.'''
import sys
import os
import shapefile
import lxml.etree as ET
import argparse
# Argument to provide help with Command Line arguments
parser = argparse.ArgumentParser(description='Creates a LandXML surface file ' \
'based on reference polygon shapefile. ' \
'\nPolgon shapefile should be created from' \
'ESRI TIN Triangle tool.')
parser.add_argument('TIN_shp',help='Input polygon shapefile based on ESRI TIN')
parser.add_argument('Units',
choices=['ft','m','ft-int'],
help=('Input polygon shapefile units. Options are \n' \
'us survey feet (ft), meters (m) or ' \
'international feet (ft-int).'))
parser.parse_args()
tin_shp = sys.argv[1]
unit_len = sys.argv[2]
# Outputs
out_xml = os.path.splitext(tin_shp)[0]+'_Surface.xml'
surf_name = os.path.splitext(os.path.basename(tin_shp))[0]
# Reading input TIN shapefile using PyShp
in_shp = shapefile.Reader(tin_shp)
shapeRecs = in_shp.shapeRecords()
# Initializing landxml surface items
namespace = {'xsi' : "http://www.w3.org/2001/XMLSchema"}
landxml = ET.Element('LandXML',
nsmap=namespace,
xmlns="http://www.landxml.org/schema/LandXML-1.2",
language = 'English',
readOnly = 'false',
time = '08:00:00',
date = '2019-01-01',
version="1.2")
units = ET.SubElement(landxml, 'Units')
surfaces = ET.SubElement(landxml, 'Surfaces')
surface = ET.SubElement(surfaces, 'Surface', name=surf_name)
definition = ET.SubElement(surface, 'Definition',
surfType="TIN")
pnts = ET.SubElement(definition, 'Pnts')
faces = ET.SubElement(definition, 'Faces')
# Dictionary to define correct units based on input
unit_opt = {'ft':('Imperial', 'squareFoot', 'USSurveyFoot',
'cubicFeet', 'fahrenheit', 'inHG'),
'm': ('Metric', 'squareMeter', 'meter',
'cubicMeter', 'celsius', 'mmHG'),
'ft-int': ('Imperial', 'squareFoot', 'foot',
'cubicFeet', 'fahrenheit', 'inHG')}
# Define units here. Has not been tested with metric.
unit = ET.SubElement(units,
unit_opt[unit_len][0],
areaUnit=unit_opt[unit_len][1],
linearUnit=unit_opt[unit_len][2],
volumeUnit=unit_opt[unit_len][3],
temperatureUnit=unit_opt[unit_len][4],
pressureUnit=unit_opt[unit_len][5])
# Initializing output variables
pnt_dict = {}
face_list = []
cnt = 0
print('Processing...')
# Creating reference point dictionary/id for each coordinate
# As well as LandXML points, and list of faces
for sr in shapeRecs:
shape_pnt_ids = [] # id of each shape point
# Each shape should only have 3 points
for pnt in range(3):
# Coordinate with y, x, z format
coord = (sr.shape.points[pnt][1],
sr.shape.points[pnt][0],
sr.shape.z[pnt])
# If element is new, add to dictionary and
# write xml point element
if coord not in pnt_dict:
cnt+=1
pnt_dict[coord] = cnt
shape_pnt_ids.append(cnt) # Add point id to list
# Individual point landxml features
pnt_text = f'{coord[0]:.5f} {coord[1]:.5f} {coord[2]:.3f}'
pnt = ET.SubElement(pnts, 'P', id=str(cnt)).text = pnt_text
# If point is already in the point dictionary, append existing point id
else:
shape_pnt_ids.append(pnt_dict[coord])
# Check if too many or too few points created
if len(shape_pnt_ids) != 3:
print('Error - check input shapefile. '\
'Must be a polygon with only three nodes for each shape.')
sys.exit(0)
# Reference face list for each shape
face_list.append(shape_pnt_ids)
# Writing faces to landxml
for face in face_list:
ET.SubElement(faces, 'F').text = f'{face[0]} {face[1]} {face[2]}'
# Writing output
tree = ET.ElementTree(landxml)
tree.write(out_xml, pretty_print=True, xml_declaration=True, encoding="iso-8859-1")
print(f'Successfully created {out_xml}. Check file output.')
NOTE: An update was made to the script based on user feedback about the validation scheme.
There are comments.