1""" 2Generate DTM in Blender 3 4Using the Displace and Decimate modifiers to generate a terrain 3D model from 5(geo)images (png, tif, geotif, etc) 6 7usage: blender -b -P blend_dtm.py -- dsm.tif image.jpg [terrain_resolution [strength]] 8 9The DSM must be Float32, Blender does not support Float64, If you get the error: 10 imb_loadtiff: Failed to read tiff image. 11 IMB_ibImageFromMemory: unknown fileformat 12Convert first using: gdal_translate -ot Float32 dsm.tif dsm32.tif 13 14You can safely ignore the "unknown field with tag" TIFFReadDirectory Warning. 15""" 16import os 17import sys 18import json 19import subprocess 20import bpy # Blender Python API 21 22TERRAIN_RESOLUTION=0.50 # in meter 23 24def new_texture(): 25 bpy.ops.texture.new() 26 return bpy.data.textures[-1] 27 28def new_material(): 29 bpy.ops.material.new() 30 return bpy.data.materials[-1] 31 32def open_image(filepath): 33 bpy.ops.image.open(filepath=filepath) 34 return bpy.data.images[-1] 35 36def new_image_texture(material, image, name='img'): 37 """ Create a new texture from an image (which we open) 38 39 we need a material here, to set the texture active 40 otherwise we can't change its type (to image in this case) 41 """ 42 material.active_texture = new_texture() 43 material.active_texture.name = name 44 material.active_texture.type = 'IMAGE' 45 material.active_texture.image = image 46 return material.active_texture 47 48def new_grid(name='Grid', x_subdivisions=10, y_subdivisions=10, radius=1, scale=(1,1,1)): 49 bpy.ops.object.select_all(action='DESELECT') 50 # if blender version < 2.68 : radius -> size 51 if bpy.app.version < (2,68): 52 bpy.ops.mesh.primitive_grid_add(x_subdivisions=x_subdivisions, \ 53 y_subdivisions=y_subdivisions, size=radius) 54 else: 55 bpy.ops.mesh.primitive_grid_add(x_subdivisions=x_subdivisions, \ 56 y_subdivisions=y_subdivisions, radius=radius) 57 bpy.context.object.scale = scale 58 bpy.context.object.name = name 59 return bpy.context.object 60 61def select_only(obj): 62 bpy.ops.object.select_all(action='DESELECT') 63 obj.select = True 64 bpy.context.scene.objects.active = obj 65 66def subdivide(obj, number_cuts=1): 67 select_only(obj) 68 bpy.ops.object.mode_set(mode='EDIT') 69 bpy.ops.mesh.subdivide(number_cuts=number_cuts) 70 bpy.ops.object.mode_set(mode='OBJECT') 71 return obj 72 73def new_cube(name='Cube', scale=(1,1,1)): 74 bpy.ops.object.select_all(action='DESELECT') 75 bpy.ops.mesh.primitive_cube_add() 76 bpy.context.object.scale = scale 77 bpy.context.object.name = name 78 return bpy.context.object 79 80def set_viewport(viewport_shade='WIREFRAME', clip_end=1000): 81 """ Set the default view mode 82 83 :param viewport_shade: enum in ['BOUNDBOX', 'WIREFRAME', 'SOLID', 'TEXTURED'], default 'WIREFRAME' 84 """ 85 for area in bpy.context.window.screen.areas: 86 if area.type == 'VIEW_3D': 87 for space in area.spaces: 88 if space.type == 'VIEW_3D': 89 space.viewport_shade = viewport_shade 90 space.clip_end = clip_end 91 space.grid_scale = 10 92 space.grid_lines = 50 93 94def save(filepath=None, check_existing=False, compress=True): 95 """ Save .blend file 96 97 :param filepath: File Path 98 :type filepath: string, (optional, default: current file) 99 :param check_existing: Check and warn on overwriting existing files 100 :type check_existing: boolean, (optional, default: False) 101 :param compress: Compress, Write compressed .blend file 102 :type compress: boolean, (optional, default: True) 103 """ 104 if not filepath: 105 filepath = bpy.data.filepath 106 bpy.ops.wm.save_mainfile(filepath=filepath, check_existing=check_existing, 107 compress=compress) 108 109def setup(): 110 """ Setup the scene """ 111 # Delete the default objects 112 bpy.ops.object.select_all(action='SELECT') 113 bpy.ops.object.delete() 114 # Add light (sun) 115 bpy.ops.object.lamp_add(type='SUN', location=(0, 0, 300)) 116 # Set Game mode 117 bpy.context.scene.render.engine = 'BLENDER_GAME' 118 # make sure OpenGL shading language shaders (GLSL) is the 119 # material mode to use for rendering 120 bpy.context.scene.game_settings.material_mode = 'GLSL' 121 # Set the unit system to use for button display (in edit mode) to metric 122 bpy.context.scene.unit_settings.system = 'METRIC' 123 # Select the type of Framing to Extend, 124 # Show the entire viewport in the display window, 125 # viewing more horizontally or vertically. 126 bpy.context.scene.game_settings.frame_type = 'EXTEND' 127 # see the texture in realtime (in the 3D viewport) 128 set_viewport('TEXTURED', 10000) 129 # Set the color at the horizon to dark azure 130 bpy.context.scene.world.horizon_color = (0.05, 0.22, 0.4) 131 # Display framerate and profile information of the simulation 132 bpy.context.scene.game_settings.show_framerate_profile = True 133 134def usage(): 135 sys.stderr.write(__doc__) 136 137def ext_exec(cmd, python=None): 138 if not python: 139 python = 'python%i.%i'%(sys.version_info.major, sys.version_info.minor) 140 return subprocess.getoutput('%s -c"%s"' % (python, cmd) ) 141 142def fix_python_path(python=None): 143 pythonpath = ext_exec("import os,sys;print(os.pathsep.join(sys.path))") 144 sys.path.extend(pythonpath.split(os.pathsep)) 145 146def get_gdalinfo(filepath): 147 """ Get the GeoTransform from the :param:filepath DEM using gdal 148 149 In a north up image, padfTransform[1] is the pixel width, 150 and padfTransform[5] is the pixel height. The upper left corner of the 151 upper left pixel is at position (padfTransform[0],padfTransform[3]). 152 153 The default transform is (0,1,0,0,0,1) 154 """ 155 # Run externally since `gdal.Open` crashes Blender :-/ 156 meta = ext_exec("import json,gdal,gdalconst;" 157 "g=gdal.Open('%s',gdalconst.GA_ReadOnly);" 158 "print(json.dumps({'transform':g.GetGeoTransform()," 159 "'minmax':g.GetRasterBand(1).ComputeRasterMinMax()," 160 "'meta':g.GetMetadata()}))"%filepath, 'python2') 161 return json.loads(meta) # {'transform':(0,1,0,0,0,1)} 162 163def main(argv=[]): 164 args = [] 165 if '--' in argv: 166 args = argv[argv.index('--')+1:] 167 168 if len(args) < 2: 169 usage() 170 return 1 171 172 strength = 1 173 if len(args) > 2: 174 terrain_resolution = float(args[2]) 175 if len(args) > 3: 176 strength = float(args[3]) 177 else: 178 terrain_resolution = TERRAIN_RESOLUTION 179 180 image_dem = open_image(args[0]) 181 image_img = open_image(args[1]) 182 if image_dem.size[:] != image_img.size[:]: 183 print('[WARN] The size of the 2 image differ') 184 185 fix_python_path() 186 gdalinfo = get_gdalinfo(args[0]) 187 geot = gdalinfo['transform'] 188 meta = gdalinfo['meta'] 189 xsize = image_dem.size[0] * abs(geot[1]) # in meters 190 ysize = image_dem.size[1] * abs(geot[5]) # in meters 191 192 translation = [0.0, 0.0, 0.0] 193 if 'CUSTOM_X_ORIGIN' in meta and 'CUSTOM_Y_ORIGIN' in meta: 194 custom_x_origin = float(meta['CUSTOM_X_ORIGIN']) 195 custom_y_origin = float(meta['CUSTOM_Y_ORIGIN']) 196 print("[gdal] got custom (%f, %f)"%(custom_x_origin, custom_y_origin)) 197 center_x_utm = geot[0] + image_dem.size[0] * geot[1] / 2 198 center_y_utm = geot[3] + image_dem.size[1] * geot[5] / 2 199 translation = [ center_x_utm - custom_x_origin, 200 center_y_utm - custom_y_origin, 0.0 ] 201 if gdalinfo['minmax']: 202 translation[2] = - round(gdalinfo['minmax'][0] + 1) 203 setup() 204 205 ######################################################################### 206 # Add our ground object 207 ground = new_grid('Ground', xsize/terrain_resolution, \ 208 ysize/terrain_resolution, 1, (xsize/2.0, ysize/2.0, 1)) 209 210 ######################################################################### 211 # Add material 212 ground.active_material = new_material() 213 ground.active_material.specular_intensity = 0 214 ground.active_material.name = 'Ground' 215 material = ground.active_material 216 217 ######################################################################### 218 # Add texture 219 texture_img = new_image_texture(material, image_img, 'img') 220 221 material.active_texture_index += 1 222 223 if bpy.app.version > (2,63): 224 # since 2.64, default colorspace is sRGB which breaks decimate modifier 225 image_dem.colorspace_settings.name = 'Linear' 226 texture_dem = new_image_texture(material, image_dem, 'dem') 227 # do not show on the material (use only later for the terrain modifier) 228 material.use_textures[material.active_texture_index] = False 229 230 #bpy.ops.object.shade_smooth() 231 # displace from image_dem (gdal) 232 add_displace_modifier(ground, texture_dem, strength=strength, apply=True) 233 # unlink dem after apply (reduce size) 234 material.texture_slots.clear(1) 235 image_dem.user_clear() 236 bpy.data.images.remove(image_dem) 237 texture_dem.user_clear() 238 bpy.data.textures.remove(texture_dem) 239 #add_decimate_modifier(ground) 240 #add_water_cube((xsize+1, ysize+1, max(image_dem.pixels)/2)) 241 242 # move to: custom origin - center (utm from gdalinfo) 243 ground.location = translation 244 245 bpy.ops.file.pack_all() # bpy.ops.image.pack(as_png=True) 246 save('/tmp/dtm.blend') 247 # bpy.ops.view3d.view_selected(use_all_regions=True) # XXX no poll() 248 print("\n----\nsaved in /tmp/dtm.blend\n----\n") 249 return 0 250 251def add_water_cube(scale=(1,1,1)): 252 """ Add cube 'water' to hide the bottom 253 254 (save memory, don't render many vertices) 255 """ 256 cube = new_cube('Water', scale) 257 cube.active_material = new_material() 258 cube.active_material.name = 'Water' 259 cube.active_material.use_shadeless = True 260 cube.active_material.diffuse_color = (0.3, 0.4, 0.8) 261 return cube 262 263def add_displace_modifier(obj, texture, strength=1, direction='Z', apply=False): 264 """ Add displace modifier 265 266 http://wiki.blender.org/index.php/Doc:2.6/Manual/Modifiers/Deform/Displace 267 """ 268 select_only(obj) 269 bpy.ops.object.modifier_add(type='DISPLACE') 270 modifier = bpy.context.object.modifiers[-1] 271 modifier.texture = texture 272 modifier.strength = strength 273 modifier.direction = direction 274 if apply: 275 bpy.ops.object.modifier_apply(modifier=modifier.name) 276 return modifier 277 278def add_decimate_modifier(obj, ratio=0.5, apply=False, show_viewport=False): 279 """ Add displace modifier 280 281 http://wiki.blender.org/index.php/Doc:2.6/Manual/Modifiers/Generate/Decimate 282 """ 283 select_only(obj) 284 bpy.ops.object.modifier_add(type='DECIMATE') 285 modifier = bpy.context.object.modifiers[-1] 286 modifier.show_viewport = show_viewport # no preview : faster 287 modifier.ratio = ratio 288 if apply: 289 bpy.ops.object.modifier_apply(modifier=modifier.name) 290 return modifier 291 292def add_hillshade(demfile, material): 293 """ Add hillshade """ 294 hillshade_png = '%s.hillshade.png'%demfile 295 if not os.path.isfile(hillshade_png): 296 subprocess.getstatusoutput("gdaldem hillshade %s %s -of png"%(demfile, hillshade_png)) 297 image_hsh = open_image(hillshade_png) 298 material.active_texture_index += 1 299 texture_hsh = new_image_texture(material, image_hsh, 'hsh') 300 texture_hsh.use_normal_map = True 301 texture_slot_hsh = material.texture_slots[material.active_texture_index] 302 texture_slot_hsh.use_map_color_diffuse = False 303 texture_slot_hsh.use_map_normal = True 304 texture_slot_hsh.normal_factor = -0.5 305 306if __name__ == '__main__': 307 main(sys.argv) 308