将netCDF雷达数据转化为VTK格式进行3D可视化

导语：之所以写这些实在是因为太冷门，国内能找到的资料很少，再加上我是对这些概念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()

如果缺引用看提示引入就好。我这里是生产的文件在ParaView中打开，也可以直接加上代码预览。

这是单层的毛坯效果，因为是点云，效果不太好，美化方法我也还在找，希望有人知道的指导一下，哈哈。

另外学这些新东西，特别是还冷门的一定要看官方文档！自己乱查效果很差，我自己就是浪费了两周才看官方文档做出来的。