Files
osmWorkflow/blender/generate_nantaizi.py
2026-07-24 16:03:15 +08:00

1193 lines
48 KiB
Python

"""Build a lightweight 3D park model from the Nantaizi OSM export.
Run from Blender 4.x:
blender --python generate_nantaizi.py -- \
--osm "/Users/que01/Desktop/南台子湖创新谷OSM.osm"
The OSM bounds element is used deliberately. The export contains subway
relation members far outside the park, so using every node for the scene
extent would produce a misleadingly large model.
"""
import json
import math
import os
import sys
import xml.etree.ElementTree as ET
from collections import defaultdict
import bpy
from mathutils import Matrix, Vector
DEFAULT_OSM = "/Users/que01/Desktop/南台子湖创新谷OSM.osm"
DEFAULT_GEOJSON = (
"/Users/que01/osm2streets-qgis-workflow/outputs/"
"nantaizi-lake-innovation-valley/osm2streets_web_out"
)
DEFAULT_OUTPUT = (
"/Users/que01/osm2streets-qgis-workflow/outputs/"
"nantaizi-lake-innovation-valley/nantaizi_lake_innovation_valley.blend"
)
DEFAULT_RENDER = (
"/Users/que01/osm2streets-qgis-workflow/outputs/"
"nantaizi-lake-innovation-valley/nantaizi_lake_innovation_valley.png"
)
TEXTURE_ROOT = os.path.abspath(os.path.join(
os.path.dirname(__file__), "..", "assets", "textures", "polyhaven"
))
MODEL_ROOT = os.path.abspath(os.path.join(
os.path.dirname(__file__), "..", "assets", "models", "polyhaven"
))
TREE_MODEL_PATH = os.path.join(
MODEL_ROOT, "78-hazelnutbush", "Hazelnut.obj"
)
# These two footprints are ordinary office buildings despite their current
# OSM building=industrial tags. Keep the correction explicit and traceable.
OFFICE_OVERRIDE_WAY_IDS = {"117753521", "117753535"}
def cli_args():
values = {"osm": DEFAULT_OSM, "geojson": DEFAULT_GEOJSON,
"output": DEFAULT_OUTPUT, "render": DEFAULT_RENDER}
argv = sys.argv[sys.argv.index("--") + 1:] if "--" in sys.argv else []
i = 0
while i < len(argv):
if argv[i].startswith("--") and i + 1 < len(argv):
values[argv[i][2:]] = argv[i + 1]
i += 2
else:
i += 1
return values
def tags(element):
return {t.attrib.get("k", ""): t.attrib.get("v", "")
for t in element.findall("tag")}
def parse_osm(path):
root = ET.parse(path).getroot()
bounds_node = root.find("bounds")
if bounds_node is None:
raise RuntimeError("OSM file does not contain a bounds element")
bounds = {"min_lon": float(bounds_node.attrib["minlon"]),
"min_lat": float(bounds_node.attrib["minlat"]),
"max_lon": float(bounds_node.attrib["maxlon"]),
"max_lat": float(bounds_node.attrib["maxlat"])}
nodes = {}
point_features = []
for node in root.findall("node"):
try:
node_id = int(node.attrib["id"])
coord = (float(node.attrib["lon"]), float(node.attrib["lat"]))
node_tags = tags(node)
nodes[node_id] = coord
if node_tags:
point_features.append({"id": node.attrib.get("id", ""),
"coord": coord, "tags": node_tags})
except (KeyError, ValueError):
continue
ways = []
for way in root.findall("way"):
if way.attrib.get("action") == "delete":
continue
refs = []
for ref in way.findall("nd"):
try:
refs.append(int(ref.attrib["ref"]))
except (KeyError, ValueError):
pass
coords = [nodes[r] for r in refs if r in nodes]
if len(coords) >= 2:
ways.append({"id": way.attrib.get("id", ""),
"coords": coords, "tags": tags(way)})
return bounds, ways, point_features
class Projector:
def __init__(self, bounds):
self.bounds = bounds
self.lon0 = (bounds["min_lon"] + bounds["max_lon"]) / 2
self.lat0 = (bounds["min_lat"] + bounds["max_lat"]) / 2
self.m_per_lat = 111320.0
self.m_per_lon = 111320.0 * math.cos(math.radians(self.lat0))
def xy(self, lon_lat):
lon, lat = lon_lat
return ((lon - self.lon0) * self.m_per_lon,
(lat - self.lat0) * self.m_per_lat)
def inside(self, lon_lat, pad=0.00035):
lon, lat = lon_lat
b = self.bounds
return (b["min_lon"] - pad <= lon <= b["max_lon"] + pad and
b["min_lat"] - pad <= lat <= b["max_lat"] + pad)
def ring(self, coords):
return [self.xy(c) for c in coords]
def new_collection(name):
collection = bpy.data.collections.new(name)
bpy.context.scene.collection.children.link(collection)
return collection
def make_material(name, color, roughness=0.8, metallic=0.0):
material = bpy.data.materials.get(name) or bpy.data.materials.new(name)
material.diffuse_color = (*color, 1.0)
material.use_nodes = True
bsdf = material.node_tree.nodes.get("Principled BSDF")
if bsdf:
bsdf.inputs["Base Color"].default_value = (*color, 1.0)
bsdf.inputs["Roughness"].default_value = roughness
bsdf.inputs["Metallic"].default_value = metallic
return material
def add_procedural_surface(material, colors, scale=2.0, detail=2.0, bump_strength=0.08):
"""Add small-scale color and normal variation without external textures."""
nodes = material.node_tree.nodes
links = material.node_tree.links
bsdf = nodes.get("Principled BSDF")
if not bsdf:
return
noise = nodes.new("ShaderNodeTexNoise")
noise.inputs["Scale"].default_value = scale
noise.inputs["Detail"].default_value = detail
noise.inputs["Roughness"].default_value = 0.65
texcoord = nodes.new("ShaderNodeTexCoord")
ramp = nodes.new("ShaderNodeValToRGB")
ramp.color_ramp.elements[0].color = (*colors[0], 1.0)
ramp.color_ramp.elements[1].color = (*colors[1], 1.0)
bump = nodes.new("ShaderNodeBump")
bump.inputs["Strength"].default_value = bump_strength
bump.inputs["Distance"].default_value = 0.12
links.new(texcoord.outputs["Generated"], noise.inputs["Vector"])
links.new(noise.outputs["Fac"], ramp.inputs["Fac"])
links.new(ramp.outputs["Color"], bsdf.inputs["Base Color"])
links.new(noise.outputs["Fac"], bump.inputs["Height"])
links.new(bump.outputs["Normal"], bsdf.inputs["Normal"])
def make_textured_material(name, diffuse_file, normal_file, roughness,
scale, normal_is_bump=False, metallic=0.0,
tint=None, tint_factor=0.0):
diffuse_path = os.path.join(TEXTURE_ROOT, diffuse_file)
normal_path = os.path.join(TEXTURE_ROOT, normal_file)
if not os.path.exists(diffuse_path) or not os.path.exists(normal_path):
return make_material(name, (0.5, 0.5, 0.5), roughness, metallic)
material = make_material(name, (0.5, 0.5, 0.5), roughness, metallic)
nodes = material.node_tree.nodes
links = material.node_tree.links
bsdf = nodes.get("Principled BSDF")
texcoord = nodes.new("ShaderNodeTexCoord")
mapping = nodes.new("ShaderNodeMapping")
mapping.inputs["Scale"].default_value = (scale, scale, scale)
diffuse = nodes.new("ShaderNodeTexImage")
diffuse.image = bpy.data.images.load(diffuse_path, check_existing=True)
diffuse.extension = "REPEAT"
normal = nodes.new("ShaderNodeTexImage")
normal.image = bpy.data.images.load(normal_path, check_existing=True)
normal.image.colorspace_settings.name = "Non-Color"
normal.extension = "REPEAT"
links.new(texcoord.outputs["Generated"], mapping.inputs["Vector"])
links.new(mapping.outputs["Vector"], diffuse.inputs["Vector"])
links.new(mapping.outputs["Vector"], normal.inputs["Vector"])
if tint and tint_factor > 0.0:
tint_node = nodes.new("ShaderNodeRGB")
tint_node.outputs["Color"].default_value = (*tint, 1.0)
mix = nodes.new("ShaderNodeMixRGB")
mix.blend_type = "MIX"
mix.inputs["Fac"].default_value = tint_factor
links.new(diffuse.outputs["Color"], mix.inputs[1])
links.new(tint_node.outputs["Color"], mix.inputs[2])
links.new(mix.outputs["Color"], bsdf.inputs["Base Color"])
else:
links.new(diffuse.outputs["Color"], bsdf.inputs["Base Color"])
if normal_is_bump:
bump = nodes.new("ShaderNodeBump")
bump.inputs["Strength"].default_value = 0.22
bump.inputs["Distance"].default_value = 0.12
links.new(normal.outputs["Color"], bump.inputs["Height"])
links.new(bump.outputs["Normal"], bsdf.inputs["Normal"])
else:
normal_map = nodes.new("ShaderNodeNormalMap")
normal_map.inputs["Strength"].default_value = 0.52
links.new(normal.outputs["Color"], normal_map.inputs["Color"])
links.new(normal_map.outputs["Normal"], bsdf.inputs["Normal"])
return material
def make_image_material(name, image_path, roughness=0.8, metallic=0.0,
alpha_path=None, alpha_clip=0.33,
cull_backface=False, tint=None, tint_factor=0.0,
saturation=1.0, value=1.0, emission_color=None,
emission_strength=0.0):
material = make_material(name, (0.5, 0.5, 0.5), roughness, metallic)
nodes = material.node_tree.nodes
links = material.node_tree.links
bsdf = nodes.get("Principled BSDF")
if not bsdf or not os.path.exists(image_path):
return material
image = nodes.new("ShaderNodeTexImage")
image.image = bpy.data.images.load(image_path, check_existing=True)
color_output = image.outputs["Color"]
if saturation != 1.0 or value != 1.0:
hsv = nodes.new("ShaderNodeHueSaturation")
hsv.inputs["Saturation"].default_value = saturation
hsv.inputs["Value"].default_value = value
links.new(color_output, hsv.inputs["Color"])
color_output = hsv.outputs["Color"]
if tint and tint_factor > 0.0:
tint_node = nodes.new("ShaderNodeRGB")
tint_node.outputs["Color"].default_value = (*tint, 1.0)
mix = nodes.new("ShaderNodeMixRGB")
mix.blend_type = "MIX"
mix.inputs["Fac"].default_value = tint_factor
links.new(color_output, mix.inputs[1])
links.new(tint_node.outputs["Color"], mix.inputs[2])
color_output = mix.outputs["Color"]
links.new(color_output, bsdf.inputs["Base Color"])
if "Specular IOR Level" in bsdf.inputs:
bsdf.inputs["Specular IOR Level"].default_value = 0.2
elif "Specular" in bsdf.inputs:
bsdf.inputs["Specular"].default_value = 0.2
if emission_color and emission_strength > 0.0:
if "Emission Color" in bsdf.inputs:
bsdf.inputs["Emission Color"].default_value = (*emission_color, 1.0)
elif "Emission" in bsdf.inputs:
bsdf.inputs["Emission"].default_value = (*emission_color, 1.0)
if "Emission Strength" in bsdf.inputs:
bsdf.inputs["Emission Strength"].default_value = emission_strength
alpha_source = image
if alpha_path and os.path.exists(alpha_path):
alpha_source = nodes.new("ShaderNodeTexImage")
alpha_source.image = bpy.data.images.load(alpha_path, check_existing=True)
alpha_source.image.colorspace_settings.name = "Non-Color"
if "Alpha" in bsdf.inputs:
links.new(alpha_source.outputs["Alpha"], bsdf.inputs["Alpha"])
if hasattr(material, "blend_method"):
material.blend_method = "CLIP" if alpha_path else "OPAQUE"
if hasattr(material, "shadow_method"):
material.shadow_method = "CLIP" if alpha_path else "OPAQUE"
if hasattr(material, "alpha_threshold"):
material.alpha_threshold = alpha_clip
if hasattr(material, "surface_render_method") and alpha_path:
material.surface_render_method = "DITHERED"
if hasattr(material, "use_backface_culling"):
material.use_backface_culling = cull_backface
return material
def import_model(filepath):
extension = os.path.splitext(filepath)[1].lower()
if extension in {".gltf", ".glb"}:
bpy.ops.import_scene.gltf(filepath=filepath)
return
if extension == ".obj":
if hasattr(bpy.ops.wm, "obj_import"):
bpy.ops.wm.obj_import(filepath=filepath)
return
bpy.ops.import_scene.obj(filepath=filepath)
return
raise RuntimeError("Unsupported tree model format: " + extension)
def assign_tree_model_materials(objects):
model_dir = os.path.dirname(TREE_MODEL_PATH)
bark = make_image_material(
"Hazelnut Bark",
os.path.join(model_dir, "HazelnutBark.png"),
roughness=0.86,
)
leaves = make_image_material(
"Hazelnut Leaves",
os.path.join(model_dir, "HazelnutLeaves.png"),
roughness=0.82,
alpha_path=os.path.join(model_dir, "HazelnutLeavesMask.png"),
alpha_clip=0.18,
cull_backface=False,
tint=(0.25, 0.52, 0.18),
tint_factor=0.62,
saturation=1.45,
value=1.42,
emission_color=(0.21, 0.36, 0.15),
emission_strength=0.26,
)
for obj in objects:
if obj.type != "MESH":
continue
obj.data.materials.clear()
lowered = obj.name.lower()
if "leaf" in lowered:
obj.data.materials.append(leaves)
else:
obj.data.materials.append(bark)
for polygon in obj.data.polygons:
polygon.material_index = 0
polygon.use_smooth = True
class MeshBatch:
def __init__(self, name, collection, material):
self.name = name
self.collection = collection
self.material = material
self.vertices = []
self.faces = []
def add_polygon(self, ring, z):
if len(ring) < 3:
return
if ring[0] == ring[-1]:
ring = ring[:-1]
if len(ring) < 3:
return
start = len(self.vertices)
self.vertices.extend((x, y, z) for x, y in ring)
self.faces.append(tuple(range(start, start + len(ring))))
def add_prism(self, ring, base, height):
if len(ring) < 3:
return
if ring[0] == ring[-1]:
ring = ring[:-1]
if len(ring) < 3:
return
start = len(self.vertices)
self.vertices.extend((x, y, base) for x, y in ring)
self.vertices.extend((x, y, base + height) for x, y in ring)
n = len(ring)
self.faces.append(tuple(range(start, start + n)))
self.faces.append(tuple(range(start + n, start + 2 * n)))
for i in range(n):
j = (i + 1) % n
self.faces.append((start + i, start + j, start + n + j, start + n + i))
def finish(self):
if not self.vertices:
return None
mesh = bpy.data.meshes.new(self.name + "Mesh")
mesh.from_pydata(self.vertices, [], self.faces)
mesh.materials.append(self.material)
if self.name.startswith("Tree_") or self.name.startswith("Scrub_"):
for polygon in mesh.polygons:
polygon.use_smooth = True
mesh.update()
obj = bpy.data.objects.new(self.name, mesh)
self.collection.objects.link(obj)
return obj
def make_prism(name, ring, base, height, material, collection):
batch = MeshBatch(name, collection, material)
batch.add_prism(ring, base, height)
return batch.finish()
def add_roof(name, ring, z, material, collection):
batch = MeshBatch(name + "_Roof", collection, material)
batch.add_polygon(ring, z)
return batch.finish()
def add_wall_panel(batch, start, end, base, height, thickness=0.045, inset=0.08):
dx, dy = end[0] - start[0], end[1] - start[1]
length = math.hypot(dx, dy)
if length < 3.0:
return
ux, uy = dx / length, dy / length
a = (start[0] + dx * inset, start[1] + dy * inset)
b = (end[0] - dx * inset, end[1] - dy * inset)
nx, ny = -uy * thickness / 2, ux * thickness / 2
panel = [(a[0] + nx, a[1] + ny), (b[0] + nx, b[1] + ny),
(b[0] - nx, b[1] - ny), (a[0] - nx, a[1] - ny)]
batch.add_prism(panel, base, height)
def add_building_details(name, ring, height, industrial, materials, collection):
"""Add restrained facade and roof detail without changing OSM massing."""
footprint = ring[:-1] if len(ring) > 1 and ring[0] == ring[-1] else ring
if len(footprint) < 3:
return
glass_mat = materials["factory_glass"] if industrial else materials["glass"]
glass_batch = MeshBatch(name + "_Windows", collection, glass_mat)
edges = list(zip(footprint, footprint[1:] + footprint[:1]))
if industrial:
band_height = min(1.8, max(0.75, height * 0.16))
band_base = max(0.9, height * 0.52)
for start, end in edges:
add_wall_panel(glass_batch, start, end, band_base, band_height,
thickness=0.055, inset=0.12)
else:
floor_height = 3.25
floor_count = max(1, int((height - 0.7) / floor_height))
for floor in range(floor_count):
band_base = 0.55 + floor * floor_height + 0.95
if band_base + 1.25 > height - 0.18:
break
for start, end in edges:
add_wall_panel(glass_batch, start, end, band_base, 1.25,
thickness=0.045, inset=0.10)
glass_batch.finish()
def geometry_rings(geometry):
if not geometry:
return []
kind = geometry.get("type")
coordinates = geometry.get("coordinates", [])
if kind == "Polygon":
return coordinates[:1]
if kind == "MultiPolygon":
return [polygon[0] for polygon in coordinates if polygon]
return []
def feature_in_bounds(feature, projector):
def walk(value):
if isinstance(value, list) and value and isinstance(value[0], (int, float)):
return projector.inside(value)
return any(walk(v) for v in value) if isinstance(value, list) else False
return walk(feature.get("geometry", {}).get("coordinates", []))
def clip_polygon(ring, xmin, xmax, ymin, ymax):
"""Clip a projected polygon to the explicit OSM scene bounds."""
if len(ring) < 3:
return []
def clip_edge(points, inside, intersection):
if not points:
return []
result = []
previous = points[-1]
previous_inside = inside(previous)
for current in points:
current_inside = inside(current)
if current_inside != previous_inside:
result.append(intersection(previous, current))
if current_inside:
result.append(current)
previous = current
previous_inside = current_inside
return result
ring = clip_edge(
ring, lambda p: p[0] >= xmin,
lambda a, b: (xmin, a[1] + (b[1] - a[1]) * (xmin - a[0]) /
(b[0] - a[0]) if b[0] != a[0] else a[1]))
ring = clip_edge(
ring, lambda p: p[0] <= xmax,
lambda a, b: (xmax, a[1] + (b[1] - a[1]) * (xmax - a[0]) /
(b[0] - a[0]) if b[0] != a[0] else a[1]))
ring = clip_edge(
ring, lambda p: p[1] >= ymin,
lambda a, b: (a[0] + (b[0] - a[0]) * (ymin - a[1]) /
(b[1] - a[1]) if b[1] != a[1] else a[0], ymin))
ring = clip_edge(
ring, lambda p: p[1] <= ymax,
lambda a, b: (a[0] + (b[0] - a[0]) * (ymax - a[1]) /
(b[1] - a[1]) if b[1] != a[1] else a[0], ymax))
return ring
def add_geojson_layer(path, layer, projector, collection, material, z):
if not os.path.exists(path):
return 0
with open(path, "r", encoding="utf-8") as handle:
data = json.load(handle)
batch = MeshBatch("Road_" + layer, collection, material)
b = projector.bounds
xmin, ymin = projector.xy((b["min_lon"], b["min_lat"]))
xmax, ymax = projector.xy((b["max_lon"], b["max_lat"]))
count = 0
for feature in data.get("features", []):
if not feature_in_bounds(feature, projector):
continue
for ring in geometry_rings(feature.get("geometry")):
points = [projector.xy(pair) for pair in ring]
points = clip_polygon(points, xmin, xmax, ymin, ymax)
if len(points) >= 3:
batch.add_polygon(points, z)
count += 1
batch.finish()
return count
def add_polyline(name, coords, projector, collection, material, width, z):
points = [projector.xy(c) for c in coords]
if len(points) < 2:
return
curve = bpy.data.curves.new(name, "CURVE")
curve.dimensions = "3D"
curve.resolution_u = 1
curve.bevel_depth = width / 2
curve.bevel_resolution = 1
spline = curve.splines.new("POLY")
spline.points.add(len(points) - 1)
for point, (x, y) in zip(spline.points, points):
point.co = (x, y, z, 1)
obj = bpy.data.objects.new(name, curve)
collection.objects.link(obj)
obj.data.materials.append(material)
def parse_height(feature_tags, default):
try:
return max(0.5, float(feature_tags.get("height", default)))
except ValueError:
return default
def sample_tree_row(points, spacing, height):
if len(points) < 2:
return []
samples = [(points[0][0], points[0][1], height)]
distance_until_next = spacing
for start, end in zip(points, points[1:]):
dx = end[0] - start[0]
dy = end[1] - start[1]
segment_length = math.hypot(dx, dy)
if segment_length == 0:
continue
while distance_until_next <= segment_length:
ratio = distance_until_next / segment_length
samples.append((start[0] + dx * ratio, start[1] + dy * ratio, height))
distance_until_next += spacing
distance_until_next -= segment_length
last = points[-1]
if math.hypot(samples[-1][0] - last[0], samples[-1][1] - last[1]) > spacing * 0.45:
samples.append((last[0], last[1], height))
return samples
def add_tree_batch(positions, collection, trunk_material, leaf_material):
trunk = MeshBatch("Tree_Trunks", collection, trunk_material)
leaves = MeshBatch("Tree_Crowns", collection, leaf_material)
sides = 10
def add_blob(batch, cx, cy, cz, rx, ry, rz, phase):
rings = 5
start = len(batch.vertices)
for ring in range(rings):
latitude = -math.pi / 2 + math.pi * ring / (rings - 1)
ring_radius = math.cos(latitude)
for side in range(sides):
angle = math.tau * side / sides
variation = 1.0 + 0.09 * math.sin(phase + side * 1.73 + ring * 0.91)
batch.vertices.append((cx + rx * ring_radius * math.cos(angle) * variation,
cy + ry * ring_radius * math.sin(angle) * variation,
cz + rz * math.sin(latitude)))
for ring in range(rings - 1):
for side in range(sides):
next_side = (side + 1) % sides
batch.faces.append((start + ring * sides + side,
start + ring * sides + next_side,
start + (ring + 1) * sides + next_side,
start + (ring + 1) * sides + side))
for index, (x, y, height) in enumerate(positions):
base = len(trunk.vertices)
radius = max(0.12, height * 0.035)
trunk_top = height * 0.62
for z, ring_radius in ((0.0, radius), (trunk_top, radius * 0.68)):
for i in range(sides):
a = math.tau * i / sides
trunk.vertices.append((x + ring_radius * math.cos(a),
y + ring_radius * math.sin(a), z))
trunk.faces.append(tuple(base + i for i in range(sides - 1, -1, -1)))
for i in range(sides):
j = (i + 1) % sides
trunk.faces.append((base + i, base + j, base + sides + j, base + sides + i))
trunk.faces.append(tuple(base + sides + i for i in range(sides)))
crown_r = max(0.85, height * 0.30)
crown_z = height * 0.82
# Three overlapping blobs read as a natural crown at close range.
add_blob(leaves, x, y, crown_z, crown_r * 0.70,
crown_r * 0.62, crown_r * 0.72, index * 1.41)
add_blob(leaves, x - crown_r * 0.42, y + crown_r * 0.08,
crown_z * 0.98, crown_r * 0.52, crown_r * 0.48,
crown_r * 0.58, index * 2.17 + 0.7)
add_blob(leaves, x + crown_r * 0.40, y - crown_r * 0.05,
crown_z * 1.02, crown_r * 0.50, crown_r * 0.46,
crown_r * 0.55, index * 2.63 + 1.3)
trunk.finish()
leaves.finish()
def add_tree_model_instances(positions, collection):
"""Use the configured tree model for OSM tree nodes and rows.
The source asset is imported once as hidden template geometry. Every OSM
tree becomes a linked duplicate sharing the same mesh data instead of
copying the source geometry repeatedly. If the model is missing or import
fails, return False so the procedural tree builder can be used as a safe
fallback.
"""
if not positions or not os.path.exists(TREE_MODEL_PATH):
return False
before = set(bpy.data.objects)
try:
import_model(TREE_MODEL_PATH)
except Exception as exc:
print("TREE_MODEL_IMPORT_FAILED", exc)
return False
template_objects = [obj for obj in bpy.data.objects if obj not in before]
template_meshes = [obj for obj in template_objects if obj.type == "MESH"]
if not template_meshes:
for obj in template_objects:
bpy.data.objects.remove(obj, do_unlink=True)
return False
assign_tree_model_materials(template_meshes)
for obj in template_objects:
link_object_to_collection(obj, collection)
min_z = min(
(obj.matrix_world @ Vector(corner)).z
for obj in template_meshes
for corner in obj.bound_box
)
max_z = max(
(obj.matrix_world @ Vector(corner)).z
for obj in template_meshes
for corner in obj.bound_box
)
source_height = max(0.1, max_z - min_z)
# Normalize the hidden template so its base sits on z=0. The source model
# may not match the intended street-tree height, so every instance scales
# to the OSM height/default height below.
base_shift = Matrix.Translation((0.0, 0.0, -min_z))
for obj in template_objects:
obj.matrix_world = base_shift @ obj.matrix_world
obj.hide_viewport = True
obj.hide_render = True
obj.name = "Tree_Template_" + obj.name
for index, (x, y, height) in enumerate(positions):
target_height = max(4.6, min(9.2, height * 1.12))
scale = target_height / source_height
yaw = Matrix.Rotation((index * 1.61803398875) % math.tau, 4, "Z")
transform = Matrix.Translation((x, y, 0.0)) @ yaw @ Matrix.Diagonal(
(scale, scale, scale, 1.0)
)
for template in template_objects:
inst = template.copy()
if template.data:
inst.data = template.data
inst.animation_data_clear()
inst.matrix_world = transform @ template.matrix_world
inst.hide_viewport = False
inst.hide_render = False
inst.name = "Tree_Model_" + str(index)
collection.objects.link(inst)
return True
def polygon_area(ring):
if len(ring) < 3:
return 0.0
area = 0.0
for (x1, y1), (x2, y2) in zip(ring, ring[1:] + ring[:1]):
area += x1 * y2 - x2 * y1
return abs(area) * 0.5
def point_in_polygon(point, ring):
x, y = point
inside = False
j = len(ring) - 1
for i, (xi, yi) in enumerate(ring):
xj, yj = ring[j]
crosses = ((yi > y) != (yj > y))
if crosses:
x_at_y = (xj - xi) * (y - yi) / (yj - yi) + xi
if x < x_at_y:
inside = not inside
j = i
return inside
def add_scrub_patch(name, ring, material, collection):
"""Render OSM natural=scrub as textured, uneven shrub cover.
Important: do not use Poly Haven shrub model atlases as the surface
material here. Model atlases are laid out for a specific plant mesh, not
for tiling across an OSM polygon, and they appear as large patchwork blocks
in Cesium. Use a tileable foliage/grass texture for the ground-cover
surface, then add low deterministic domes for shrub volume.
"""
if len(ring) < 3:
return None
if ring[0] == ring[-1]:
ring = ring[:-1]
if len(ring) < 3:
return None
batch = MeshBatch(name, collection, material)
batch.add_polygon(ring, 0.055)
xmin = min(x for x, _ in ring)
xmax = max(x for x, _ in ring)
ymin = min(y for _, y in ring)
ymax = max(y for _, y in ring)
width = max(0.1, xmax - xmin)
depth = max(0.1, ymax - ymin)
area = polygon_area(ring)
clump_count = max(8, min(70, int(area / 24.0) + 6))
sides = 10
rings = 4
def add_dome(cx, cy, rx, ry, height, phase):
start = len(batch.vertices)
for ring_index in range(rings):
t = ring_index / (rings - 1)
z = 0.055 + height * math.sin(t * math.pi / 2)
radius_scale = math.cos(t * math.pi / 2)
for side in range(sides):
angle = math.tau * side / sides
wobble = 1.0 + 0.12 * math.sin(phase + side * 1.37 + ring_index * 0.73)
batch.vertices.append((
cx + rx * radius_scale * math.cos(angle) * wobble,
cy + ry * radius_scale * math.sin(angle) * wobble,
z,
))
for ring_index in range(rings - 1):
for side in range(sides):
next_side = (side + 1) % sides
batch.faces.append((
start + ring_index * sides + side,
start + ring_index * sides + next_side,
start + (ring_index + 1) * sides + next_side,
start + (ring_index + 1) * sides + side,
))
added = 0
attempts = 0
while added < clump_count and attempts < clump_count * 8:
attempts += 1
# Deterministic low-discrepancy sampling: stable between runs, but not
# grid-like. This avoids random scene churn while keeping natural spread.
u = (attempts * 0.61803398875) % 1.0
v = (attempts * 0.41421356237) % 1.0
x = xmin + u * width
y = ymin + v * depth
if not point_in_polygon((x, y), ring):
continue
scale = 0.65 + 0.55 * ((attempts * 0.754877666) % 1.0)
add_dome(x, y, 0.82 * scale, 0.64 * scale,
0.18 + 0.14 * scale, attempts * 0.91)
added += 1
obj = batch.finish()
if obj:
obj["scrub_texture"] = "Poly Haven leafy_grass, scrub-tinted"
obj["scrub_style"] = "tileable foliage ground cover with low shrub domes"
return obj
def link_object_to_collection(obj, collection):
for current in list(obj.users_collection):
current.objects.unlink(obj)
collection.objects.link(obj)
def add_fountain(name, x, y, collection, materials):
"""Create a restrained fountain from an explicit amenity=fountain node."""
def cylinder(part_name, radius, depth, z, material, vertices=48):
bpy.ops.mesh.primitive_cylinder_add(
vertices=vertices, radius=radius, depth=depth,
location=(x, y, z))
obj = bpy.context.object
obj.name = name + "_" + part_name
link_object_to_collection(obj, collection)
obj.data.materials.append(material)
for polygon in obj.data.polygons:
polygon.use_smooth = True
return obj
basin = cylinder("Basin", 3.0, 0.32, 0.16,
materials["fountain_stone"])
basin["osm_feature"] = "amenity=fountain"
cylinder("Water", 2.52, 0.045, 0.335,
materials["fountain_water"])
cylinder("Pedestal", 0.30, 0.78, 0.72,
materials["fountain_stone"], vertices=32)
# A compact central spray reads at overview distance without creating a
# high-polygon particle system that would be expensive in Cesium.
bpy.ops.mesh.primitive_uv_sphere_add(
segments=20, ring_count=10, radius=0.22,
location=(x, y, 1.30))
crown = bpy.context.object
crown.name = name + "_Water_Crown"
link_object_to_collection(crown, collection)
crown.data.materials.append(materials["fountain_spray"])
for index in range(8):
angle = math.tau * index / 8.0
radius = 0.50
bpy.ops.mesh.primitive_uv_sphere_add(
segments=12, ring_count=6, radius=0.075,
location=(x + math.cos(angle) * radius,
y + math.sin(angle) * radius,
1.02 + 0.10 * math.sin(angle * 2.0)))
droplet = bpy.context.object
droplet.name = name + "_Droplet_" + str(index + 1)
link_object_to_collection(droplet, collection)
droplet.data.materials.append(materials["fountain_spray"])
def look_at(obj, target):
obj.rotation_euler = (Vector(target) - obj.location).to_track_quat("-Z", "Y").to_euler()
def clear_scene():
bpy.ops.object.select_all(action="SELECT")
bpy.ops.object.delete(use_global=False)
for collection in list(bpy.data.collections):
if collection.name != "Collection" and collection.users == 0:
bpy.data.collections.remove(collection)
def configure_scene():
scene = bpy.context.scene
scene.render.engine = "BLENDER_EEVEE_NEXT"
scene.render.resolution_x = 1200
scene.render.resolution_y = 900
scene.render.resolution_percentage = 100
scene.render.image_settings.file_format = "PNG"
scene.render.film_transparent = False
scene.world.color = (0.055, 0.075, 0.095)
scene.view_settings.look = "AgX - Medium High Contrast"
def configure_default_viewport():
"""Make the saved Layout workspace show the same rendered camera view."""
workspace = bpy.data.workspaces.get("Layout")
if workspace:
try:
bpy.context.window.workspace = workspace
except (AttributeError, RuntimeError):
pass
for screen in bpy.data.screens:
for area in screen.areas:
if area.type != "VIEW_3D":
continue
space = area.spaces.active
space.shading.type = "MATERIAL"
space.shading.light = "STUDIO"
space.shading.color_type = "MATERIAL"
space.overlay.show_floor = False
if space.region_3d:
space.region_3d.view_perspective = "CAMERA"
space.region_3d.view_camera_zoom = 0.0
def build(args):
bounds, ways, point_features = parse_osm(args["osm"])
projector = Projector(bounds)
clear_scene()
configure_scene()
ground_c = new_collection("00_Ground")
water_c = new_collection("01_Water")
green_c = new_collection("02_Green")
roads_c = new_collection("03_Roads")
buildings_c = new_collection("04_Buildings")
props_c = new_collection("05_Props")
ground_mat = make_material("Ground", (0.27, 0.32, 0.24))
water_mat = make_material("Lake Water", (0.035, 0.22, 0.30), 0.18, 0.05)
grass_mat = make_textured_material(
"Grass", "leafy_grass_diff_1k.jpg", "leafy_grass_nor_gl_1k.jpg",
roughness=0.92, scale=7.0, tint=(0.12, 0.48, 0.08), tint_factor=0.72)
scrub_mat = make_textured_material(
"Scrub Ground Cover", "leafy_grass_diff_1k.jpg",
"leafy_grass_nor_gl_1k.jpg", roughness=0.96, scale=13.0,
tint=(0.05, 0.34, 0.08), tint_factor=0.42)
fountain_mats = {
"fountain_stone": make_material("Fountain Stone", (0.42, 0.45, 0.43), 0.72),
"fountain_water": make_material("Fountain Water", (0.03, 0.32, 0.42), 0.16, 0.05),
"fountain_spray": make_material("Fountain Spray", (0.20, 0.70, 0.78), 0.12, 0.02),
}
building_mats = {
"default": make_textured_material(
"Office White Plaster Facade", "white_plaster_02_diff_1k.jpg",
"white_plaster_02_nor_gl_1k.jpg", roughness=0.82,
scale=4.2, metallic=0.0, tint=(0.92, 0.94, 0.92),
tint_factor=0.38),
"industrial": make_textured_material(
"Industrial White Ribbed Facade", "corrugated_iron_03_diff_1k.jpg",
"corrugated_iron_03_nor_gl_1k.jpg", roughness=0.56,
scale=2.4, metallic=0.16, tint=(0.86, 0.92, 0.94),
tint_factor=0.68),
"office_roof": make_textured_material(
"Office Light Flat Roof", "concrete_floor_02_diff_1k.jpg",
"concrete_floor_02_bump_1k.jpg", roughness=0.84,
scale=5.0, normal_is_bump=True, tint=(0.82, 0.86, 0.88),
tint_factor=0.35),
"industrial_roof": make_textured_material(
"Factory Blue Metal Roof", "blue_metal_plate_diff_1k.jpg",
"blue_metal_plate_nor_gl_1k.jpg", roughness=0.48,
scale=3.4, metallic=0.28, tint=(0.03, 0.42, 0.78),
tint_factor=0.45),
"glass": make_material("Office Blue Gray Glass", (0.12, 0.20, 0.24), 0.22, 0.10),
"factory_glass": make_material("Factory Dark Windows", (0.10, 0.14, 0.15), 0.28, 0.08),
}
road_mats = {
"road_surface": make_material("Road Asphalt", (0.055, 0.065, 0.070)),
"intersection_surface": make_material("Intersection Asphalt", (0.065, 0.075, 0.080)),
"sidewalks": make_material("Sidewalk", (0.49, 0.51, 0.49)),
"sidewalk_corners": make_material("Sidewalk Corner", (0.49, 0.51, 0.49)),
"lane_separators": make_material("Lane Separator", (0.85, 0.84, 0.72)),
"center_lines": make_material("Center Line", (0.94, 0.58, 0.06)),
"crosswalks": make_material("Crosswalk", (0.95, 0.94, 0.82)),
"vehicle_stop_lines": make_material("Stop Line", (0.95, 0.94, 0.82)),
"lane_arrows_webscale": make_material("Lane Arrow", (0.95, 0.94, 0.82)),
}
b = bounds
scene_xmin, scene_ymin = projector.xy((b["min_lon"], b["min_lat"]))
scene_xmax, scene_ymax = projector.xy((b["max_lon"], b["max_lat"]))
ground_ring = [projector.xy((b["min_lon"] - 0.0012, b["min_lat"] - 0.0012)),
projector.xy((b["max_lon"] + 0.0012, b["min_lat"] - 0.0012)),
projector.xy((b["max_lon"] + 0.0012, b["max_lat"] + 0.0012)),
projector.xy((b["min_lon"] - 0.0012, b["max_lat"] + 0.0012))]
ground_batch = MeshBatch("Ground Plane", ground_c, ground_mat)
ground_batch.add_polygon(ground_ring, -0.35)
ground_batch.finish()
grass_rings = []
tree_rows = []
lake_count = 0
grass_count = 0
scrub_count = 0
fountain_count = 0
building_count = 0
industrial_count = 0
focus_points = []
for way in ways:
coords = way["coords"]
if not any(projector.inside(c) for c in coords):
continue
ring = projector.ring(coords)
tag = way["tags"]
if tag.get("natural") == "water" or tag.get("water") == "lake":
ring = clip_polygon(ring, scene_xmin, scene_xmax,
scene_ymin, scene_ymax)
batch = MeshBatch("Lake Surface", water_c, water_mat)
if len(ring) >= 3:
batch.add_polygon(ring, 0.10)
batch.finish()
lake_count += 1
elif tag.get("landuse") == "grass":
ring = clip_polygon(ring, scene_xmin, scene_xmax,
scene_ymin, scene_ymax)
grass_rings.append(ring)
focus_points.extend(ring)
batch = MeshBatch("Grass_" + str(way["id"]), green_c, grass_mat)
if len(ring) >= 3:
batch.add_polygon(ring, 0.015)
batch.finish()
grass_count += 1
elif tag.get("natural") == "scrub" and len(ring) >= 3:
ring = clip_polygon(ring, scene_xmin, scene_xmax,
scene_ymin, scene_ymax)
focus_points.extend(ring)
if len(ring) >= 3:
add_scrub_patch("Scrub_" + str(way["id"]), ring,
scrub_mat, green_c)
scrub_count += 1
elif tag.get("natural") == "tree_row":
tree_rows.append((ring, tag))
focus_points.extend(ring)
elif "building" in tag and len(ring) >= 3:
way_id = str(way["id"])
industrial = (tag.get("building") == "industrial" and
way_id not in OFFICE_OVERRIDE_WAY_IDS)
source_height = max(3.0, parse_height(tag, 12.0))
# Ordinary park offices are represented as three floors plus roof;
# retain explicit high-rise massing and all industrial heights.
height = source_height if industrial or source_height >= 30.0 else 11.4
material = building_mats["industrial"] if industrial else building_mats["default"]
building_name = "Building_" + way_id
building_obj = make_prism(building_name, ring, 0.08, height,
material, buildings_c)
if building_obj:
building_obj["osm_height"] = source_height
building_obj["render_height"] = height
building_obj["building_kind"] = "industrial" if industrial else "office"
building_obj["osm_building_tag"] = tag.get("building", "")
building_obj["office_override"] = way_id in OFFICE_OVERRIDE_WAY_IDS
bevel = building_obj.modifiers.new("Soft facade edges", "BEVEL")
bevel.width = 0.16
bevel.segments = 2
roof_mat = (building_mats["industrial_roof"] if industrial
else building_mats["office_roof"])
add_roof(building_name, ring, height + 0.095, roof_mat, buildings_c)
add_building_details(building_name, ring, height, industrial,
building_mats, buildings_c)
focus_points.extend(ring)
building_count += 1
industrial_count += int(industrial)
# Fine road surfaces and markings from the osm2streets GeoJSON output.
layer_z = {"road_surface": 0.03, "intersection_surface": 0.035,
"sidewalks": 0.065, "sidewalk_corners": 0.067,
"lane_separators": 0.090, "center_lines": 0.092,
"crosswalks": 0.094, "vehicle_stop_lines": 0.096,
"lane_arrows_webscale": 0.098}
road_counts = {}
for layer, z in layer_z.items():
road_counts[layer] = add_geojson_layer(
os.path.join(args["geojson"], layer + ".geojson"), layer,
projector, roads_c, road_mats[layer], z)
# If no road GeoJSON is available, retain a useful OSM-only fallback.
if road_counts.get("road_surface", 0) == 0:
for way in ways:
highway = way["tags"].get("highway")
if highway and len(way["coords"]) >= 2:
width = {"secondary": 7.0, "residential": 5.5, "service": 3.5}.get(highway, 4.0)
add_polyline("OSM_Road_" + str(way["id"]), way["coords"], projector,
roads_c, road_mats["road_surface"], width, 0.03)
# Trees come only from explicit OSM tree nodes and tree_row ways.
trees = []
individual_tree_count = 0
for feature in point_features:
if feature["tags"].get("natural") != "tree":
continue
if not projector.inside(feature["coord"]):
continue
x, y = projector.xy(feature["coord"])
trees.append((x, y, parse_height(feature["tags"], 5.5)))
individual_tree_count += 1
row_tree_count = 0
for row, row_tags in tree_rows:
row_samples = sample_tree_row(row, spacing=5.0,
height=parse_height(row_tags, 5.0))
trees.extend(row_samples)
row_tree_count += len(row_samples)
if trees:
if not add_tree_model_instances(trees, props_c):
tree_trunk = make_textured_material(
"Tree Trunk", "bark_brown_01_diff_1k.jpg",
"bark_brown_01_nor_gl_1k.jpg", roughness=0.92, scale=5.0)
tree_leaf = make_material("Tree Crown", (0.08, 0.30, 0.09), 0.88)
add_procedural_surface(tree_leaf,
((0.035, 0.16, 0.045), (0.12, 0.42, 0.13)),
scale=2.8, detail=3.2, bump_strength=0.10)
add_tree_batch(trees, props_c, tree_trunk, tree_leaf)
for feature in point_features:
if feature["tags"].get("amenity") != "fountain":
continue
if not projector.inside(feature["coord"]):
continue
fx, fy = projector.xy(feature["coord"])
add_fountain("Fountain_" + str(feature["id"]), fx, fy,
props_c, fountain_mats)
fountain_count += 1
# A simple sun/area-light rig keeps the model readable in viewport and render.
bpy.ops.object.light_add(type="SUN", location=(0, 0, 500))
sun = bpy.context.object
sun.name = "Sun"
sun.data.energy = 3.0
sun.rotation_euler = (math.radians(28), math.radians(-22), math.radians(-32))
bpy.ops.object.light_add(type="AREA", location=(0, -220, 420))
area = bpy.context.object
area.name = "Fill Light"
area.data.energy = 1700
area.data.shape = "DISK"
area.data.size = 260
look_at(area, (0, 0, 0))
width = (b["max_lon"] - b["min_lon"]) * projector.m_per_lon
height = (b["max_lat"] - b["min_lat"]) * projector.m_per_lat
if focus_points:
min_fx = min(point[0] for point in focus_points)
max_fx = max(point[0] for point in focus_points)
min_fy = min(point[1] for point in focus_points)
max_fy = max(point[1] for point in focus_points)
focus_x = (min_fx + max_fx) / 2
focus_y = (min_fy + max_fy) / 2
focus_span = max(max_fx - min_fx, (max_fy - min_fy) * 1.25)
cam_location = (focus_x + focus_span * 0.78,
focus_y - focus_span * 0.92,
focus_span * 1.22)
camera_target = (focus_x, focus_y, 3)
else:
cam_location = (width * 0.78, -height * 1.15, max(width, height) * 1.22)
camera_target = (0, 0, 3)
bpy.ops.object.camera_add(location=cam_location)
camera = bpy.context.object
camera.name = "Park Overview Camera"
camera.data.lens = 48
camera.data.clip_start = 0.1
camera.data.clip_end = 5000.0
look_at(camera, camera_target)
bpy.context.scene.camera = camera
configure_default_viewport()
scene = bpy.context.scene
scene.render.filepath = args["render"]
scene["source_osm"] = args["osm"]
scene["source_geojson"] = args["geojson"]
scene["osm_bounds"] = json.dumps(bounds, ensure_ascii=True)
scene["building_count"] = building_count
scene["industrial_building_count"] = industrial_count
scene["office_override_way_ids"] = json.dumps(sorted(OFFICE_OVERRIDE_WAY_IDS))
scene["lake_count"] = lake_count
scene["grass_count"] = grass_count
scene["scrub_count"] = scrub_count
scene["fountain_count"] = fountain_count
scene["tree_node_count"] = individual_tree_count
scene["tree_row_count"] = row_tree_count
scene["tree_count"] = len(trees)
scene["road_feature_counts"] = json.dumps(road_counts, ensure_ascii=True)
os.makedirs(os.path.dirname(args["output"]), exist_ok=True)
os.makedirs(os.path.dirname(args["render"]), exist_ok=True)
bpy.ops.file.pack_all()
bpy.ops.wm.save_as_mainfile(filepath=args["output"])
bpy.ops.render.render(write_still=True)
print("NANTAIZI_DONE", json.dumps({"output": args["output"],
"render": args["render"],
"buildings": building_count,
"industrial_buildings": industrial_count,
"lake": lake_count,
"grass": grass_count,
"scrub": scrub_count,
"fountains": fountain_count,
"tree_nodes": individual_tree_count,
"tree_row_instances": row_tree_count,
"trees": len(trees),
"road_features": road_counts}, ensure_ascii=True))
if __name__ == "__main__":
build(cli_args())