【小沐学GIS】基于Openstreetmap创建Sionna RT场景（Python）

文章目录

• 1、简介
• • 1.1 blender
• 2、下载和安装
• • 2.1 Python
  • 2.2 jupyter
• 3、运行
• 结语

1、简介

1.1 blender

https://www.blender.org/

Blender 是一款免费开源的3D创作套件。

使用 Blender，您可以创建3D可视化效果，例如静态图像、3D动画、VFX（视觉特效）快照和视频编辑。它非常适合那些受益于其统一解决方案和响应式开发过程的独立和小型工作室。

Blender 是一款跨平台的应用工具，可以在 Linux、macOS 以及 Windows 系统下运行。与其他三维建模工具相比，Blender 对内存和驱动的需求更低。其界面使用 OpenGL，在所有支持的硬件与平台都能提供一致的用户体验。

Blender 是一个完整集成的 3D 创作套件，提供了大量的基础工具，包括建模、渲染、动画 & 绑定、视频编辑、视觉效果、合成、贴图，以及多种类型的模拟。

跨平台，使用了 OpenGL 的 GUI 可以在所有主流平台上都表现出一致的显示效果（并且可通过 Python 脚本来自定义界面）。

2、下载和安装

2.1 Python

https://www.python.org/
Python是一种广泛使用的高级编程语言，它以其清晰的语法和代码可读性而闻名。Python支持多种编程范式，包括面向对象、命令式、函数式和过程式编程。

安装后打印Python版本信息如下：

python -v

pip -v

2.2 jupyter

Jupyter Notebook 是一个基于 Web 的交互式计算环境，支持多种编程语言，包括 Python、R、Julia 等。它的主要功能是将代码、文本、数学方程式、可视化和其他相关元素组合在一起，创建一个动态文档，用于数据分析、机器学习、科学计算和数据可视化等方面。Jupyter Notebook 提供了一个交互式的界面，使用户能够以增量和可视化的方式构建和执行代码，同时支持 Markdown 格式的文本和 LaTeX 数学符号。

安装 Jupyter Notebook：在命令提示符中输入以下命令，使用 pip 安装 Jupyter Notebook。

pip install jupyter notebook

启动 Jupyter Notebook：在命令提示符中输入以下命令，启动 Jupyter Notebook。

jupyter notebook

接下来Jupyter Notebook 会在默认的浏览器中打开，如果没有自动打开，可以在浏览器中输入 http://localhost:8888/tree 来访问。

http://localhost:8888/tree

新建文件如下：

界面如下：

3、运行

https://github.com/manoj-kumar-joshi/sionna_osm_scene

从 Openstreetmap 数据创建 Sionna 光线追踪的场景文件

这是一个实用程序项目，它使用 Openstreetmap 数据生成与 Mitsuba 兼容的 3D 场景文件以及建筑物、地面和道路的 3D 对象。使用 OSM_to_Sionna 生成的文件可以直接用作 Sionna Ray 追踪应用程序的输入。

打开文件：
https://github.com/manoj-kumar-joshi/sionna_osm_scene/blob/main/OSM_to_Sionna.ipynb

测试代码如下：

import ipyleaflet
import IPython.display
import ipyvolume.pylab as p3
import pyproj
import shapely
from shapely.geometry import shape
from shapely.ops import transform
import math
import pyvista as pv
import numpy as np
import osmnx as ox
from shapely.geometry import Polygon, Point, LineString
import os
from pyproj import Transformer
import open3d as o3d
import xml.etree.ElementTree as ET
import xml.dom.minidom as minidom

初始化Sionna Scene XML对象并添加默认值：

# Set the center position lat and lon as a starting point. Use any string of your choice for LOCATION_STR
center_lat = 14.557097311312177
center_lon = 121.02000072883868
LOCATION_STR = "PHILLIPINES

# Set up default values for resolution
spp_default = 4096
resx_default = 1024
resy_default = 768

# Define camera settings
camera_settings = {
    "rotation": (0, 0, -90),  # Assuming Z-up orientation
    "fov": 42.854885
}

# Define material colors. This is RGB 0-1 formar https://rgbcolorpicker.com/0-1
material_colors = {
    "mat-itu_concrete": (0.539479, 0.539479, 0.539480),
    "mat-itu_marble": (0.701101, 0.644479, 0.485150),
    "mat-itu_metal": (0.219526, 0.219526, 0.254152),
    "mat-itu_wood": (0.043, 0.58, 0.184),
    "mat-itu_wet_ground": (0.91,0.569,0.055),
}
transformer = Transformer.from_crs("EPSG:4326", "EPSG:26915")
center_26915 = transformer.transform(center_lat,center_lon)
sionna_center_x = center_26915[0]
sionna_center_y = center_26915[1]
sionna_center_z = 0

scene = ET.Element("scene", version="2.1.0")
# Add defaults
ET.SubElement(scene, "default", name="spp", value=str(spp_default))
ET.SubElement(scene, "default", name="resx", value=str(resx_default))
ET.SubElement(scene, "default", name="resy", value=str(resy_default))
# Add integrator
integrator = ET.SubElement(scene, "integrator", type="path")
ET.SubElement(integrator, "integer", name="max_depth", value="12")

# Define materials
for material_id, rgb in material_colors.items():
    bsdf_twosided = ET.SubElement(scene, "bsdf", type="twosided", id=material_id)
    bsdf_diffuse = ET.SubElement(bsdf_twosided, "bsdf", type="diffuse")
    ET.SubElement(bsdf_diffuse, "rgb", value=f"{rgb[0]} {rgb[1]} {rgb[2]}", name="reflectance")

# Add emitter
emitter = ET.SubElement(scene, "emitter", type="constant", id="World")
ET.SubElement(emitter, "rgb", value="1.000000 1.000000 1.000000", name="radiance")

# Add camera (sensor)
sensor = ET.SubElement(scene, "sensor", type="perspective", id="Camera")
ET.SubElement(sensor, "string", name="fov_axis", value="x")
ET.SubElement(sensor, "float", name="fov", value=str(camera_settings["fov"]))
ET.SubElement(sensor, "float", name="principal_point_offset_x", value="0.000000")
ET.SubElement(sensor, "float", name="principal_point_offset_y", value="-0.000000")
ET.SubElement(sensor, "float", name="near_clip", value="0.100000")
ET.SubElement(sensor, "float", name="far_clip", value="10000.000000")
sionna_transform = ET.SubElement(sensor, "transform", name="to_world")
ET.SubElement(sionna_transform, "rotate", x="1", angle=str(camera_settings["rotation"][0]))
ET.SubElement(sionna_transform, "rotate", y="1", angle=str(camera_settings["rotation"][1]))
ET.SubElement(sionna_transform, "rotate", z="1", angle=str(camera_settings["rotation"][2]))
camera_position = np.array([0, 0, 100])  # Adjust camera height
ET.SubElement(sionna_transform, "translate", value=" ".join(map(str, camera_position)))
sampler = ET.SubElement(sensor, "sampler", type="independent")
ET.SubElement(sampler, "integer", name="sample_count", value="$spp")
film = ET.SubElement(sensor, "film", type="hdrfilm")
ET.SubElement(film, "integer", name="width", value="$resx")
ET.SubElement(film, "integer", name="height", value="$resy")

打开交互式地图，选择要使用的区域：

m = ipyleaflet.Map(center=(center_lat, center_lon), zoom=16)
dc = ipyleaflet.DrawControl()
m.add(dc)
m

# Get coordinates in meter for the area of interst polygon (This will be used in next steps)
wsg84 = pyproj.CRS("epsg:4326")
lambert = pyproj.CRS("epsg:26915")
transformer = pyproj.Transformer.from_crs(wsg84, lambert, always_xy=True)
coords = [transformer.transform(x, y) for x, y in dc.last_draw['geometry']['coordinates'][0]]
print(coords)
print(dc.last_draw['geometry']['coordinates'][0])
aoi_polygon = shapely.geometry.Polygon(coords)

# Store center of the selected area to be used in calculations later on
center_x = aoi_polygon.centroid.x
center_y = aoi_polygon.centroid.y

# Set Location of the directory where scene and objects will be stored
LOCATION_DIR = f"{LOCATION_STR}_{center_x}_{center_y}"
print(LOCATION_DIR)

# Create Directories
os.mkdir(f"d:/simple_scene")
os.mkdir(f"d:/simple_scene/mesh")  

# Utility Function
def points_2d_to_poly(points, z):
    """Convert a sequence of 2d coordinates to a polydata with a polygon."""
    faces = [len(points), *range(len(points))]
    poly = pv.PolyData([p + (z,) for p in points], faces=faces)
    return poly

创建地面网格并添加到场景中：

# Utility Function
def points_2d_to_poly(points, z):
    """Convert a sequence of 2d coordinates to a polydata with a polygon."""
    faces = [len(points), *range(len(points))]
    poly = pv.PolyData([p + (z,) for p in points], faces=faces)
    return poly

wsg84 = pyproj.CRS("epsg:4326")
lambert = pyproj.CRS("epsg:26915")
transformer = pyproj.Transformer.from_crs(wsg84, lambert, always_xy=True)
coords = [transformer.transform(x, y) for x, y in dc.last_draw['geometry']['coordinates'][0]]

ground_polygon = shapely.geometry.Polygon(coords)
z_coordinates = np.full(len(ground_polygon.exterior.coords), 0)  # Assuming the initial Z coordinate is zmin
exterior_coords = ground_polygon.exterior.coords
oriented_coords = list(exterior_coords)
# Ensure counterclockwise orientation
if ground_polygon.exterior.is_ccw:
    oriented_coords.reverse()
points = [(coord[0]-center_x, coord[1]-center_y) for coord in oriented_coords]
# bounding polygon
boundary_points_polydata = points_2d_to_poly(points, z_coordinates[0])
edge_polygon = boundary_points_polydata
footprint_plane = edge_polygon.delaunay_2d()
footprint_plane.points[:] = (footprint_plane.points - footprint_plane.center)*1.5 + footprint_plane.center
pv.save_meshio(f"d:/simple_scene/mesh/ground.ply",footprint_plane)

material_type = "mat-itu_wet_ground"
sionna_shape = ET.SubElement(scene, "shape", type="ply", id=f"mesh-ground")
ET.SubElement(sionna_shape, "string", name="filename", value=f"mesh/ground.ply")
bsdf_ref = ET.SubElement(sionna_shape, "ref", id=material_type, name="bsdf")
ET.SubElement(sionna_shape, "boolean", name="face_normals",value="true")

创建建筑网格并添加到场景中：

import osmnx as ox
wsg84 = pyproj.CRS("epsg:4326")
lambert = pyproj.CRS("epsg:4326")
transformer = pyproj.Transformer.from_crs(wsg84, lambert, always_xy=True)
coords = [transformer.transform(x, y) for x, y in dc.last_draw['geometry']['coordinates'][0]]

osm_polygon = shapely.geometry.Polygon(coords)
# Query the OpenStreetMap data
buildings = ox.geometries.geometries_from_polygon(osm_polygon, tags={'building': True})

# Filter buildings that intersect with the polygon
filtered_buildings = buildings[buildings.intersects(osm_polygon)]

filtered_buildings.head(5)

以下代码使用建筑足迹并拉伸它们来创建三角形网格，并逐一添加Sionna场景。

buildings_list = filtered_buildings.to_dict('records')
source_crs = pyproj.CRS(filtered_buildings.crs)
target_crs = pyproj.CRS('EPSG:26915')
transformer = pyproj.Transformer.from_crs(source_crs, target_crs, always_xy=True).transform
for idx, building in enumerate(buildings_list):
    # Convert building geometry to a shapely polygon
    building_polygon = shape(building['geometry'])
    if building_polygon.geom_type != 'Polygon':
        continue
    building_polygon = transform(transformer, building_polygon)
    if math.isnan(float(building['building:levels'])):
        building_height = 3.5
    else:
        building_height = int(building['building:levels']) * 3.5
    z_coordinates = np.full(len(building_polygon.exterior.coords), 0)  # Assuming the initial Z coordinate is zmin
    exterior_coords = building_polygon.exterior.coords
    oriented_coords = list(exterior_coords)
    # Ensure counterclockwise orientation
    if building_polygon.exterior.is_ccw:
        oriented_coords.reverse()
    points = [(coord[0]-center_x, coord[1]-center_y) for coord in oriented_coords]
    # bounding polygon
    boundary_points_polydata = points_2d_to_poly(points, z_coordinates[0])
    edge_polygon = boundary_points_polydata
    footprint_plane = edge_polygon.delaunay_2d()
    footprint_plane = footprint_plane.triangulate()
    footprint_3D = footprint_plane.extrude((0, 0, building_height), capping=True)
    footprint_3D.save(f"d:/simple_scene/mesh/building_{idx}.ply")
    local_mesh = o3d.io.read_triangle_mesh(f"d:/simple_scene/mesh/building_{idx}.ply")
    o3d.io.write_triangle_mesh(f"d:/simple_scene/mesh/building_{idx}.ply", local_mesh)
    material_type = "mat-itu_marble"
    # Add shape elements for PLY files in the folder
    sionna_shape = ET.SubElement(scene, "shape", type="ply", id=f"mesh-building_{idx}")
    ET.SubElement(sionna_shape, "string", name="filename", value=f"mesh/building_{idx}.ply")
    bsdf_ref = ET.SubElement(sionna_shape, "ref", id= material_type, name="bsdf")
    ET.SubElement(sionna_shape, "boolean", name="face_normals",value="true")

创建道路网格并添加到场景中:

def convert_lane_to_numeric(lane):
    try:
        return int(lane)
    except ValueError:
        try:
            return float(lane)
        except ValueError:
            return None

        # Helper function to calculate edge geometry if missing
def calculate_edge_geometry(u, v, data):
    u_data = graph.nodes[u]
    v_data = graph.nodes[v]
    return LineString([(u_data['x'], u_data['y']), (v_data['x'], v_data['y'])])

G = ox.graph_from_polygon(polygon = osm_polygon, simplify= False, retain_all=True,truncate_by_edge=True,network_type = 'all_private')
graph = ox.project_graph(G, to_crs='epsg:26915')
ox.plot_graph(graph)

现在，使用车道作为参数将每条线段转换为道路网格，以设置道路宽度:

# Create a list to store GeoDataFrames for each road segment
gdf_roads_list = []
# Set the fixed Z coordinate for the buffer polygons
Z0 = .25  # You can adjust this value based on the desired elevation of the roads
# Create a list to store the meshes
mesh_list = []
mesh_collection = pv.PolyData()
# Iterate over each edge in the graph
for u, v, key, data in graph.edges(keys=True, data=True):
    # Check if the edge has geometry, otherwise create geometries from the nodes
    if 'geometry' not in data:
        data['geometry'] = calculate_edge_geometry(u, v, data)

    # Get the lanes attribute for the edge
    lanes = data.get('lanes', 1)  # Default to 1 lane if lanes attribute is not available

    if not isinstance(lanes, list):
        lanes = [lanes]
        
    # Convert lane values to numeric (integers or floats) using the helper function
    num_lanes = [convert_lane_to_numeric(lane) for lane in lanes]

    # Filter out None values (representing non-numeric lanes) and calculate the road width
    num_lanes = [lane for lane in num_lanes if lane is not None]
    road_width = num_lanes[0] * 3.5
    # Buffer the LineString with the road width and add Z coordinate
    line_buffer = data['geometry'].buffer(road_width)
    # Convert the buffer polygon to a PyVista mesh
    exterior_coords = line_buffer.exterior.coords
    z_coordinates = np.full(len(line_buffer.exterior.coords), Z0)
    oriented_coords = list(exterior_coords)
    # Ensure counterclockwise orientation
    if line_buffer.exterior.is_ccw:
        oriented_coords.reverse()
    points = [(coord[0]-center_x, coord[1]-center_y) for coord in oriented_coords]
    # bounding polygon
    boundary_points_polydata = points_2d_to_poly(points, z_coordinates[0])
    mesh = boundary_points_polydata.delaunay_2d()
    # Add the mesh to the list
    mesh_collection = mesh_collection + mesh
    mesh_list.append(mesh)
output_file = f"d:/simple_scene/mesh/road_mesh_combined.ply"
pv.save_meshio(output_file,mesh_collection)
material_type = "mat-itu_concrete"
# Add shape elements for PLY files in the folder
sionna_shape = ET.SubElement(scene, "shape", type="ply", id=f"mesh-roads_{idx}")
ET.SubElement(sionna_shape, "string", name="filename", value=f"mesh/road_mesh_combined.ply")
bsdf_ref = ET.SubElement(sionna_shape, "ref", id= material_type, name="bsdf")
ET.SubElement(sionna_shape, "boolean", name="face_normals",value="true")

最后保存场景文件:

# Create and write the XML file
tree = ET.ElementTree(scene)
xml_string = ET.tostring(scene, encoding="utf-8")
xml_pretty = minidom.parseString(xml_string).toprettyxml(indent="    ")  # Adjust the indent as needed

with open(f"d:/simple_scene/simple_OSM_scene.xml", "w", encoding="utf-8") as xml_file:
    xml_file.write(xml_pretty)

生成模型文件如下：


在blender加载上面结果文件如下：

结语

如果您觉得该方法或代码有一点点用处，可以给作者点个赞，或打赏杯咖啡；╮(￣▽￣)╭
如果您感觉方法或代码不咋地//(ㄒoㄒ)//，就在评论处留言，作者继续改进；o_O???
如果您需要相关功能的代码定制化开发，可以留言私信作者；(✿◡‿◡)
感谢各位大佬童鞋们的支持！( ´ ▽´ )ﾉ ( ´ ▽´)っ！！！