VTK Python Example: Looking at satellite data

For this example we are going to be looking at the data from the SCIVIS Contest 2014. All data can be downloaded from the following url: http://www.viscontest.rwth-aachen.de/download.html

Test that vtk is correctly installed



In [ ]:

    
import vtk

1: Loading the data

All the data comes in vtk format, therefore it is straightforward to read it in vtk.



In [1]:

    
import os
os.chdir(r"D:\viscontest\viscontest14-data\mipas")
reader = vtk.vtkPolyDataReader()
reader.SetFileName("2011-MIPAS-VisContest.vtk")
reader.Update()
mipas=reader.GetOutput()

This contains the data for several orbits of the satellite, we are going to extract just the first orbit to makes thing easier



In [2]:

    
c1 = mipas.GetCell(0)
pids = c1.GetPointIds()
nps = pids.GetNumberOfIds()
# we are going to store the points in this new array
points = vtk.vtkPoints()
points.SetNumberOfPoints(nps)
for i in xrange(nps):
    pid = pids.GetId(i)
    p = mipas.GetPoint(pid)
    points.SetPoint(i,(p[0],p[1],20))  #  20 will be a fixed height for now

# and we are going to encapsulate them in this polydata
pd = vtk.vtkPolyData()
pd.SetPoints(points)

Now lets create a filter to make it easier to visualize the points



In [3]:

    
vfilter = vtk.vtkVertexGlyphFilter()
vfilter.SetInputData(pd)
vfilter.Update()

Lets setup the vtk environment



In [4]:

    
ren_win = vtk.vtkRenderWindow()
ren = vtk.vtkRenderer()
iren = vtk.vtkRenderWindowInteractor()
ren_win.SetInteractor(iren)
ren_win.AddRenderer(ren)
iren.SetInteractorStyle(vtk.vtkInteractorStyleTrackballCamera())
ren.GradientBackgroundOn()
ren.SetBackground2((0.5, 0.5, 0.5))
ren.SetBackground((0.2, 0.2, 0.2))
iren.Initialize()

Finally, lets add the points to the viewer



In [5]:

    
points_mapper = vtk.vtkPolyDataMapper()
points_mapper.SetInputConnection(vfilter.GetOutputPort())
points_ac = vtk.vtkActor()
points_ac.SetMapper(points_mapper)

#make points bigger
p = points_ac.GetProperty()
p.SetPointSize(2)
ren.AddActor(points_ac)

The vtk window should now be open but frozen, with the following command we can make it interactive. To exit interactive mode, press q while having focus on the vtk window



In [6]:

    
iren.Start()
#press q to stop interaction

You should see the following image in the 3d viewer

2: Adding some context

The above image is hard to interpret, lets add the map of the earth to make it easier to understand



In [7]:

    
os.chdir( r"D:\viscontest\viscontest14-data\support")
reader = vtk.vtkPolyDataReader()
reader.SetFileName("coastlines.vtk")
reader.Update()
cl = reader.GetOutput()
cl_mapper = vtk.vtkPolyDataMapper()
cl_mapper.SetInputData(cl)
cl_actor = vtk.vtkActor()
cl_actor.SetMapper(cl_mapper)
ren.AddActor(cl_actor)
ren_win.Render()

Notice in this case we don't need the vertex filter because the polydata is made of lines (instead of points). Now lets start the interactive window again



In [8]:

    
iren.Start()
#remember: q to stop interaction

You should see the following image in the 3d viewer

3: Mapping the samples attitude

Now we are going to map the attitude of the points over the map using the attitude data from the satellites. First lets look at this data



In [9]:

    
pdata = mipas.GetPointData()
print pdata









    



vtkPointData (000000000417BDD0)
  Debug: Off
  Modified Time: 215
  Reference Count: 2
  Registered Events: (none)
  Number Of Arrays: 5
  Array 0 name = time
  Array 1 name = altitude
  Array 2 name = orbit_id
  Array 3 name = profile_id
  Array 4 name = detection
  Number Of Components: 5
  Number Of Tuples: 1295837
  Copy Tuple Flags: ( 1 1 1 1 1 0 1 1 )
  Interpolate Flags: ( 1 1 1 1 1 0 0 1 )
  Pass Through Flags: ( 1 1 1 1 1 1 1 1 )
  Scalars: (none)
  Vectors: (none)
  Normals: (none)
  TCoords: (none)
  Tensors: (none)
  GlobalIds: (none)
  PedigreeIds: (none)
  EdgeFlag: (none)

We can see the altitude is stored in Array 1, now lets find its range



In [10]:

    
altitude = pdata.GetArray(1)
print altitude.GetRange()









    



(4.515999794006348, 30.0)

In the above drawing we used a constant height of 20, so this range looks appropriate. Now lets recreate the points polydata using this new information



In [11]:

    
for i in xrange(nps):
    pid = pids.GetId(i)
    alt = altitude.GetValue(pid)
    p = mipas.GetPoint(pid)
    points.SetPoint(i,(p[0],p[1],alt))
pd = vtk.vtkPolyData()
pd.SetPoints(points)
vfilter.SetInputData(pd)

And finally let start the viewer again in order to see the changes



In [12]:

    
ren_win.Render()
iren.Start()
# remember, q to stop

You should see the following image in the 3d viewer

4: Scalar data

Now we are going to look at some other scalars included in the data.

4.1: Time data

Lets start by looking at the time to get a better feeling of the direction in which the satellite travels



In [13]:

    
# array to use inside our new polydata
time_array = vtk.vtkFloatArray()
time_array.SetNumberOfComponents(1)
time_array.SetNumberOfTuples(nps)

#original time array
time_array_big = pdata.GetArray(0)

for i in xrange(nps):
    pid = pids.GetId(i)
    t = time_array_big.GetValue(pid)
    time_array.SetComponent(i,0,t)

time_array.SetName("time")
pdata2 = pd.GetPointData()
pdata2.SetScalars(time_array)
print pd









    



vtkPolyData (00000000042E57C0)
  Debug: Off
  Modified Time: 79460
  Reference Count: 3
  Registered Events: (none)
  Information: 00000000042EB8E0
  Data Released: False
  Global Release Data: Off
  UpdateTime: 79049
  Field Data:
    Debug: Off
    Modified Time: 78913
    Reference Count: 1
    Registered Events: (none)
    Number Of Arrays: 0
    Number Of Components: 0
    Number Of Tuples: 0
  Number Of Points: 1326
  Number Of Cells: 0
  Cell Data:
    Debug: Off
    Modified Time: 78916
    Reference Count: 1
    Registered Events: (none)
    Number Of Arrays: 0
    Number Of Components: 0
    Number Of Tuples: 0
    Copy Tuple Flags: ( 1 1 1 1 1 0 1 1 )
    Interpolate Flags: ( 1 1 1 1 1 0 0 1 )
    Pass Through Flags: ( 1 1 1 1 1 1 1 1 )
    Scalars: (none)
    Vectors: (none)
    Normals: (none)
    TCoords: (none)
    Tensors: (none)
    GlobalIds: (none)
    PedigreeIds: (none)
    EdgeFlag: (none)
  Point Data:
    Debug: Off
    Modified Time: 79460
    Reference Count: 2
    Registered Events: (none)
    Number Of Arrays: 1
    Array 0 name = time
    Number Of Components: 1
    Number Of Tuples: 1326
    Copy Tuple Flags: ( 1 1 1 1 1 0 1 1 )
    Interpolate Flags: ( 1 1 1 1 1 0 0 1 )
    Pass Through Flags: ( 1 1 1 1 1 1 1 1 )
    Scalars: 
      Debug: Off
      Modified Time: 79457
      Reference Count: 2
      Registered Events: (none)
      Name: time
      Data type: float
      Size: 1326
      MaxId: 1325
      NumberOfComponents: 1
      Information: 0000000000000000
      Name: time
      Number Of Components: 1
      Number Of Tuples: 1326
      Size: 1326
      MaxId: 1325
      LookupTable: (none)
      Array: 00000000042C28E0
    Vectors: (none)
    Normals: (none)
    TCoords: (none)
    Tensors: (none)
    GlobalIds: (none)
    PedigreeIds: (none)
    EdgeFlag: (none)
  Bounds: 
    Xmin,Xmax: (-148.527, 145.782)
    Ymin,Ymax: (-87.2112, 89.4011)
    Zmin,Zmax: (4.878, 29.963)
  Compute Time: 79462
  Number Of Points: 1326
  Point Coordinates: 0000000002CB3560
  Locator: 0000000000000000
  Number Of Vertices: 0
  Number Of Lines: 0
  Number Of Polygons: 0
  Number Of Triangle Strips: 0
  Number Of Pieces: 1
  Piece: 0
  Ghost Level: 0

We can see the time array is now located inside the point data. Now lets update the viewer and look at the changes



In [14]:

    
vfilter.SetInputData(pd)



In [15]:

    
ren_win.Render()
iren.Start()
# remember, q to stop

You should see the following image in the 3d viewer

In order to look properly at scalar data we need a look up table. Lets create one



In [16]:

    
t0,tf = time_array.GetRange()
lut = vtk.vtkColorTransferFunction()
lut.AddRGBPoint(t0 ,0,1,0) # green
lut.AddRGBPoint(tf ,1,0,0) # red
points_mapper.SetLookupTable(lut)

And now lets look at the changes



In [17]:

    
ren_win.Render()
iren.Start()
# remember, q to stop

You should see the following image in the 3d viewer

4.2: Detection data

Now lets look at the more interesting data about the detections of the satellite. From the web page we see the posible values are as follows:

0: A value of 0 designates a clear air measurement, i.e. the sample detected nothing. These measurements are included for providing context only, i.e. in order to show the satellite's track.
1: The least significant bit is used for ice detections.
2: The second bit stands for ash detections.
4: Finally, the third bit designates sulfate aerosol.

As before we start by loading the scalars



In [18]:

    
detection_array = vtk.vtkFloatArray()
detection_array.SetNumberOfComponents(1)
detection_array.SetNumberOfTuples(nps)
detection_array_big = pdata.GetArray(4)

for i in xrange(nps):
    pid = pids.GetId(i)
    d = detection_array_big.GetValue(pid)
    detection_array.SetComponent(i,0,d)

detection_array.SetName("detection")
pdata2 = pd.GetPointData()
pdata2.SetScalars(detection_array)
vfilter.SetInputData(pd)

Now lets create a more appropriate lookuptable



In [19]:

    
lut = vtk.vtkColorTransferFunction()
lut.AddRGBPoint(0 ,0,0.2,0) # green
lut.AddRGBPoint(1 ,0,0,0.9) # blue
lut.AddRGBPoint(2 ,0,0.9,0.9) # cyan
lut.AddRGBPoint(4 ,0.9,0,0) # red
points_mapper.SetLookupTable(lut)
points_mapper.UseLookupTableScalarRangeOn()

And now lets look at the changes



In [20]:

    
ren_win.Render()
iren.Start()
# remember, q to stop

You should see the following image in the 3d viewer

5: Multple orbits

Now we are going to repeat what we did above for a random set of 20 orbits.



In [21]:

    
import random
norbits = mipas.GetNumberOfCells()
sample = [random.randrange(0,norbits) for _ in xrange(20)]

#remove original orbit
ren.RemoveActor(points_ac)
actors_dict = {}

#We will reuse the lut from the last section
lut = vtk.vtkColorTransferFunction()
lut.AddRGBPoint(0 ,0,0.2,0) # green
lut.AddRGBPoint(1 ,0,0,0.9) # blue
lut.AddRGBPoint(2 ,0,0.9,0.9) # cyan
lut.AddRGBPoint(4 ,0.9,0,0) # red

for orb in sample:
    c1 = mipas.GetCell(orb)
    pids = c1.GetPointIds()
    nps = pids.GetNumberOfIds()
    points = vtk.vtkPoints()
    points.SetNumberOfPoints(nps)
    for i in xrange(nps):
        pid = pids.GetId(i)
        alt = altitude.GetValue(pid)
        p = mipas.GetPoint(pid)
        points.SetPoint(i,(p[0],p[1],alt))
    
    pd = vtk.vtkPolyData()
    pd.SetPoints(points)
    detection_array = vtk.vtkFloatArray()
    detection_array.SetNumberOfComponents(1)
    detection_array.SetNumberOfTuples(nps)

    for i in xrange(nps):
        pid = pids.GetId(i)
        d = detection_array_big.GetValue(pid)
        detection_array.SetComponent(i,0,d)

    detection_array.SetName("detection")
    pdata2 = pd.GetPointData()
    pdata2.SetScalars(detection_array)

    vfilter = vtk.vtkVertexGlyphFilter()
    vfilter.SetInputData(pd)
    #add points to viewer
    points_mapper = vtk.vtkPolyDataMapper()
    points_mapper.SetInputConnection(vfilter.GetOutputPort())
    points_mapper.SetLookupTable(lut)
    # to keep the mapper from adjusting the range of the lut
    points_mapper.UseLookupTableScalarRangeOn()
    points_ac = vtk.vtkActor()
    points_ac.SetMapper(points_mapper)

    #make points bigger
    p = points_ac.GetProperty()
    p.SetPointSize(2)
    ren.AddActor(points_ac)
    actors_dict[points_ac] = orb

When showing multiple datasets it is important that the lookuptables for all of them matches. Here we are looking at 20 orbits with altitude and detection scalars. Lets look at the results



In [22]:

    
ren_win.Render()
iren.Start()
# remember, q to stop

You should see the following image in the 3d viewer

6: Some interaction

It is easy to get lost when there is so much data in the screen. Lets add some interaction in order to help fix this. We are going to make a small example by adding tooltips when you hover the mouse over an orbit.



In [23]:

    
balloon = vtk.vtkBalloonWidget()
balloon.SetInteractor(iren)
balloon.On()
for ac,orb in actors_dict.iteritems():
    message = "orbit number: %d"%orb
    balloon.AddBalloon(ac,message)

Lets test this, after starting the interaction leave the mouse over one of the points.



In [24]:

    
ren_win.Render()
iren.Start()
# remember, q to stop

You should see the following image in the 3d viewer

7: What's next?

We have seen the basics of vtk. What follows is a list of possible improvements to the visualization you could try:

Add more complete messages to the tooltips
Show time and detection information at the same time
You can use vtkGlyph3D insted of GlyphFilter to get more complex representations for the points
You could embed the application inside a GUI (GTK, QT, TK) and add controls for all the operations
You could map the data to a sphere
You could filter the data and only show certain times or data with detections
You could add the rest of the data from the contest

Feel free to try any other ideas you have. VTK is a very flexible and powerful platform. The documentation can be found here; and the wiki contains several useful examples



In [24]: