导语:之所以写这些实在是因为太冷门,国内能找到的资料很少,再加上我是对这些概念0基础,所以像无头苍蝇一样查了很久才勉强能有一点效果,希望能帮到有需要的人。
项目目的:以台风的雷达数据为基础,生产点云进行3D可视化(效果较一般)。
源数据类型:netCDF,如下图,雷达数据为77层,每层1689 x 2520个坐标点。
转换这些数据我是使用Python写的,所以这里的例子就是py的,官方文档里支持C#、Java、Cpp和Py。
第一步,引入相关的包。
① 引入读取netCDF相关的库 pip install netCDF4
② 引入vtk相关的库 pip install vtk
代码
import netCDF4 as nc
import vtk as vtk
import vtkmodules.vtkInteractionStyle
# noinspection PyUnresolvedReferences
import vtkmodules.vtkRenderingOpenGL2
from vtkmodules.vtkCommonColor import (
vtkColorSeries,
vtkNamedColors
)
from vtkmodules.vtkCommonCore import (
vtkLookupTable,
vtkUnsignedCharArray
)
from vtkmodules.vtkCommonTransforms import vtkTransform
from vtkmodules.vtkFiltersCore import vtkElevationFilter
from vtkmodules.vtkFiltersGeneral import vtkTransformPolyDataFilter
from vtkmodules.vtkFiltersModeling import vtkBandedPolyDataContourFilter
from vtkmodules.vtkFiltersSources import (
vtkConeSource,
vtkCubeSource
)
from vtkmodules.vtkInteractionWidgets import vtkOrientationMarkerWidget
from vtkmodules.vtkRenderingAnnotation import (
vtkAnnotatedCubeActor,
vtkAxesActor
)
from vtkmodules.vtkRenderingCore import (
vtkActor,
vtkPolyDataMapper,
vtkPropAssembly,
vtkRenderWindow,
vtkRenderWindowInteractor,
vtkRenderer
)
def FileConvert():
try:
file_path = "F:/3D/SourceData/Data/wrfout_d03_2017-08-23_03" #file path
ds = nc.Dataset(file_path)
points = vtk.vtkPoints()
cells = vtk.vtkCellArray()
polydata = vtk.vtkPolyData()
obj = ds.variables["REFL_10CM"]#Rander Data Object
latObj = ds.variables["XLAT"] #lat
longObj = ds.variables["XLONG"] #long
latData = latObj[:][0]
longData = longObj[:][0]
nx,ny,nz = obj.shape[2],obj.shape[3],obj.shape[1]
#nx,ny,nz = 500,500,1
#nz = 1
pointsNum = nx*ny*nz
data = obj[:][0]
for lel in range(0,nz):
print(lel)
for width in range(0,nx):
for height in range(0,ny):
if data[lel][width][height] != -35 and data[lel][width][height] != 0:# Filter
id = points.InsertNextPoint(longData[width][height]*100,latData[width][height]*100,data[lel][width][height])#Add Point
cells.InsertNextCell(1)
cells.InsertCellPoint(id)
print("Point Done")
pointsNum = points.GetNumberOfPoints()
# Create the color map
colorLookupTable = vtkLookupTable()
colorLookupTable.SetTableRange(-50.0,100.0)
colorLookupTable.Build()
# Generate the colors for each point based on the color map
#colors = vtkNamedColors()
Colors = vtk.vtkUnsignedCharArray()
Colors.SetNumberOfComponents(3)
Colors.SetName('Colors')
for i in range(0, pointsNum):
p = 3 * [0.0]
points.GetPoint(i, p)
dcolor = 3 * [0.0]
colorLookupTable.GetColor(p[2], dcolor)
color = 3 * [0.0]
for j in range(0, 3):
color[j] = int(255.0 * dcolor[j])
Colors.InsertNextTuple(color)
polydata.SetPoints(points)
#polydata.SetPolys(cells)
polydata.SetVerts(cells)
polydata.GetPointData().SetScalars(Colors)
polydata.Modified()
writer = vtk.vtkXMLPolyDataWriter()
writer.SetFileName("polyData15.vtp")
writer.SetInputData(polydata)
writer.Write()
except(Exception, ArithmeticError) as e:
template = "An exception of type {0} occurred. Arguments:\n{1!r}"
message = template.format(type(e).__name__, e.args)
print (message)
else:
print("the code is no problem")
finally:
print("Done")
FileConvert()