Example: Create Terrain Elevation Raster from a NASA PDS Table

This is a long, multistep example showing many powerful features of Manifold in various steps.

We import terrain elevation raster data from a NASA PDS archive and style the images with an elevation palette and hill shading. We notice one of the images was wrongly georegistered by NASA, so we import the original LiDAR point data from a table in the PDS archive. We transform the imported table from 0 to 360 degree longitudes into +/-180 longitudes, construct a geometry field, and then we create a drawing from the table that shows the LiDAR point data in the table. Next, we use Kriging to create a terrain elevation raster from the LiDAR point drawing, and then we style that and hill shade it so it can be exactly overlaid on a Bing satellite layer to enhance the satellite photography with enhanced terrain relief.

This topic covers many steps, but an experienced Manifold user can do them in just a few minutes. The Create Terrain Elevation from a NASA PDS Table video version of this topic shows how in under six minutes.

Import and Use PDS Images

We will use data from the NASA PDS Geosciences Node page at https://pds-geosciences.wustl.edu/missions/earth_lidar/index.htm The data set contains Wallops/GSFC Airborne Topographic Mapper digital elevation data of explosion craters at the Nevada Test Site in the US. The large craters excavated by open air nuclear explosions at the Test Site provide pristine examples of craters created by very high energy events, and are useful examples for modeling hypervelocity impact cratering on Earth and other worlds.

We use a PDS archive provided as a zipfile at https://pds-geosciences.wustl.edu/earth/wff-e-atm-1_5-v1/atm_9001.zip

Unzipping the zipfile and then drilling down into the data folder as seen above, we find a typical group of PDS files, with images in .img files accompanied by .hdr header files, , tables in .tab files and metadata in .lbl files. This is an older PDS archive that does not use XML files for metadata and cataloging, as do more recent PDS file sets. See the PDS, XML topic for more info on PDS formats.

We import the .img files using File - Import and then double-clicking on the image file desired. We style them to use Channel 0 (these are single-channel terrain elevation files), using a palette and hill shading to generate hill shaded images like that seen below.

The image shows the 280 meter (920 feet) diameter crater created in hard, dry rock by the Schooner 30 kiloton nuclear explosive detonated at a depth of 111 meters (365 feet) on 8 December 1968, on Pahute Mesa in the Nevada Test Site.

We can use the hill shaded images to enhance satellite photography. We create a map using a Bing satellite image server layer as a base layer, and then we drag and drop our hill shaded imported images as layers into the map, using the Layers pane to set 50% opacity for the hill shaded imported images. This is the same technique shown in the Example: Enhance Terrain with Curvatures topic.

The illustration above shows the map panned and zoomed to the Schooner crater image, seen overlaid on a Bing satellite photo. Manifold uses the .hdr information when importing .img images accompanied by a header file to georeference the images, automatically assigning the correct coordinate system, so the imported images line up with "known good" background layers like Bing.

Zooming further into the Schooner crater, we can see how overlaying the partially-transparent hill shaded terrain elevation image brings out terrain relief not easily visible in the Bing satellite photo.

Without the partially-transparent schooner2x2 image layer, the Bing satellite photo does not show the topographic terrain undulations that are easily visible when the schooner2x2 layer is turned on.

The above illustration shows the Sedan crater, a massive, 410 meter (1345 feet) diameter nuclear crater, the largest explosion crater on Earth, created by the Sedan 104 kiloton nuclear detonated at a depth of 194 meters (635 feet) on 6 July, 1962. Other craters visible in the view are subsidence craters where the surface has collapsed into a large cavity created by an underground nuclear explosion.

The Nevada Test Site has hundreds of subsidence craters from the hundreds of underground tests conducted at the site, so many that some parts of the site look like the surface of the Moon, but craters excavated by tests designed to vent into the open air are rare. Relatively few tests deliberately created craters, to minimize the amount of fallout spewed into the atmosphere over the continental United States. See a short Department of Energy video of the Sedan explosion, providing a close-up view of dome formation, venting and the initial stages of material ejection.

Our third image example shows the swath of terrain elevation data acquired over the Buckboard Mesa section of the test site.

Zooming into the Buckboard Mesa image, we see that the PDS terrain elevation layer is offset from the Bing image. The other PDS images align very well with Bing, but in this case the craters in the PDS image are offset from the crater circles seen in Bing. Given ambiguities in the PDS specification, at times georeferencing of images can be off. The magenta arrow added to the illustration shows how the PDS image should be offset.

Import a PDS Table and Build a New Terrain Image

We could tinker with coordinate system settings for the PDS image in a trial and error way until we get it to line up with Bing. A more precise approach, though, would be to create a new terrain image from the bbmesa.tab table that is provided in the data set. Manifold in the hands of an experienced user can do that very rapidly, much more rapidly than messing about with trial and error.

Opening the bbmesa.tab file in Notepad we see it contains plain text listing the X (longitude), Y (latitude), and Z (height) values for the LiDAR data points acquired by the Airborne Topographic Mapper laser altimeter instrument. There are no field names at the beginning of the file, and the spaces between the values in each row are multiple space characters, and not tab characters as we might expect in a .tab file. In this case, the .tab extension is a mnemonic for table, and not for tab character.

Tables in PDS data sets can be huge. It is not a good idea to attempt to open very large .tab files in Notepad. Use a more capable text editor.

The format used in the .tab file is similar enough to CSV that we can get it into Manifold using the CSV data source dataport.

We choose File - Create - New Data Source from the main menu. The dropdown menu provides a list of favorites to choose from as well as a More... option.

We choose More... to launch the New Data Source dialog.

Using the New Data Source dialog gives us more options than a simple import. We choose File: csv in the Type box and then we click the [...] browse button to choose the file.

We navigate to the bbmesa.tab file, click it to choose it, and press Open.

We uncheck the First line contains field names box, and for the List delimiter we choose (space). There is more than one space between each value, but we can clean that up later after we get the table into our project.

We press Create Data Source.

A new data source appears in the project. We drill down into it, click onto the bbmesa table and press Ctrl-C to Copy the table. Next, we click outside of the data source to move the focus onto the main part of the project and we press Ctrl-V to Paste a copy of the table. That pastes the copy into the project, as opposed to pasting another copy of the table within the data source.

A new table appears in the main part of the project. We now have a copy of the table imported into the project. We can now delete the data source, since we do not need it. We Click on the data source to highlight it, and then we press the Delete button on the toolbar.

The data source disappears. Next, we double-click the bbmesa 2 table to open it.

The table appears with gray background indicating it is read-only. It is read-only because it does not yet have any index. It also has extra, blank, columns, created as result of the extra space characters between values in the original .tab text file. Those are easy to delete.

To clean up the table we choose Edit - Schema from the main menu, to launch the Schema dialog.

We Ctrl-click the columns we wish to delete. That selects them.

We press the Delete button on the toolbar to delete the selected columns.

Next, we rename the columns to something more mnemonic. We double-click each column name to rename them.

Manifold shows the new column names in provisional, blue color. Changes will not be committed to the table until we press the Save Changes button.

Next, we press the Add button to add a new field. This will be a new longitude field that will host GIS-style +/-180 degree longitude values instead of the 0 to 360 degree longitude values originally in the table.

We provide a name, Longitude, and choose float64 as the Type. We press OK.

We press the Add button again to add another new field. This will be a new geometry field that will host point geometry, from which we can create a drawing.

We provide a name, Geom and choose geom as the Type. We also change the coordinate system to Latitude / Longitude, since the table provides latitude and longitude coordinate values. We press OK.

The schema shows the changes made in blue, provisional color. We now need to add a key field and index to the table, to make it read/write. That is easy to do with one click.

We press the Add Identity button.

That adds a standard mfd_id identity field and mfd_id_x index to the table. We press Save Changes to commit the changes to the table and to exit the Schema dialog.

Much better! With a quick use of the Schema dialog we have cleaned up the table, provided sensible field names, added new fields to host data we will create, and have added a key field and index to make the table writable. Not bad for a few minutes with one dialog.

The new fields we created, Longitude and Geom have NULLs in them since we have yet to populate them. To do that, we use the Transform pane.

With the focus on the open bbmesa 2 table window, in the Transform pane we choose the Longitude field. Double-click on the Expression template to launch it.

In the Expression template, press the Edit Expression button to launch the expression builder dialog.

In the expression builder dialog, we write the expression

[Lon360] - 360

That is the simple expression, a hack, to convert the 0 to 360 values in the Lon360 field using NASA's system to +/-180 longitude values for the Western hemisphere, where the Nevada Test Site is located.

We can get away with such a simple expression for the Nevada test site, because we know all the values will end up being negative numbers. A more complicated expression would be required for a world-wide conversion.

Press OK to accept the expression and to return to the Transform pane.

The Transform pane parses the expression we wrote, and if it does not contain syntactic errors loads it into the Expression template. For the Result destination, we leave the default choice of Same Field, meaning we will save the results of the expression into the Longitude field, the same field chosen in the Field box.

If we would like to see what the expression does before putting the transform into action, we can press the Preview button to get a preview.

The preview display adds a virtual preview column in blue preview color, with the column titled the name of the transform template that is used. The destination column, the Longitude column, is marked with a dot icon to show where the previewed data will go. In the illustration above, we have dragged the preview column next to the Longitude column. This provides a useful check our expression is correct in this case, since we can tell by mousing over the map window that Buckboard Mesa is approximately at -116 longitude.

In the Transform pane, we press the Transform button, to put the transform into action.

Pressing the Transform button commits the change and populates the Longitude field with new values.

We now will populate the Geom field with geometry calculated from the other fields.

With the focus on the open bbmesa 2 table window, in the Transform pane we choose the Geom field. Double-click on the Compose template to launch it.

In the Compose template we choose point with z as the Compose option.

In the X box, from the pull down menu we choose the Longitude field.

In the Y box, from the pull down menu we choose the Latitude field.

In the Z box, from the pull down menu we choose the Z field.

For the Result destination, we choose the default Same Field meaning that we will update the Geom field that has been chosen in the Field box.

Press Transform.

The table now has a geometry field called Geom, which is populated with geometry for a point for each record. The points also include Z values. That is not required for the next steps, but it seems like good form with point data sets that represent X,Y,Z data.

Create a Drawing from the Table

In the Project pane we right-click on the bbmesa 2 table and choose Create - New Drawing from the context menu.

We choose a name for the drawing, change the coordinate system to Latitude / Longitude, and leave other defaults as is. We press Create Drawing.

A new drawing, called bbmesa Drawing, appears in the project. We can drag and drop it into our map.

The drawing shows the large number of LiDAR acquired X,Y,Z data points.

The LiDAR data set was published in Latitude / Longitude coordinate system, which uses degrees as a unit of measure, so when we created a drawing from the LiDAR points table it was created in Latitude / Longitude. We will take a moment to reproject it into Pseudo-Mercator, so we can use sensible units of measure, like meters, in subsequent steps.

With the focus on the bbmesa Drawing layer in the map, in the Info pane, click the coordinate picker button and choose Reproject Component.

In the Reproject Component dialog, click the coordinate system picker button, and then click the Pseudo-Mercator choice, one of the default Favorites in the resulting menu.

Back in the Reproject Component dialog, press Update Component. That will reproject the bbmesa Drawing layer into Pseudo-Mercator.

There will be no visible difference in the map, since the map window automatically reprojects whatever coordinate system the layer uses into the Pseudo-Mercator projection used by the map window.

Create a Terrain Elevation Image using Kriging

A mass of points in a vector drawing do not provide a very useful display, so next we will create a terrain elevation raster image using the Interpolate : Kriging transform operation.

With the focus on the map window, in the Transform pane we choose the bbmesa Drawing layer and the Geom field. We double-click the Interpolate template to launch it. In the illustrations below, we have resized the Transform pane to provide more room for the many controls used in Interpolation.

In the Interpolate template, from the pull down menu for Interpolation we choose Kriging as the interpolation method.

For the Z parameter, from the pull down menu we choose the Z field.

Fro the Resolution parameter, we use the default value of 1 which we can see from the Unit field will be a value in meters. This is one reason we reprojected the drawing into Pseudo-Mercator, so we could sensible (meters) units of measure, instead of having to specify a number like 0.00001 in degrees, that being approximately the size of one meter in degrees at the latitude of the area of interest.

For the number of Neighbors we choose 0. , which means to use all Voronoi neighbors. That tends to be a safe choice in terms of computation time required, even with a Radius set to 0, meaning automatic Radius.

Increasing the number of neighbors when using automatic Radius can dramatically increase the time required for the computation. With zero neighbors, the computation takes about 29 seconds on a Ryzen 9 3900x (24 hypercore) CPU. Changing that to 2 neighbors without manually specifying the Radius to some number other than 0 (that is, not automatic) can take many minutes. See the discussion in the Notes section at the end of this topic.

The Result destination for Interpolation operations is always New Table, meaning to create a new table and a new image to visualize the tiles in the new table. We enter bbmesa Kriged terrain as the name of the New image. As we enter the name for the drawing, the pane will automatically fill in an analogous name for the table. We can change that if we like. We can use whatever names we want, but it makes sense to use names that will later remind us what those components are supposed to be.

If we would like to see a preview, we press the Preview button.

The preview shows the proposed result in shades of blue preview color, with a blue preview caption bar displayed at the top of the window giving the name of the template used for the preview.

Previews on big data or complex tasks can take a while on smaller machines, just as long as calculating the final result, so we may want to skip doing a preview and go straight to pressing Transform to compute the final result.

Press Transform.

Manifold thinks for a while, 20 to 30 seconds on a faster machine and perhaps a minute or more on a slow machine, applying all the CPU cores we have to the Kriging task, to do in a few tens of seconds what might require many minutes or even hours in a non-parallel system. The result is that a new bbmesa Kriged terrain image and its table appears in the project.

We drag and drop the new bbmesa Kriged terrain image into the map.

We drag and drop the new image into our map. We also turn off the drawing layer and the original PDS terrain image layer so those will not interfere with the visualization we create.

The new image appears in default, grayscale formatting, with no hill shading, since it has not yet been styled. That is easy to do.

With the focus on our new image layer in the map, we choose the Style pane. It automatically configures itself with controls relevant to images. We begin by changing the (use RGBA channels) setting to Channel 0, the one channel in the image.

We next choose 10 breaks, Tally, and then we add a palette. We press the palette button and choose the Spectral palette from the Color Brewer collection. See the Style: Images topic for details.

Press Update Style.

That results in a more interesting image, but we need to hill shade to bring out terrain relief.

Back in the Style pane, we press the Options tab and check the Use shading box, leaving default settings as is. We press Update Style.

Much better.

We pan the map to get a better look at the big craters on Buckboard Mesa. To see if our new terrain elevation image is correctly georeferenced, we will use the Layers pane to make it partially opaque, so the Bing layer can show through from underneath.

In the Layers pane we double-click into the opacity setting for our terrain image layer and change opacity to 50%.

Right away, the layer changes to partial opacity, so the Bing layer shows through. We can see the georegistration is perfect.

If we double-click our interpolated layer's tab to turn off the layer, we can see how overlaying a partially transparent terrain relief image layer brings out terrain relief that is difficult to see in the Bing satellite photo.

Notes

GPU and Kriging - Surprisingly, most of the Kriging task involves very little mathematics and will be preferentially dispatched to parallel, manycore CPU over GPU. Parallel GPU is used at the very beginning of the job to compute the mathematically intensive Kriging model and then for rest of the calculation parallel CPU tends to do the job, especially if we have a manycore processor like a 24 hypercore Ryzen 9 (12 physical cores) or other modern manycore processor.

Manually setting Radius - Using a Radius of 0 tells Manifold to autocompute the Radius to be used. Manifold uses algorithms to match ESRI computations for Radius, to provide results that will be no surprise to government and other users where ESRI is popular. However, using the automatic value can result in very large Radius, and thus slower performance when manually specifying the number of Neighbors.

Consider the illustration above, where we have measured outward from the mass of points using the Measurements tools to measure a distance of about 15 meters. Within the mass of points there are very many points within 15 meters of each other, but beyond the edge of the mass of points there are none. We gain nothing by using a larger Radius than that, but the ESRI style algorithms can generate a larger Radius.

Therefore, it makes sense to use the Path tool to note what a reasonable distance is and to then manual specify that distance as the Radius when manually specifying the number of Neighbors. For example, if we used 10 as the number of Neighbors we could use 15 as the Radius. Using those values, the computation takes about 33 seconds on a Ryzen 9 3900x (24 hypercore) CPU.

Videos

Manifold Viewer - Speed Demo with 1280 GPU Cores - 2 Minutes vs 5 Days - Watch the free Manifold Viewer run CPU parallel and GPU parallel with 8 CPU cores and 1280 GPU cores to absolutely crush a job, doing in 2 minutes what would take non-GPU-parallel software 5 days. The video shows Viewer instantly pop open a 4 GB project that contains a huge, multi-gigabyte terrain elevation surface for Crater Lake, Oregon. With a point and click - no parallel code required - we compute the mean curvature at each pixel of the surface using a 7x7 matrix in under two minutes. We combine that with the original surface for enhanced hill shaded effects to better see details. Using an 11x11 matrix takes just over two minutes, a huge computation that takes days in non-GPU-parallel GIS packages.

Create Terrain Elevation from a NASA PDS Table - Do in six minutes what can take hours in other GIS packages! We import terrain elevation raster data from a NASA PDS archive and style the images with an elevation palette and hill shading. One of the images was wrongly georegistered by NASA, so we import the original LiDAR point data from a table in the PDS archive. We transform the imported table from 0 to 360 degree longitudes into +/-180 longitudes, construct a geometry field, and then we create a drawing from the table that shows the LiDAR point data in the table. Next, we use Kriging to create a terrain elevation raster from the LiDAR point drawing, and then we style that and hill shade it so it can be exactly overlaid on a Bing satellite layer to enhance the satellite photography with terrain relief. See how Manifold's wonderfully efficient user interface and blisteringly fast parallel performance makes GIS work fast and easy!