Are there any tools that can take a DTM, in asc or geotiff format and generate a mesh from this that can used in meshlab or similar software.

I have asc DTM file that I need to convert to a ply mesh file. Is this possible?

In response to the comment. Then these DTM asc files are not the ones that meshlab understands.

The files starts with:

ncols        6250
nrows        6250
xllcorner    630000.000000000000
yllcorner    6070000.000000000000
cellsize     1.600000000000
NODATA_value  -9999
  • 1
    MeshLab can read .asc grids – gene Nov 10 '14 at 15:07
  • Do you need the resulting mesh to consist of triangles or quadrilaterals? In the case of triangles, you're looking for a way to convert a raster to a TIN (triangular irregular network). Do you want a naive triangulation (two triangles per raster cell), or do you want to approximate the terrain by using fewer triangles? – Jake Nov 10 '14 at 21:06
  • I will take the first tool that get the job done and then continue searching for others if it dont work. I want to import it into a 3th party tool that imports obj,ply and a bunch of other mesh formats. I assume that the format will be containing the specifications, and just load it, without me actually having to make a decision on the question you ask. – Poul K. Sørensen Nov 10 '14 at 23:09
  • @pksorensen: Fair enough. I think you might not get away with completely avoiding decisions, however. Even though they're used for similar things, raster and mesh formats have some very fundamental differences (similar to raster and vector images), so there's no "default" way for converting between the two. What are you going to do with the mesh (3D printing, visualisation, finite element simulation, ...)? – Jake Nov 11 '14 at 6:02
  • A orthomap generation where rays are intersected with the terrain model. But this is done in a tool that just import a large set of mesh files, so I assume that the tool will handle the import and internal figure out how to read the mesh to its internal storage that enable to to do the orthomap. – Poul K. Sørensen Nov 11 '14 at 12:30
up vote 10 down vote accepted

A DTM raster can be represented by triangle meshes by finding a set of non-overlapping triangles that covers the entire mesh and approximates the elevation field. There are two different types of triangle meshes that can be used for this purpose:

  • a triangulated regular network (TRN), in which every pixel of the raster is represented by a vertex, and all triangles have the same size and shape. All the original information of the DTM raster is present in the TRN, but the memory required for storing the mesh is typically quite high.
  • a triangulated irregular network (TIN), in which there are fewer vertices than raster pixels and the triangles have different shapes and sizes. The vertices and the triangulation are chosen in such a way that the resulting surface approximates the original DTM raster up to a specified error. This typically results in much smaller files, since plane or nearly plane areas can be represented using only a couple of vertices.

In most applications, if you need to deal with elevation meshes, you'd go with a TIN since throwing out redundant or nearly redundant information allows for more efficient computations. However, creating TINs from rasters isn't straightforward, since there are many different triangulations that approximate a grid with the same error, but using different vertex sets.

Software for creating TINs

  • Michael Garland's Terra software.
  • ArcGIS: Raster to TIN function from the 3D Analyst toolbox.
  • SAGA GIS: Grid to TIN (Surface Specific Points) function, followed by Export TIN to Stereo Lithography File (STL) function to export the TIN to a mesh format readable by Meshlab.

Software for creating TRNs

  • VTBuilder, which is part of the Virtual Terrain Project: load the DTM raster using "Layer | Import Layer" and then convert it to a TRN using "Elevation | Convert grid to TIN". Then you select the newly generated layer in the Layer overview and select "Elevation | Export To...". VTBuilder can read all raster formats that GDAL supports, and exports the TRN to OBJ, PLY, GMS, DXF, DAE or WRL formats.
  • SAGA GIS: Grid to TIN function.
  • Roll your own solution: Implementing it isn't particularly hard. Here's a Python script that uses the GDAL library to read a raster DTM, and then writes out a binary PLY mesh.

    Save the script as gdal_rastertotrn.py, then call it using python gdal_rastertotrn.py <inputraster> <outputply>.

    Here's an example. Converting a raster DTM of Crater Lake called crater_lake.tif...

    ... by calling python gdal_rastertotrn.py crater_lake.tif crater_lake.ply, and opening the resulting crater_lake.ply in Meshlab:

    Here's the (unpolished) script. Since it uses the GDAL library, it can convert all raster types supported by GDAL. It only writes PLY files.

    #!/usr/bin/python
    
    import sys
    import numpy as np
    from osgeo import gdal
    
    def write_ply(filename, coordinates, triangles, binary=True):
        template = "ply\n"
        if binary:
            template += "format binary_" + sys.byteorder + "_endian 1.0\n"
        else:
            template += "format ascii 1.0\n"
        template += """element vertex {nvertices:n}
    property float x
    property float y
    property float z
    element face {nfaces:n}
    property list int int vertex_index
    end_header
    """
    
        context = {
         "nvertices": len(coordinates),
         "nfaces": len(triangles)
        }
    
        if binary:
            with  open(filename,'wb') as outfile:
                outfile.write(template.format(**context))
                coordinates = np.array(coordinates, dtype="float32")
                coordinates.tofile(outfile)
    
                triangles = np.hstack((np.ones([len(triangles),1], dtype="int") * 3,
                    triangles))
                triangles = np.array(triangles, dtype="int32")
                triangles.tofile(outfile)
        else:
            with  open(filename,'w') as outfile:
                outfile.write(template.format(**context))
                np.savetxt(outfile, coordinates, fmt="%.3f")
                np.savetxt(outfile, triangles, fmt="3 %i %i %i")
    
    def readraster(filename):
        raster = gdal.Open(filename)
        return raster
    
    
    def createvertexarray(raster):
        transform = raster.GetGeoTransform()
        width = raster.RasterXSize
        height = raster.RasterYSize
        x = np.arange(0, width) * transform[1] + transform[0]
        y = np.arange(0, height) * transform[5] + transform[3]
        xx, yy = np.meshgrid(x, y)
        zz = raster.ReadAsArray()
        vertices = np.vstack((xx,yy,zz)).reshape([3, -1]).transpose()
        return vertices
    
    
    def createindexarray(raster):
        width = raster.RasterXSize
        height = raster.RasterYSize
    
        ai = np.arange(0, width - 1)
        aj = np.arange(0, height - 1)
        aii, ajj = np.meshgrid(ai, aj)
        a = aii + ajj * width
        a = a.flatten()
    
        tria = np.vstack((a, a + width, a + width + 1, a, a + width + 1, a + 1))
        tria = np.transpose(tria).reshape([-1, 3])
        return tria
    
    
    def main(argv):
        inputfile = argv[0]
        outputfile = argv[1]
    
        raster = readraster(inputfile)
        vertices = createvertexarray(raster)
        triangles = createindexarray(raster)
    
        write_ply(outputfile, vertices, triangles, binary=True)
    
    if __name__ == "__main__":
        main(sys.argv[1:])
    

Your Answer

 
discard

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Not the answer you're looking for? Browse other questions tagged or ask your own question.