2
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ If nobody has an idea, I might try to just delete all polygons whose x,y are outside of the bounding polygon. This is of course less elegant than a clean crop, but could give a first approximation. $\endgroup$ – gelonida Dec 2 '19 at 19:35
  • $\begingroup$ ...it's quite possible.. just a couple of questions for a programmer that might come along.. Are the cropping polygon points ordered? Is it always a rectangle? If not, do you want the crop to follow the polygon's edges, or the edges of its bounding box? $\endgroup$ – Robin Betts Dec 2 '19 at 21:08
  • 1
    $\begingroup$ One approach might be to make repeated use of the bmesh bisect plane operator, once for each edge of the poly. I can't get to it right now, maybe someone else will... $\endgroup$ – Robin Betts Dec 2 '19 at 21:37
  • 1
    $\begingroup$ 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. Have to see if that is viable. $\endgroup$ – gelonida Dec 3 '19 at 0:55
  • 1
    $\begingroup$ still looking for an answer. $\endgroup$ – gelonida Dec 11 '19 at 23:47
1
$\begingroup$

Boolean intersection modifier

enter image description here

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

enter image description here

$\endgroup$
  • 1
    $\begingroup$ @RobinBetts bisect plane would be easy to implement as long as the poly is convex. $\endgroup$ – batFINGER Dec 12 '19 at 8:28
  • $\begingroup$ thanks a lot. Well this is what I tried, but something must be wrong either with my polygon or my grid. Also (probably not related to the issue, but who knows) how do you have the surface non transparent, but the polygon in wireframe? For me either both are with surfaces or both are displayed as wireframes. Will continue investigating when back from work $\endgroup$ – gelonida Dec 12 '19 at 8:50
  • $\begingroup$ Properties > Object Properties > Viewport Display $\endgroup$ – batFINGER Dec 12 '19 at 9:00
  • $\begingroup$ ...@batFINGER ..I was thinking some kind of point-in-polygon test first.. external algorithm, or maybe by using Blender's own BVH, and ray-casting down an orthogonal view? ..then bisecting .. then tagging new geometry if on edges of the clipping polygon. All a bit scrappy, though. $\endgroup$ – Robin Betts Dec 12 '19 at 11:10

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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