Dividing terrain mesh up into multiple tiles with script - border gap issue 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.
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.