Example: Enhance Terrain with Curvatures

In this example, we enhance a terrain showing Crater Lake, Oregon, by using mean curvature calculation to bring out details.   The example uses a 4 GB project containing a large terrain elevation surface.  Using a point-and-click dialog with no SQL, we apply automatic CPU parallelism and GPU parallelism to absolutely crush a task in two and a half minutes that would take non-parallel software days.   This topic is a documentation version, with a few additional notes of the Speed Demo with 1280 GPU Cores video.

Example Data

The data we use was downloaded from the USGS website at https://pubs.usgs.gov/ds/716/, a page titled High-Resolution Digital Elevation Dataset for Crater Lake National Park and Vicinity, Oregon, based on LiDAR Survey of August-September 2010 and Bathymetric Survey of July 2000.     The data was downloaded in a zip file containing an ESRI ArcGrid of elevations in ADF format.   The zip file is  2.26 GB compressed and 4.18 GB expanded.

 

 

Importing into Manifold, we get an image that is 28359 pixels wide and 36416 pixels tall, using float32 data type for the elevation of each pixel.   We have used the Style pane to color the image by elevation using a thematic format, and we have also applied hill shading.

 

 

Zooming in to take a closer view, we see the data set is a composite of higher resolution LiDAR data at one meter resolution for terrains on land, and lower resolution sonar data for underwater portions of Crater Lake, which have been re-sampled to the same 1 meter pixel size as the land data and then combined with the LiDAR data to form a single data set.

 

 

Zooming further in (and resizing the Map window slightly into horizontal format) we see there are portions of the data set where there is not a clean join between the underwater, sonar data and the land surface LiDAR data.   Invisible pixels appear where there are gaps, allowing the white background layer of the map to show through.     That is mildly distracting, so we will fix that by assigning a green color to the background of the map.

 

 

In the Layers pane, we double-click into the white color well of the Background layer.  

 

 

In the color dialog that pops open we click the Color PIcker choice.

 

 

The cursor switches to a color picker eyedropper and a slight pale overlay appears everywhere we can click to pick the color at that spot.  That covers all the space on all monitors we have, so we can click anywhere we like to grab the color at that location.   We click onto a green portion of the display near a region of missing data.   If we set the background to that color, it will seem to fill in the region of missing data.

 

 

The color at the location we clicked appears in the Layers pane.

 

 

Immediately, in the map window the missing pixels appear to be filled in with adjacent color.  They are not actually filled in, since the pixels where data is missing are still transparent.  However, now those transparent pixels allow the color of the background layer to show through.   We picked the color for the background layer from immediately adjacent terrain so the color of the background layer is such a good visual match to what the missing data pixels should be that they appear to be filled in.

Calculating Curvature

Curvature is a measure of how curved the surface is at a particular location, with different Curvature templates made available within the Transform pane to allow computing curvature in different ways.   The Curvature, Mean template calculates an overall measure of curvature as discussed in the Speed Demo with 1280 GPU Cores  video.  Other forms of curvature we might use, such as Curvature, Plan or Curvature, Profile, are discussed in a fine ESRI blog entry and links therein on how curvature calculations can be used to enhance the presentation of terrain elevation data.

 

With the focus on the open map window, in the Transform pane we choose the craterlake_1m layer and then we chose the Tile field.  We double-click the Slope template to launch that template.

 

 

In the Slope template we choose channel 0 for the Channel, and we choose mean curvature for the Operation.

 

The system we use has a GPU in it (of course), so we have no hesitation about specifying 3 for the Radius, thus specifying a 7x7 convolution matrix, as discussed in the How Matrix Filters Work topic.   We press Add Component.

 

With a GPU and Manifold CPU parallelism + GPU parallelism, the calculation on a terrain surface of this time takes a minute or two, much faster than the hours or days it would take with non-parallel software.  That is why non-parallel software is usually limited a much more easily computed small matrix size, 3x3, with a Radius of 1.  

 

Using a Ryzen 9 3900x (24 hypercores) with a GTX 1060 6GB card (1280 CUDA cores), the computation takes a minute and a half.  Using a very old Core i7 960 (8 hypercores) with a really old GT 710 (192 CUDA cores), the computation takes four minutes.  Even with a really old, virtually free, GPU the computation goes much faster than not having a GPU, although having a manycore CPU even without a GPU will go very much faster than can be accomplished with non-parallel software that can run only a single CPU core.

 

See the GPGPU topic for Windows Task Manager illustrations showing CPU and GPU activity while computing mean curvature in this example topic.

 

For the Result destination we choose New Table and then we enter Mean Curvature as the name of the New Image and Mean Curvature Tiles for the name of the corresponding New Table.  We can specify whatever names we want, but it is wise to use short names that will remind us what the new image and table are about.

 

The default table name created automatically as we enter a name for the New image ends in ...Table.  This topic uses a table name that ends in ...Tiles, which is how some users prefer to name tables for images.  Others prefer to have all of their tables, both for drawings and images, end in ...Table.  We can adjust the name automatically created for the table however we like.

 

If our data is not too large and our hardware is fast, we might want to see a preview.  It takes as long to compute a preview as doing the actual transform, so when working with larger data or slower machines we may as well just apply the transform to see what happens.

 

However, if we like we can press Preview to see a preview.

 

 

The preview will appear in shades of blue preview color.

 

 

The preview layer is opaque by default, hiding layers below.  To see the layers below, we can right-click onto the preview caption bar and choose partial opacity, for 75%, 50%, or 25% opacity.

 

 

Setting 50% opacity allows the layer below to be seen, visually combining with the preview layer.

 

 

Another way to see layers below the preview layer is to right-click the blue caption bar and then choose either Left or Right to show the preview only for the left or the right half of the screen.  We can then drag the blue vertical dividing line between preview and no preview to the left or right as desired.   This is a convenient way of seeing both the preview and parts of the layers underneath the preview.   As seen in the illustration above we can combine split screen with partial opacity to use both effects at the same time.

 

If we like what we see in the preview, or if we simply do not have the patience to do a preview first, we can apply the transform by pressing Transform.

 

 

The results appear as a new image and table appear in the Project pane, named as we specified.

Visually Combining Curvature with Terrain

We drag and drop the new image into the map.

 

 

It appears georegistered within the map (of course), using default styling, where Manifold has set the range of grayscale values based on a quick estimation of what the maximum and minimum floating point pixel values are in the data.    That is a reasonable visualization, but for our purposes in blending this layer with our main layer, we would prefer less contrast.  We can do that using the Style pane.

 

With the focus on the Mean Curvature layer in the map window, in the Style pane we click into the channel rows to move the focus there.

 

 

We press Ctrl-A to select all of the channels, and then we right-click into the 0-255 range portion of any of the three R, G, or B channels to launch the context menu.  We choose Full Range.   The system takes a moment to compute new values.  

 

When the new values appear, press Update Style.

 

 

Next, we again right-click into the 0-255 range portion of any of the three R, G, or B channels to launch the context menu.  We choose Autocontrast and then the 2.0 std autocontrast value, as discussed in the Style: Autocontrast topic.

 

Press Update Style.

 

 

Immediately, the mean curvature image is styled to appear in medium contrast grayscale, the gray coloring showing regions of greater curvature as lighter color.   We will now visually combine that curvature image with the original image.

 

 

We drag the craterlake_1m tab to the left to move it up in the layer stack.  It now completely covers the curvature image.

 

 

In the Layers pane, we double-click into the opacity number of the craterlake_1 layer and change it to 80, for 80%.  

 

 

80% opacity allows some of the curvature image layer to show through, which visually enhances the terrain layer as seen above, and as seen below in a more zoomed-in view.

 

 

Comparing the image above with the un-enhanced terrain display below (shown with 100% opacity, so none of the curvature layer below blends with the upper layer), we see that adding curvature brings out more details, especially terracing in the lava flow at upper left in the map window.

 

 

In this case, we have used only one layer of curvature.  We could have computed profile curvature and used that with partial transparency as well, to combine visual effects from different curvatures, each adding a slight enhancement to the original surface layer.

Videos

Manifold Viewer - How Matrix Filters Work - The easy, simple way to learn how filters work! Watch this action-packed video using Manifold Viewer that illustrates how matrix filters, also known as convolution filters, work. In addition to explaining filters, the video provides a real-life look at simple Manifold techniques for moving objects around in drawings using the Shift transform, and fast and easy use of Selection and tables to quickly put desired values into attributes. Sound technical? Sure, but in a very easy and simple way.

 

Manifold Viewer - Create Custom GPU Accelerated Filters in Seconds - A technical video using the free Viewer showing how to create your own, fully custom, fully GPU-parallel, convolution matrix filters, including Emboss, Sobel, Prewitt, and Kirsch edge detection and many more, for use in Viewer or Release 9. Modify the spatial SQL examples in the downloadable example project to specify a custom matrix and in seconds your custom filter can do image processing at GPU-parallel speeds. Viewer is read-only, but you can copy and paste the query text for custom filters to and from Notepad or any other text editor. Download the Custom_Filter_Examples.mxb sample project from the Examples page on the Manifold website to try out the video in Viewer or Release 9.

 

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 hillshading. 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 hillshade 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!

Other Links

The Quick and Dirty Introduction to the Curvature of Surfaces - Mean curvature explained in pictures and a few very easy equations.

 

See the ESRI blog entry on curvatures and the ESRI blog entry on enhancing symbology for a great discussion of how curvature calculations can be used to enhance the presentation of terrain elevation data.

 

For very large, very detailed images of the Crater Lake data set enhanced with mean curvature and profile curvature, See the Eleven Views of Crater Lake images in the Gallery page.

 

The data set used:  Robinson, J.E., 2012, High-resolution digital elevation dataset for Crater Lake National Park and vicinity, Oregon, based on LiDAR survey of August-September 2010 and bathymetric survey of July 2000: U. S. Geological Survey Data Series 716. (Available at https://pubs.usgs.gov/ds/716/.)

See Also

GPGPU

 

Maps

 

Images

 

Data Types

 

Layers Pane

 

Transform Pane

 

Style: Images

 

Style: Autocontrast

 

Transform - Tiles

 

How Matrix Filters Work

 

Command Window

 

SQL Functions

 

ADF, ESRI ArcGrid

 

Example: Display an NAIP Four Band Image as Color Infrared (CIR) - How to use the Style pane for images to re-assign channels in a four band NAIP image to produce a Color Infrared (CIR) image display.

 

SQL Example: Process Images with 3x3 Filters -  Shows a step-by-step example of developing an SQL query that takes a query written by the Edit Query button and then modifies that query into a general purpose query that can apply any 3x3 filter.   This makes it easy to use matrix filters we find on the web for custom image processing.   We extend the query by using parameters and adding a function, and then show how it can be adapted to use a 5x5 filter.

 

SQL Example: Process Images using Dual 3x3 Filters  - A continuation of the above topic, extending the example query to utilize two filters for processing, as commonly done with Sobel and Prewitt two filter processing.

 

SQL Example: Process RGB Images using Matrix Filters - A continuation of the above two topics, extending the example query to process three channel, RGB images.

 

SQL Example: Create NDVI Displays - How to create a query that creates an NDVI display from a four-band NAIP image, with tips and tricks on how to copy and paste existing information to get the result we want.

 

SQL Example: Create Topographic Position Index TPI Displays - In this example, we use a few short lines of SQL to create a Topographic Position Index (TPI) display.  TPI characterizes the undulations of a terrain elevation surface.  TPI value above zero show locations that are higher than then average of immediately surrounding terrain, and thus tend to show ridges.   TPI values below zero show locations that are lower than the average of immediately surrounding terrain, and thus tend to show valleys.  TPI values that are zero show areas of constant slope.

 

Example: Create Terrain Elevation Raster from a NASA PDS Table - 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 hillshade it so it can be exactly overlaid on a Bing satellite layer to enhance the satellite photography with enhanced terrain relief.

 

Example: Rearrange Channels using an Expression - We use a simple expression in the Transform pane to rearrange the order of channels within the data.