Downloading Scientific Datasets into TOPs (Part 03)

We left off having acquired our data from NOAA’s servers. We now need to convert our downloaded data into geometry that we can put into a renderer, and make sure that our geometry is mapped correctly onto a sphere. Let’s start by converting data to curves.

Creating the Initial Geometry

What we acquired in the last post are what are called ESRI Shapefiles. While we could reinvent the wheel and write our own shapefile parser, I (and likely you) don’t have that kind of free time. Instead, let’s make use of the Python Shapefile Library, created by Joel Lawhead. Make sure that you install this and make it available to Houdini.

https://pypi.org/project/pyshp

After doing that, set up a ROP Geometry Output to save out the created geometry to bgeo.sc format. I decided in this one to do the actual geometry processing in an external SOP network, which is represented in the screenshot by the node titled “process_shapefile_geo”.

Make sure that the “Use External SOP” checkbox is enabled, and also make sure that your bgeo files will be saving out somewhere logical and with an appropriate filename. As you can see, I embedded a small expression to end each saved file according to the @pdg_index attribute, with three digits of numpadding. This will save out files numbered 001 to 366, one for each of the shapefiles extracted in the Python Processor.

Inside of the SOP network, all we need is a Python node and an out null for organization.

Breaking the code down into sections, we start by importing the python modules we need. We want os.path to build the search path for the shapefiles.

import os.path
import shapefile

If you forgot to install the shapefile library earlier, this node will immediately start throwing an error.

Next, we build the file path to the shapefiles. The files are named consistently and are placed in a folder that is defined on a parameter on the HDA. The defining difference between each file is the three-digit numpadding before the date range, so we fetch the pdg_index by creating a parameter on the python node that references @pdg_index, evaluate it as a string, and use the zfill method to pad. We then insert that string into the file name, and append the final filename onto the folder path.

root_dir = hou.evalParm("../../shapefile_directory")
filenum = str(node.evalParm("./pdg_index")).zfill(3)
filename = "median_extent_N_%s_1981-2010_polyline_v3.0.shp" % (filenum)
shape_path = os.path.join(root_dir, filename)

Next, we create a shapefile Reader object, providing the created path to the current shapefile. The reader object provides methods for interpreting the shapefile data. In this case, we want to specifically read shape data so we call the shapes method fetch that. With that shape data, we can extract the points. The points array holds a list of coordinates for each and every point. The return type of the shapes method is an array, however in this case there is only one element. We also want the parts which, in Houdini terms, could be considered analogous to the primitive definitions.

We’ll want to loop through the points and parts in order to build this data into Houdini points and prims, so we’ll store keep the lengths of the returned arrays on hand for easy access.

# Shapefile Data
sf = shapefile.Reader(shape_path)
shapes = sf.shapes()
points = shapes[0].points # This data set seems to only have one shape per file.
# List of the first point numbers of each independent segment
parts = shapes[0].parts 
numpt = len(points)
numpart = len(parts)
next_part = 0

Finally, we build the loops necessary to build the geo. The encasing loop runs per-part, creating polylines per part. We need to cut up the points array according to the point ranges of each part. The part array specifically holds the point index of the first point of each part, so we can use the part value and the value of the next part to slice the points array. We then run the sliced array through a for loop to create a point on the polyline for each item of the sliced point array.

for index, part in enumerate(parts):
    if (index < numpart - 1):
        next_part = parts[index + 1]
    else:
        next_part = numpt - 1
    # Create new polyline per part
    line = geo.createPolygon(is_closed=False)    
    part_points = points[part:next_part]
    for point in part_points:

        pos = (point[0], 0, point[1])
        point = geo.createPoint()
        point.setPosition(pos)
        line.addVertex(point)

And with that, we should have what we need to save out the curves representing the mean sea ice extent for each day of the year. At this point you can dirty and cook the TOP chain, saving out each shapefile as a set of curves in a sequence of bgeo.sc files.

Mapping our Geometry to a Sphere

At this point, we can now take our work into SOPs for final presentation. After pulling in the curves with a File node and adjusting your camera with space+f, you may come to the sudden realization that these curves are at extremely realistic scale, i.e. thousands of kilometers across. For my own sanity, I use a transform node here to scale everything down to a more comfortable size

In the pos_to_angle wrangle, we need to start converting our coordinates from a rectangular grid centered on the North Pole to a polar coordinate so that we can position them to a polar-unwrapped sphere. In this wrangle, we’re just acquiring the rotation component of the coordinate. We use dihedral to figure out how much we need to rotate away from the z-axis around y. This will give us values from -270 -> 90, which we can refit to the 0 to 1 range to map to the UV grid of the sphere.

// pos_to_angle

vector4 qangle = dihedral({0,0,1}, normalize(v@P));
vector angle = degrees(quaterniontoeuler(qangle, XFORM_XYZ));
float rotation = angle.y;
f@nrot = fit(rotation, -270, 90, 0, 1);

In the wrap_to_sphere wrangle, we take our curves and actually map them against the sphere. The sphere has been unwrapped using a UV Texture node set to Polar mode, and the uv coordinates have been scaled by -1 on the x axis and transformed by 1 on x on the UV Transform node.

In the wrangle itself, we need to start by figuring out the distance component of the polar coordinate. First, we want to grab the scaling factor we used in the transform1 null. Then, we take the circumference of the earth in meters and scale it by our scaling factor. The reason for doing this is that our sphere’s UV coordinates go from 0 to 1 from top to bottom of the sphere, so we can remap the length of our P coordinates of our curve points from an input range of 0 to half of our circumference (the distance from north pole to south pole over the surface) to 0 to 1.

Next, we need to make sure we have the w-coordinate of the sphere’s UVs. When you run the UV Texture node in Polar mode, the radius of the sphere gets placed in the third component of the UV vector. As much as I’d love to ignore it, we need this value for the uvsample() function to work. To fetch it, I simply promoted the UV to a detail attribute in Maximum mode. In the wrangle we can then fetch the value of the UV detail attribute and isolate the z/w coordinate.

Finally, we assemble the vector for the uvsample() using our direction, length, and radius values. We then set the P of each point of the curve to P location of the matching spot on the sphere’s UV map from the assembled uv coordinate.

// wrap_to_sphere

float scaling_factor = chf("scaling_factor");
float polar_circumference = 40007863 * scaling_factor;
float length = fit(length(v@P), 0, polar_circumference / 2, 1, 0);

vector sphere_uvw = detail(1, "uv", 0);
float wmax = sphere_uvw.z;

vector uvloc = set(f@nrot, length, wmax);
v@P = uvsample(1, "P", "uv", uvloc);

You can then finishing this off by matching the rotation of the curves to your texture with a Transform node, in this case I had to rotate by -135 around Y to match correctly. Besides that, just add whatever you want to make the curve stand out in viewport or render, and render it as a curve object in your renderer of choice.