根据地图数据程序化城市建模 Houdini

花了大概一个礼拜左右的时间把城市建模的插件做出来了，还有一些瑕疵，但是已经达到预期想要的效果了。说在最前面的是，做这个东西完全没有任何商业目的，只是希望在houdini里面掌握更多程序化的思维。

先看效果：

 

整体思路分为：

1；抓取地图数据

2：分析地图数据，并根据坐标位置来生成最基本的建筑和道路的几何形状

3：根据道路和现有建筑布局，用点填充其他空余空间

4：程序化生成建筑楼，并将上一步的点提换成建筑

 

下面将每一步的具体方法和思路讲解一下：

1：抓取地图数据：

现在市面上能够免费得到地图所有元素坐标和类别信息的机构我只找到了openstreetmap.org这一家，他们这个项目的目的也是为了能够向世界所有人创造并提供免费的地理数据。我们可以直接在他的地图网页上框选一定的地图空间然后导出相应经纬度的地图坐标数据。非常重要的一点事他们目前只提供osm格式文件，其实就是xml的一种，如果想用编辑器直接查看里面的内容，可以吧后后缀名.osm改为.xml就可以打开了。

地图相关资源：

地图链接：https://www.openstreetmap.org

关于地图到处的xml数据结构：http://wiki.openstreetmap.org/wiki/OSM_XML

关于xml几个重要标记的解释：http://wiki.openstreetmap.org/wiki/Elements

地图中的元素一览表：http://wiki.openstreetmap.org/wiki/Map_Features

 

这里简单讲解一下.osm格式的结构:

<osm>
	<bounds minlat="***" minlon="***" maxlat="***" maxlon="***"/>
	<node id="***" ******** lat="40.4444343" lon="-79.9486641"/>
	<node/>
	....
	<node/>
	<way id="27574365" ....>
  		<nd ref="1705867277"/>
  		<nd ref="302440336"/>
  		<tag k="highway" v="service"/>
 	</way>
 	....
 	<relation/>
 	...
 </osm>

 

bounds定义了整个数据的上下左右边界经纬度，Node相当于是地图中记录下来了的每一个点的数据，他可以是道路上的也可以是建筑上的。way是把所有点集结起来的一种方法，通过对他的tag做标签能够确定单个way下面所有的点是某一条街道还是建筑或者其他物体。relation是node和way之间的所有关系，我们这里没有用上所以不做介绍了。

 

2：分析地图数据，并根据坐标位置来生成最基本的建筑和道路的几何形状

地图数据分析到这就基本上差不多了，下面我们需要一种机制能够将地图的xml数据传入到houdini里面。houdini自身的程序语言Vex实现不了这个，但是好在houdini已经非常好的整合了Python进来，通过Python那问题就容易多了。

首先科普一下python查询xml的参考：

1：http://www.cnblogs.com/fnng/p/3581433.html

2：http://www.cnblogs.com/xuxm2007/archive/2011/01/16/1936610.html

 

针对于houdini里面建立python节点并抓取数据我放了源代码，这里我只讲一讲方法：

1：根据bounds做好边界限定（源数据有很多超出边界的杂点），并建立好缩放。

2：读取每个node标签，根据每一个node上面的经纬坐标增加一个对应的几何点。并把所有有效node的id存入一个数组中。

3：读取每个way标签，遍历每个way下面的所有ref中提到的id，判定该id是否为有效id，如果是调出该id对应的点，在点上增加vertex并用vertex生成多边形线。

4：根据way里面对该对象描述的不同，确定是路还是建筑，是封闭多边形还是开放多边形。

5：如果是路就是开放多边形，并增加宽度给每个点；如果是建筑则是封闭多边形，这个多边形则是之后挤出建筑的底座形状。

源码：

# This code is called when instances of this SOP cook.
geo = hou.pwd().geometry()

# Add code to modify the contents of geo.

def mapCreate():
    import xml.dom.minidom

    mapScale = hou.pwd().parm("map_scale").eval()

    mapFile = hou.pwd().parm("map_data").evalAsString()
    map = xml.dom.minidom.parse(mapFile)

    #test whether the file path is legal
    if not (mapFile.endswith((".osm",".xml"))):
        print "This node just except xml or osm file."
        return

    mapIDs = []
    mapPoints = []

    streetWithScale = 0.0001 * hou.pwd().parm("street_width_scale").eval()

    polyAttrib  = geo.addAttrib(hou.attribType.Prim, "type", 0)
    streetWidth = geo.addAttrib(hou.attribType.Prim, "width", 0.0)


    #define the boundary of this xml file
    def getBoundary(mapData, mapScale):
        boundary = mapData.getElementsByTagName("bounds")[0]
        minLat = float(boundary.getAttribute("minlat"))
        maxLat = float(boundary.getAttribute("maxlat"))
        minLon = float(boundary.getAttribute("minlon"))
        maxLon = float(boundary.getAttribute("maxlon"))

        #return [scale, minlon, minlat, maxlon, maxlat]
        return [(mapScale / max((maxLon - minLon), (maxLat - minLat))), minLon, minLat, maxLon, maxLat]

    #draw points
    def drawNodes(mapNode):
        lontitude = float(mapNode.getAttribute("lon"))
        latitude  = float(mapNode.getAttribute("lat"))

        scaledLontitude = boundary[0] * (lontitude - 0.5 *(boundary[1] + boundary[3]))
        scaledLatitude  = boundary[0] * (latitude  - 0.5 *(boundary[2] + boundary[4]))
        if (boundary[1] <= lontitude <= boundary[3] and boundary[2] <= latitude <= boundary[4]):
            point = geo.createPoint()
            point.setPosition((scaledLontitude,0,-scaledLatitude))

            pointID = int(mapNode.getAttribute("id"))
            #this two lists are prepared for create ways
            mapIDs.append(pointID)
            mapPoints.append(point)

    #create ways no matter is closed or open loop
    def createWays(mapWays):
        
        #loop for each way tag
        for mapWay in mapWays:
            nodes = mapWay.getElementsByTagName("nd")
            tags  = mapWay.getElementsByTagName("tag")

            #put all points of one way into the wayPoint
            wayPoints = []
            #flag the define this "way" is closed
            isClosed = 0

            #flag for different elements
            isHighway   = 0
            isBuilding  = 0
            #these way is the "v" value of "highway"
            isFirstWay  = 0
            isSecondWay = 0
            isThirdWay  = 0
            isLastWay   = 0


            #define the way should be closed or open, and get different way types
            #the feature list of every elements in this xml map could be found in this link:
            #http://wiki.openstreetmap.org/wiki/Map_Features
            for tag in tags:
                keyType = str(tag.getAttribute("k"))
                if keyType == "building" or keyType == "craft":
                    isBuilding = 1
                    isClosed = 1
                elif keyType == "highway":
                    isHighway = 1
                    valueType = str(tag.getAttribute("v"))
                    if valueType == "motoway" or valueType == "trunk":
                        isFirstWay = 1
                    elif valueType == "primary" or valueType == "secondary" or valueType == "tertiary":
                        isSecondWay = 1
                    elif valueType == "unclassified" or valueType == "residential" or valueType == "service":
                        isThirdWay = 1
                    else:
                        isLastWay =1


            if isBuilding == 1 or isHighway ==1:
                #find every nodes in each way
                for node in nodes:
                    ref = int(node.getAttribute("ref"))
    
                    #test wether this point could be found in the point list we created
                    #if we can find it, get the index of this point in this list,
                    #we will use this index to get the actual point
                    try:
                        index = mapIDs.index(ref)
                    except:
                        index = -1
    
                    if(index != -1):
                        point = mapPoints[index]
                        wayPoints.append(point)

                        
                #create the polygon using the way point list
                poly = geo.createPolygon()
                
                for point in wayPoints:
                    poly.addVertex(point)
                poly.setIsClosed(isClosed)

                #type : 0 - building ; 1 - first_way ; 2 - second_way ; 3 - last_way
                if isBuilding :
                    poly.setAttribValue(polyAttrib, 0)
                elif isHighway :
                    if isFirstWay:
                        poly.setAttribValue(polyAttrib, 1)
                        streetWidth_temp = streetWithScale * boundary[0] * 1
                        poly.setAttribValue(streetWidth, streetWidth_temp) 
                    elif isSecondWay: 
                        poly.setAttribValue(polyAttrib, 2) 
                        streetWidth_temp = streetWithScale * boundary[0] * 0.6
                        poly.setAttribValue(streetWidth, streetWidth_temp) 
                    elif isThirdWay: 
                        poly.setAttribValue(polyAttrib, 3) 
                        streetWidth_temp = streetWithScale * boundary[0] * 0.3
                        poly.setAttribValue(streetWidth, streetWidth_temp) 
                    else: 
                        poly.setAttribValue(polyAttrib, 4)
                        streetWidth_temp = streetWithScale * boundary[0] * 0.2
                        poly.setAttribValue(streetWidth, streetWidth_temp)


    boundary = getBoundary(map, mapScale)

    mapNodes = map.getElementsByTagName("node")
    mapWays  = map.getElementsByTagName("way")

    for mapNode in mapNodes:
        drawNodes(mapNode)

    createWays(mapWays)

mapCreate()

 

3：根据道路和现有建筑布局，用点填充其他空余空间

将地图数据提取到houdini之后，发现还有很多道路之间本来有建筑的地方是空的，主要是因为该地图的数据对于建筑还是远远不够，尤其是美国以外的国家，基本上除了大条的路其他的信息都非常不全。免费的质量就是这样没办法。

接下来就是在Houdini里面的处理了。

整体结构：

根据线生成街道：

生成原理街道的建筑点：

删除距离小于一定阈值的点：

根据每个点相互间的距离以及其到街道的距离判定每个建筑的最大宽度：

4：程序化生产成本每个不一样的建筑（把整个已经打成了一个包）：

这个自己做的节点生成的建筑其实还是有很多问题的，没有做更深入的控制了，先就做成这样了。

基本上思路就是这样了，中间遇到不不少困难，第一个就是在houdini里面Python写起来一点都没感觉，主要是因为用Python调houdini实例的方法还不熟悉，刚接触难免会这样，第二个是整个过程中每一个都步奏都是自己接触的少的，每一个环节在开始做之前需要不停的查找文献，尤其是在siggraph上发表的一些学术文章，现在算是真正体会了别人technical director平时都是在干嘛了，寻找发现方法，用艺术家的角度将方法实现，并确保别人也能够使用。我在这方面还是差太远了，总是觉得怎么做都还是不够，这种不够不是说没做好还是怎样，只是像一同水，不论怎么加水都还是看不见满，也许明年去base学习后会找到答案吧。

这是回国实习前最后做的一个东西了，下一步就是回国适应工作环境，开始实战。