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

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
yOrigin = transform
pixelWidth = transform
pixelHeight = transform

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;
zc = data[0,0]
if (zc <= 0):
zc = 340
vertices.append((xc, yc, zc))

return vertices, faces

def generate(image_name, geometry_name):
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)
obj.select = True
obj.show_bounds = True

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

generate(theFile, "theModel")

• – user1853
Sep 18, 2015 at 5:16
• Do you perhaps have one such file to share, i might try another way Sep 22, 2015 at 18:31
• that's more than 2 million vertices, even on a fast computer that'll take a while to process :) but i wouldn't give up on python just yet.. Sep 22, 2015 at 18:54
• Here's a link that shares a smaller file: drive.google.com/file/d/0BwNOkTM_fO6-ZU9PRTRUOHdGaFU/… Sep 23, 2015 at 5:21
• is it fast for you to make a csv from the large geoTIFF? the less work python has to do the better,.. if so a csv with just the heights wouldn't hurt Sep 23, 2015 at 6:36

### 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
num_y = h = img_height = img.size

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

verts = []

# 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

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) 