Triangulation using X, Y, Z data file to create an elevation model

Question

I have a raster based height map, which I can convert to X,Y,Z heights, i.e. I can get a CSV (Comma Separated Value) file, which has X, Y and Z values.

Since it is a raster, the sampling of the X,Y,Z values would be uniform. Can I create a python script that would create a rectangular mesh with the same resolution as the raster map, and then I modify the Z-coordinate of the vertices from the height map?

Edit: Coding developed using GDAL

With the help of @zelfii in one of the answers here, I have developed a code, which can read the GeoTIFF raster data (SRTM based elevation model) and can create a mesh. However, I am facing some memory problems (Blender is taking too long to process the data)

The size of the raster in question is 4500 x 4500 (approximately).

import bpy 
from osgeo import gdal 
from gdalconst import *
from math import floor
import easygui 

def getFacetIndex(c, r, cols): 
    return (c + r * cols)

def getFaceVertexList(geoData):
    faces = []
    vertices = []

    cols = geoData.RasterXSize
    rows = geoData.RasterYSize 
    transform = geoData.GetGeoTransform()
    xOrigin = transform[0]
    yOrigin = transform[1]
    pixelWidth = transform[1]
    pixelHeight = transform[5]

    band = geoData.GetRasterBand(1)

    verts_per_side_x = cols

    for i in range(cols * rows):
        x = (i % (cols))
        y = floor(i / cols)

        level = x + (y * verts_per_side_x)
        idx1 = level
        idx2 = level + 1
        idx3 = level + verts_per_side_x + 1
        idx4 = level + verts_per_side_x
        if (x < cols - 1) & (y < rows - 1):
            faces.append([idx1, idx2, idx3, idx4])

        xc = xOrigin + pixelWidth * x; 
        yc = yOrigin + pixelHeight * y; 
        data = band.ReadAsArray(x,y,1,1)
        zc = data[0,0]
        if (zc <= 0):
            zc = 340
        vertices.append((xc, yc, zc))

    return vertices, faces


def generate(image_name, geometry_name):
    theData = gdal.Open(image_name, GA_ReadOnly)
    verts, faces = getFaceVertexList(theData) #   

    mesh = bpy.data.meshes.new(geometry_name + "_mesh")
    mesh.from_pydata(verts, [], faces)
    mesh.update()

    obj = bpy.data.objects.new(geometry_name, mesh)
    bpy.context.scene.objects.link(obj)
    obj.select = True
    obj.show_bounds = True


#open the file dialog 
theFile = easygui.fileopenbox()
#open the GDAL file 
gdal.AllRegister() 

generate(theFile, "theModel")

Answer 1

raster to csv to 3d surface

in a CSV all you need to provide are the z-heights per column / row. Decide on a formula to convert pixel intensity (RGB) to the z-height. You do not need to specify x,y components associated with the z-height, because they are implied by the position of the z value in the csv row/column.

At most you need to provide the spread distance between rows and columns, once. Given you are working from a raster image, and assuming the pixels are square, usually the spread_x and spread_y are the same value.

raster direct to 3d surface

If you can import the raster as a png or other supported file format, it is possible to convert straight from the image to a mesh..

import bpy
from math import floor

scalar_z = 2.3    #  amplification in z direction
scalar_xy = 1.0   #  amplification in xy


def generate_face_keys(num_x, num_y):
    faces = []
    verts_per_side_x = num_x

    for i in range((num_x - 1) * (num_y - 1)):
        x = (i % (num_x - 1))
        y = floor(i / (num_x - 1))

        level = x + (y * verts_per_side_x)
        idx1 = level
        idx2 = level + 1
        idx3 = level + verts_per_side_x + 1
        idx4 = level + verts_per_side_x
        faces.append([idx1, idx2, idx3, idx4])

    return faces


def geometry_function(image_name):
    img = bpy.data.images[image_name]
    num_x = w = img_width = img.size[0]
    num_y = h = img_height = img.size[1]

    # work on copy only
    pxls = img.pixels[:]
    num_pixels = int(len(pxls) / 4)

    verts = []
    add_vertex = verts.append

    # generator expression
    gen_obj = (i for i in pxls)

    for idx in range(num_pixels):
        y = int(idx / w) * scalar_xy
        x = (idx % w) * scalar_xy
        r = next(gen_obj)
        g = next(gen_obj)
        b = next(gen_obj)
        a = next(gen_obj)

        # height-map, would be grayscale, so z can be set to r*scalar
        z = r * scalar_z
        add_vertex((x, y, z))

    faces = generate_face_keys(num_x, num_y)
    return verts, faces


def generate(func, image_name, geometry_name):
    verts, faces = func(image_name)

    mesh = bpy.data.meshes.new(geometry_name + "_mesh")
    mesh.from_pydata(verts, [], faces)
    mesh.update()

    obj = bpy.data.objects.new(geometry_name, mesh)
    bpy.context.scene.objects.link(obj)
    obj.select = True
    obj.show_bounds = True


generate(geometry_function, "wvb_flat_2b.png", "3d_Surface_Object")

produces: