python 科学计算三维可视化笔记（第四周 进阶实战）

---

内容来自中国大学MOOC，北京理工大学，python数据分析与展示课程，侵删。
如有错误，烦请指出。

---

python 科学计算三维可视化笔记 第四周 进阶实战

一、实例1：Dragon 绘制

• 读取 tar 压缩文件
• 读取 .ply 数据文件
• 渲染 dragon ply文件，绘制数据的 surface
• 显示可交互的结果
• 删除解压的文件夹

'''实例1：dragon'''
from mayavi import mlab
from os.path import join
import tarfile
 
# 读取tar压缩文件
dragon_tar_file = tarfile.open('dragon.tar.gz')
try:
    os.mkdir('dragon_data')
except:
    pass
dragon_tar_file.extractall('dragon_data')
dragon_tar_file.close()

# 读取数据文件
dragon_ply_file = join('dragon_data', 'dragon_recon', 'dragon_vrip.ply')
 
# 渲染dragon ply文件
mlab.pipeline.surface(mlab.pipeline.open(dragon_ply_file))

# 显示可交互的结果
mlab.show()
 
# 删除解压的文件夹
import shutil
shutil.rmtree('dragon_data')

二、实例2：Canyon 地形可视化

• 读取压缩文件
• 处理地形数据
• 渲染地形数据
• 清空内存
• 创建交互式的可视化窗口
• 显示可交互的结果

'''实例2：Canyon 地形可视化'''
import zipfile
import numpy as np
from mayavi import mlab
 
# 读取压缩文件
hgt = zipfile.ZipFile('N36W113.hgt.zip').read('N36W113.hgt')

# 处理地形数据
data = np.fromstring(hgt,'>i2')
data.shape = (3601, 3601)
data = data.astype(np.float32)
data = data[:1000, 900:1900]
data[data == -32768] = data[data > 0].min()
 
# 渲染地形数据
mlab.figure(size=(400, 320), bgcolor=(0.16, 0.28, 0.46))
mlab.surf(data, colormap='gist_earth', warp_scale=0.2,
            vmin=1200, vmax=1610)
 
# 清空内存
del data

# 创建交互式的可视化窗口
mlab.view(-5.9, 83, 570, [5.3, 20, 238])

# 显示可交互的结果
mlab.show()

三、实例3：Earth Graph

1. 数据处理

• 构造城市互联数据 connec_data
• 构造城市坐标数据 cities_data
• 建立城市互联列表 connec
• 建立城市索引 cities
• 建立城市:经纬度索引 cities_coord
• 建立经纬度列表 coords

2. 坐标转换

将经纬度的位置转换为三维坐标，lng 为经度，lat 为纬度：

• $x=\cos (lng)\ast \cos (lat)$
• $y=\cos (lng)\ast \sin (lat)$
• $z=\sin (lng)$

3. 绘制图像

• 建立窗口
• 绘制地球
• 绘制城市点位
• 显示城市互联
• 显示城市名字
• 绘制大洲边界
• 绘制赤道、南北回归线、南北极圈
• 显示可交互窗口

4. 代码

'''实例3：Earth Graph'''

import csv
import numpy as np
from mayavi import mlab
from mayavi.sources.builtin_surface import BuiltinSurface

########### 数据处理 ###########
# 城市互联数据
connec_data = """
Beijing,Shanghai
Bombay,Atlanta
Bombay,Newark
Bombay,New York
Chicago,Delhi
Doha,Houston
Dubai,Houston
Dubai,Los Angeles
Dubai,Sao Paulo
Dubai,San Francisco
Johannesburg,Atlanta
Los Angeles,Sydney
New York,Hong Kong
Newark,Hong Kong
Shanghai,Hong Kong
Toronto,Hong Kong
Vancouver,Hong Kong
"""

# 城市坐标数据，城市city，经度longtitude，纬度latitude
cities_data = """
Atlanta,       -84.42,  33.76
Beijing,       116.23,  39.54
Bombay,        72.82,   18.96
Chicago,       -87.68,  41.84
Delhi,         77.21,   28.67
Doha,          51.53,   25.29
Dubai,         55.33,   25.27
Hong Kong,     114.19,  22.38
Houston,       -95.39,  29.77
Johannesburg,  28.04,   -26.19
Los Angeles,   -118.41, 34.11
Newark,        -82.43,  40.04
New York,      -73.94,  40.67
San Francisco, -122.45, 37.77
Sao Paulo,     -46.63,  -23.53
Shanghai,      121.52,  30.91
Sydney,        151.21,  -33.87
Toronto,       -79.38,  43.65
Vancouver,     -123.13, 49.28
"""

# 城市互联列表connec
connec = list()
for city1, city2 in list(csv.reader(routes_data.split('\n')))[1:-1]:
    connec.append((cities[city1], cities[city2]))
#print('connec:', connec)
connec = np.array(connec)
#print('connec:', connec)

# 城市索引cities，城市:经纬度索引cities_coord，经纬度列表coords
cities = dict()
cities_coord = dict()
coords = list()
for line in list(csv.reader(cities_data.split('\n')))[1:-1]:
    name, lng, lat = line
    cities[name] = len(coords)    
    cities_coord[name] = (float(lng), float(lat))
    coords.append((float(lng), float(lat)))
#print('cities: ', cities)
#print('cities_coord: ', cities_coord)
#print('coords: ', coords)
coords = np.array(coords)
#print('coords: ', coords)


################# 坐标转换 #################
# 将经纬度的位置转换为三维坐标
lat, lng = coords.T * np.pi / 180
x = np.cos(lng) * np.cos(lat)
y = np.cos(lng) * np.sin(lat)
z = np.sin(lng)


################# 绘制图像 #################
# 建立窗口，背景颜色0.48，前景颜色0,0,0为黑色
mlab.figure(bgcolor=(0.48,0.48,0.48), fgcolor=(0,0,0), size=(400,400))
#mlab.figure(1, bgcolor=(0.48, 0.48, 0.48), fgcolor=(0,0,0),size=(400, 400))

# 绘制半透明球体表示地球
sphere = mlab.points3d(0, 0, 0, 
                       scale_factor=2,
                       color=(0.67,0.77,0.93),
                       resolution=50,
                       opacity=0.7,
                       name='Earth')
# 调整镜面反射参数
sphere.actor.property.specular = 0.45
sphere.actor.property.specular_power = 5
# 设置背面剔除，以更好的显示透明效果
sphere.actor.property.backface_culling = True


# 绘制城市点位
points = mlab.points3d(x, y, z, 
                       scale_mode='none',
                       scale_factor=0.03,
                       color=(0,0,1))
# 显示城市互联
points.mlab_source.dataset.lines = connec
points.mlab_source.reset()
mlab.pipeline.surface(points, color=(1,1,1),
                      representation='wireframe',
                      line_width=4,
                      name='Connections')
# 显示城市名字
for city, index in cities.items():
    label = mlab.text(x[index], y[index], city, z=z[index],
                      width=0.016*len(city), name=city)
    label.property.shadow = True


# 绘制大洲边界
continents_src = BuiltinSurface(source='earth', name='Continents')
# 设置LOD为2
continents_src.data_source.on_ratio = 2
continents = mlab.pipeline.surface(continents_src, color=(0,0,0))

# 构造经度数组，从0到360(2pi)
theta = np.linspace(0, 2*np.pi, 100)
# 绘制赤道线0，南北回归线23.5，南北极圈66.5
for angle in (-23.5*np.pi/180, -66.5*np.pi/180, 0, 23.5*np.pi/180, 66.5*np.pi/180):
    x = np.cos(theta) * np.cos(angle)
    y = np.sin(theta) * np.cos(angle)
    z = np.ones_like(theta) * np.sin(angle)
    mlab.plot3d(x, y, z, 
                color=(1, 1, 1),
                opacity=0.2, tube_radius=None)

# 显示可交互窗口
mlab.view(100, 60, 4, [-0.05, 0, 0])
mlab.show()

5. 结果

结果如下，只有大概的大洲轮廓，一些岛屿好像绘制不出来：