Table of contents

Shape models from PDS

NASA's Planetary Data System (PDS) is an archive of data from planetary missions. There are many shape models available, including some Saturnian Satellites, asteroids, and comet nuclei. A list of available shape models can be found at the PDS Small Bodies Node.

As an example, I will use one of the shape models made for comet 67P/Churyumov-Gerasimenko based on data from the Rosetta mission. The shape model is the 788,000 facet file in version 2.0 of the data set. The model was created by astronomers at the Planetary Science Institute and the Laboratoire d'Astrophysique de Marseille: cg_spc_shap5_788k_cart.wrl. PDS data all have associated labels that contain vital metadata. For our file, it is cg_spc_shap5_788k_cart.lbl.

The file is formatted using the Virtual Reality Modeling Language (VRML). This format is compatible with PDS standards. The top of the file has some VRML markup and comments, followed by a series of vertices (x, y, z) then facet definitions (vertex1, vertex2, vertex3, vrml flag). The VRML flag indicates the end of the each facet definition. The PDS shape models I've worked with assume triangles, so this is ignored.

We programmatically read in the data with guidance from the PDS label. The label is based on PDS version 3 standards. I wrote a basic PDS3 reader and use it in this example. It reads the ODL-formatted PDS3 labels into dictionary-like objects.

import pds3
lbl = pds3.read_label('cg_spc_shap5_788k_cart.lbl')

The label consists of three table objects: the "counts" table, the vertex table, and the plate (or facet) table. The counts table describes the lengths of the vertex and plate tables, but is not strictly needed since the number of vertices and plates are also specified in their respective tables. We can loop over all the vertices and facets with the following code:

n_vertices = np.empty((lbl['VERTEX_TABLE']['ROWS'], 3))
n_facets = np.empty((lbl['PLATE_TABLE']['ROWS'], 3), int)

wrl = open('cg_spc_shap5_788k_cart.wrl')
f.seek(lbl['RECORD_BYTES'] * (lbl['^VERTEX_TABLE'][1] - 1))
for i in range(n_vertices):
    x, y, z = wrl.readline().split()

f.seek(lbl['RECORD_BYTES'] * (lbl['^PLATE_TABLE'][1] - 1))
for i in range(n_facets):
    vertices = wrl.readline().split()[:-1]

wrl.close()

Note the vertices line skips the last (VRML flag) column (

[:-1]
). The above code works with the 67P/Churyumov-Gerasimenko and 9P/Tempel 1 shape models.

Format as POVRay mesh

POVRay has a similarly formatted object called

mesh2
. It consists of a list of vertices and facets. With a few more lines of code, we can save a POVRay-formatted comet nucleus with the object name "CG":

import numpy as np
import pds3

lbl = pds3.read_label('cg_spc_shap5_788k_cart.lbl')
n_vertices = lbl['VERTEX_TABLE']['ROWS']
n_facets = lbl['PLATE_TABLE']['ROWS']

pov = open('cg_spc_shap5_788k_cart.inc', 'w')
wrl = open('cg_spc_shap5_788k_cart.wrl')

pov.write('''#declare CG = mesh2 {{
  vertex_vectors {{
    {}'''.format(n_vertices))

wrl.seek(lbl['RECORD_BYTES'] * (lbl['^VERTEX_TABLE'][1] - 1))
for i in range(n_vertices):
    v = wrl.readline().split()
    pov.write(',\n    <{}>'.format(','.join(v)))

pov.write('''  }}
  face_indices {{
    {}'''.format(n_facets))

wrl.seek(lbl['RECORD_BYTES'] * (lbl['^PLATE_TABLE'][1] - 1))
for i in range(n_facets):
    vertices = wrl.readline().split()[:-1]
    pov.write(',\n    <{}>'.format(','.join(vertices)))

pov.write('''
  }
}
''')

wrl.close()
pov.close()

Running this should create a file named cg_spc_shap5_788k_cart.inc:

#declare CG = mesh2 {
  vertex_vectors {
    393963,
    <0.57083178E+00,-0.10004439E+01,0.53256863E+00>,
    <0.56595045E+00,-0.10033593E+01,0.52524394E+00>,
    <0.55785316E+00,-0.99786305E+00,0.52076197E+00>,
...
    <-0.52038443E+00,-0.10123563E+01,-0.92725754E+00>,
    <-0.53658319E+00,-0.10089127E+01,-0.93165350E+00>,
    <-0.60233432E+00,-0.10027351E+01,-0.94259828E+00>  }
  face_indices {
    787920,
    <427,624,275>,
    <1669,1840,1420>,
    <4654,4772,4344>,
...
    <336690,393550,336861>,
    <336861,393879,336565>,
    <393879,280730,336565>
  }
}

Absolute orientation and size

The 67P nucleus shape models have units of kilometers. The object is in an excited rotation state, so the pole is not fixed in space, but is close to the shape model's +z axis. See the shape model documentation for details.

To use the shape model in a POVRay scene, include our new file and declare an object named CG. Create a file named cg.pov and insert the following to include the default POVRay colors, our shape model, and make the background black:

#include "colors.inc"
#include "cg_spc_shap5_788k_cart.inc"
background { black }

Next, create the camera:

camera {
  orthographic
  location <0,0,0>
  look_at  <-30,0,0>
  right <-1.33,0,0>
  translate <30,0,0>
  rotate <90,0,0>
  angle 28
}

I am using an orthographic projection to simulate the view from a distant telescope. I also setup the camera assuming the rendered image will have an aspect ratio of 1.33, and put the scene into a right-handed coordinate system with the

right <-1.33,0,0>
statement.
translate <30,0,0>
and
rotate <90,0,0>
places the camera long the +x-axis, the +y-axis to the right, and the +z-axis to the top.

Finally, add the comet, make shadowed areas black, and illuminated areas gray:

object {
  CG
  pigment { rgb 0.5 }
  finish {
    ambient 0
  }
}

Illumination

Rendering

Animation

Feedback

Questions? Comments? Suggestions? Contact me via e-mail.