How to clip a grid of height points within a polygon base shape?

Question

I'm very new to blender, but have some scripting background.

I have some height measurements for a surface. (z values for x,y coordinates)

I was able to use Python's scipy.interpolate.griddata to convert this into a grid of Z values. ( for lazy ones here the generated grid data https://pastebin.com/4xtuJgac ) With Create model from XYZ data points I as able to import this grid into blender.

However now I'd like to clip off everything, that is not within the x, y range of a specified polygon.

The scripts below show what I've done so far.

The first part I execute to generate my grid data with a small built in data set in a standard python executable.

from scipy.interpolate import griddata
import numpy as np

values = [ # sample data for Stack Exchange
    (2.32702550990507, 82.75699786003679, 27.789999999999992),
    (11.332438979996368, 14.02555966982618, 12.860000000000014),
    (0.7804978101048619, 75.70470666978508, 26.610000000000014),
    (2.2677287301048636, 75.15228118980303, 28.039999999999992),
    (18.682603270048276, 16.836770640220493, 12.919999999999987),
    (12.09464186988771, 17.609230670146644, 13.289999999999992),
    (5.603262080112472, 47.57122483011335, 21.53),
    (26.761874669929966, 16.459352959878743, 12.400000000000006),
    (9.299435199936852, 13.942597850225866, 12.580000000000013),
    (10.332293499959633, 14.73426689999178, 12.740000000000009),
    (5.753992310026661, 40.484739230014384, 18.909999999999997),
    (4.8296503101009876, 44.45012723002583, 20.22),
    (2.676168560050428, 76.70651125023142, 27.72),
    (12.363012589979917, 0.5158127900213003, 10.300000000000011),
    (6.900217070011422, 48.241429650224745, 22.24000000000001),
    ]

px = np.array([v[0] for v in values])
py = np.array([v[1] for v in values])
pz = np.array([v[2] for v in values])
x = np.linspace(int(min(px)), int(max(px))+1, 64)
y = np.linspace(int(min(py)), int(max(py))+1, 64)
X, Y = np.meshgrid(x,y)

gdata = griddata((px, py), pz, (X, Y), method="nearest")

with open("gridpoints.csv", "w") as fout:
    for ix, row in enumerate(gdata):
        for iy, val in enumerate(row):
            fout.write("%f, %f, %f\n" % (ix, iy, val))

Here the blender code to import this grid data. ( https://pastebin.com/4xtuJgac if you don't want to create them with above script)

import os
import csv
import bpy
from operator import itemgetter

MYDIR = os.path.realpath(os.path.dirname(__file__))
csvfile = open(os.path.join(MYDIR, "gridpoints.csv"))
inFile = csv.reader(csvfile, delimiter=',', quotechar='"')

# Read and sort the vertices coordinates (sort by x and y)
vertices = sorted( [(float(r[0]), float(r[1]), float(r[2])) for r in inFile], key = itemgetter(0,1) )

# ********* Assuming we have a rectangular grid *************
xSize = next(i for i in range(len(vertices)) if vertices[i][0] != vertices[i+1][0]) + 1  # Find the first change in X
ySize = len(vertices) // xSize

polygons = [(i, i - 1, i - 1 + xSize, i + xSize) for i in range( 1, len(vertices) - xSize ) if i % xSize != 0]

# Generate the polygons (four vertices linked)
name = "grid_unclipped"
mesh = bpy.data.meshes.new(name)  # Create the mesh (inner data)
obj = bpy.data.objects.new(name, mesh)  # Create an object

obj.data.from_pydata(vertices, [], polygons)  # Associate vertices and polygons

# obj.scale = (1, 5, 0.2) #Scale it (if needed)
for p in obj.data.polygons:  # Set smooth shading (if needed)
    p.use_smooth = True

bpy.context.scene.objects.link(obj)  # Link the object to the scene
# or for 2.8
# bpy.context.scene.collection.objects.link(obj)

and here the x,y coordinates of the polygon section, that I'd like to use for cropping.

section = [
    ( 9.5500, 14.3610),
    (38.1858, 17.4571),
    (35.2338, 44.5318),
    ( 5.7485, 41.6682),
    ]

# This is where I don't know how to proceed

Addendum 2019012012: . Meanwhile I created the polygon within blender with following code.

    name2 = "bounding_poly"
    mesh2 = bpy.data.meshes.new(name2)  # Create the mesh (inner data)
    obj2 = bpy.data.objects.new(name2, mesh2)  # Create an object


    vertices = list((x, y, 0) for (x, y) in section)
    edges = []
    for idx in range(len(vertices) - 1):
        edges.append((idx, idx+1))
    edges.append((len(vertices) - 1, 0))
    polygons = [tuple(idx for idx in range(len(vertices)))]

    obj2.data.from_pydata(vertices, edges, polygons)  # Associate vertices and polygons

    bpy.context.scene.collection.objects.link(obj2)  # Link the object to the 

Afterwards I enter edit mode and extrude the polygon to create some kind of prism (almost cube shaped)

But now I still don't know how I could do the clipping.

Answer 1

Boolean intersection modifier

another thing (but have to see how to do this as script) mgiht be to create a mesh (one flat surface of the polygon), then protrude it high and low enough in z direction and intersect my grid with this polygon prism

Yes that would be my take on it. Create a prism then use the boolean intersect modifier. Example picture shows result after running on simple z displaced grid dem.

The prism is an extruded ngon.

Adding and applying modifier in the python console, grid object has context.

>>> o = C.object
>>> bm = o.modifiers.new("BM", 'BOOLEAN')
>>> bm.object = C.scene.objects.get("Circle") # the prism
>>> bm.operation = 'INTERSECT'

>>> bpy.ops.object.modifier_apply(
modifier_apply()
bpy.ops.object.modifier_apply(apply_as='DATA', modifier="")
Apply modifier and remove from the stack
>>> bpy.ops.object.modifier_apply(modifier=bm.name)
Info: Applied modifier was not first, result may not be as expected
{'FINISHED'}

Note first modifer re warning is a displacement modifier to make DEM and can be ignored.

Script to product the section. Set the z and the solidify thickness based on max / min range of z from input grid data. simply hardcoded as 40 per example

import bpy
context = bpy.context
section = [
    ( 9.5500, 14.3610),
    (38.1858, 17.4571),
    (35.2338, 44.5318),
    ( 5.7485, 41.6682),
    ]

verts = [(x, y, 40) for x, y in section]
faces = [list(range(len(section)))]

me = bpy.data.meshes.new("Section")
me.from_pydata(verts, [], faces)

ob = bpy.data.objects.new("Section", me)
ob.display_type = 'WIRE'
context.collection.objects.link(ob)
sol = ob.modifiers.new("Sol", 'SOLIDIFY')
sol.thickness = 40

And can confirm this works as described on input data