Dividing terrain mesh up into multiple tiles with script - border gap issue

Question

Having found a fantastic Script for dividing up terrain from user batFINGER

within this post - blender.stackexchange.com/questions/80460/slice-up-terrain-mesh-into-chunks/133258

import bpy, bmesh
from bpy import context
from  mathutils import Vector
# bounding box helper methods
def bbox(ob):
    return (Vector(b) for b in ob.bound_box)

def bbox_center(ob):
    return sum(bbox(ob), Vector()) / 8

def bbox_axes(ob):
    bb = list(bbox(ob))
    return tuple(bb[i] for i in (0, 4, 3, 1))

def slice(bm, start, end, segments):
    if segments == 1:
        return
    def geom(bm):
        return bm.verts[:] + bm.edges[:] + bm.faces[:]
    planes = [start.lerp(end, f / segments) for f in range(1, segments)]
    #p0 = start
    plane_no = (end - start).normalized() 
    while(planes): 
        p0 = planes.pop(0)                 
        ret = bmesh.ops.bisect_plane(bm, 
                geom=geom(bm),
                plane_co=p0, 
                plane_no=plane_no)
        bmesh.ops.split_edges(bm, 
                edges=[e for e in ret['geom_cut'] 
                if isinstance(e, bmesh.types.BMEdge)])


bm = bmesh.new()
ob = context.object
me = ob.data
bm.from_mesh(me)

o, x, y, z = bbox_axes(ob)        

x_segments = 2
y_segments = 2
z_segments = 1

slice(bm, o, x, x_segments)
slice(bm, o, y, y_segments)
slice(bm, o, z, z_segments)    
bm.to_mesh(me)

bpy.ops.object.mode_set(mode='EDIT')
bpy.ops.mesh.separate(type='LOOSE')
bpy.ops.object.mode_set()

I'm then trying to make use of it in combo with BlenderGIS

Cutting up the heightmap works well, but when I try and layer sat image tiles (already divided up within QGIS) over the top, I notice there is a tiny anomaly on all the edges - a gap to some extent. which it seems to be trying to fill with a different part of the image - which stands out.

When posting about it on BlenderGIS - github.com/domlysz/BlenderGIS/issues/172

It was noted that an additional rows + 1 and columns + 1 might be added so borders overlap.

Not being a Python wiz just yet, can anyone help me understand the tweaks required to adjust the above script to potentially handle the gap issue?

Update - 18/05/19

So one thing I also noticed was a loss of data when importing my heightmap (.asc) into blender, here is the header detail on asc -

When importing into blender I expected 1250 faces along one edge - but in fact it seemed to lose one, and only had 1249:

This might not be a huge issue, however when I then run the script to split into 4 x 4 - dividing 1249 by 4 seems to prompt the faces to be split into 3 sections of 313 x 313 (rounding up) but the edge faces are squashed:

Because I thought this might be contributing to the issue I regenerated my square heightmap with dimensions of 10020 mtrs which made following:

so losing 1 face when imported to blender gave me faces of 1252 which cleanly divides by 4, so now my edges look neat:

My next conundrum is that if I shift the objects following the dividing of heightmap using ob.location = (0, 0, 0) then I am unable to import jpg sat tiles using BlenderGIS as it returns error 'No Overlap Data' Clearly it want's the sat tiles to match with the heightmap's converted blender geo location (however it does that within BlenderGIS)?

and if I try and import the sat tiles at the heightmaps imported location: Vector((644784.5, 4979542.0, 0.0))

I still get the usual edge anomaly:

perhaps shifting all 16 tiles to the origin will resolve that - but I'm unsure how to move all of them in one batch move?

I'm also wondering if I should just try laying the sat images over the divided heightmap sections using default UV editor / texture config, rather than relying on BlenderGIS custom import feature, as maybe the Geo location element is confusing things? I think I tried one tile like that, but had other issues getting it to lay over terrain without the texture falling through the mesh.

Answer 1

Finally got this working, turns out you require an overlap on both DEM and Sat image

and I had to tile both in QGIS before bringing them into Blender

Full description of the resolution here: https://github.com/domlysz/BlenderGIS/issues/172

And summary of fix for issue:

The overlapping of 1 pixel for DEM tiles help to compensate the loose of last column and line when connecting only pixel's center, this trick allows to obtain perfectly joint tiles.

The overlapping for satellite tiles help to compensate the difference of tiles size between the 2 rasters datasource because even if the 2 rasters have the same extent, as long as they don't have the same pixel size in meters then dividing the number of pixel by 4 will not produce a cut at the same coordinate.

Figures for my assets:

DEM 1253x1253 >> retile 315x315 with 1px overlap SAT 10020x10020 >> retile 2600x2600 with 95px overlap

I used GDAL process 'Retile' in QGIS to divide tiles + overlap.

Answer 2

Transpose to origin.

I have a feeling your issue is the location of the map, shown in blender units

>>> C.object.location
Vector((644784.5, 4979532.0, 0.0))

which is a bloody long way from the origin, and could very likely be producing rounding errors. Suggest transposing such that one tile corner is at origin.

With all tiles selected, and bottom left hand corner active run this code in python console.

>>> v = C.object.location.copy()
>>> for o in C.selected_objects:
...     o.location -= v

or in the script in question,

ob = context.object
ob.location = (0, 0, 0) # <== insert this line
me = ob.data

to set the location of dem to origin before splitting.

Note: This can also be set within BlenderGIS, read the docs https://github.com/domlysz/BlenderGIS/wiki/Georeferencing-management#custom-properties and set the scene origin geolocation accordingly. Theoretically this is best solution as it will move the map to origin, and still import tiles correctly

Image from blenderGIS wiki

there is no need to emulate global coordinates in meters in blender. Huge numbers can lead to errors.

Digits of precision of vertices coordinates

Scaling X Y Z precision?

Should I model my animation assets to real-world scale?

Edit: Lastly re question edit, script above does not "squash" a face, it cuts it with a new bisecting edge. Making the grid divisible, by in this case 4, simply makes the cut at an existing edge.