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.

 

Calculating TPI is simple in principle:  for each pixel of the terrain, we take the elevation value of that pixel and we subtract the average elevation of immediately surrounding pixels.  Manifold provides matrix filter calculation functions to compute the average value over a given matrix, and they are even GPU accelerated.  

 

We use the same data set of terrain and bathymetry in the region of Crater Lake, Oregon, that is used in the Example: Enhance Terrain with Curvatures topic.

 

 

In the illustration above, we have zoomed into the Western rim of Crater lake, with the volcanic cone of Wizard Island visible to the right.   The flattened terrace near the bottom of the Western rim's slope is approximately the waterline of Crater Lake.

 

The easiest way to write SQL is to have Manifold write it for us.   We can use the Edit Query button in the Transform pane to have Manifold write a query for us that computes the average for each pixel and saves the result into a new image.  We will then modify that query to calculate TPI.   

 

This approach has the virtue that Manifold will write all SQL necessary to create a new image and to populate it with the results of our calculation, while taking care of all infrastructure issues such as declaring the right coordinate system for that new image as well as all other necessary salad dressing.

 

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

 

 

In the Filter template, we choose channel 0 for the Channel and average as the Filter.

 

We leave the Shape parameter as the default square.

 

We enter 3 as the Radius.

 

For the Result destination, we choose New Table and then enter TPI as the name of the New Image and TPI Tiles as the name of the New table for that image.  We are going to edit the resulting query to create TPI data, so we may as well have Manifold use indicative names in the components that the query text uses.

 

Press Edit Query.   

 

That launches a Command Window that is populated with an SQL query which does what the Filter : average transform operation does given settings of 3 for the Radius and square for the Shape, saving the results to newly created components called TPI for the image and TPI Tiles for the table.

 

 

The full text of the above query is:

 

-- $manifold$

--

-- Auto-generated

--

-- Filter

--   Layer: craterlake_1m

--   Field: Tile

--   Channel: channel 0

--   Filter: average

--   Shape: square

--   Radius: 3

--   Result: (new table)

--   Channel type: float64

--   New image: TPI

--   New table: TPI Tiles

--   Resources: all CPU cores, all GPU cores

--   Transform selection only: FALSE

--

 

-- prepare begin

 

CREATE TABLE [TPI Tiles] (

  [X] INT32,

  [Y] INT32,

  [mfd_id] INT64,

  [Tile] TILE,

  INDEX [mfd_id_x] BTREE ([mfd_id]),

  INDEX [X_Y_Tile_x] RTREE ([X], [Y], [Tile] TILESIZE (256, 16) TILETYPE FLOAT64),

  PROPERTY 'FieldCoordSystem.Tile' ComponentFieldCoordSystem([craterlake_1m Tiles], 'Tile'),

  PROPERTY 'FieldTileSize.Tile' '[ 256, 16 ]',

  PROPERTY 'FieldTileType.Tile' 'float64'

);

CREATE IMAGE [TPI] (

  PROPERTY 'Table' '[TPI Tiles]',

  PROPERTY 'FieldX' 'X',

  PROPERTY 'FieldY' 'Y',

  PROPERTY 'FieldTile' 'Tile',

  PROPERTY 'Rect' CAST(ComponentFieldBounds([craterlake_1m Tiles], 'Tile', FALSE) AS NVARCHAR)

);

 

-- prepare end

 

VALUE @image TABLE = CALL ComponentFieldImage([craterlake_1m Tiles], 'Tile');

VALUE @filterDef TILE = TileFilterDefSquare(3, 1);

DELETE FROM [TPI Tiles];

INSERT INTO [TPI Tiles] (

  [X], [Y],

  [Tile]

) SELECT

  [X], [Y],

  CASTV(TileRemoveBorder(TileFilter(

    TileCutBorder(@image, VectorMakeX2([X], [Y]), 3), 3, @filterDef), 3) AS FLOAT64)

FROM [craterlake_1m Tiles] THREADS SystemCpuCount();

 

TABLE CALL TileUpdateFieldPyramids([TPI Tiles], 'Tile');

ALTER IMAGE [TPI] (

  ADD PROPERTY 'Rect' Coalesce(

    CAST(ComponentFieldBounds([TPI Tiles], 'Tile', TRUE) AS NVARCHAR),

    ComponentProperty([TPI], 'Rect'))

);

 

The query is in four major parts.   

 

 

The third part is the key portion of the query, with the rest being supporting infrastructure.  We will change that key portion of the query to not just compute an average, but instead to first take the value of each pixel and to subtract the average from that.  That is extremely easy to do:

 

We change this:

 

  CASTV(TileRemoveBorder(TileFilter(

    TileCutBorder(@image, VectorMakeX2([X], [Y]), 3), 3, @filterDef), 3) AS FLOAT64)

 

Into this:

 

  [TILE] - CASTV(TileRemoveBorder(TileFilter(

    TileCutBorder(@image, VectorMakeX2([X], [Y]), 3), 3, @filterDef), 3) AS FLOAT64)

 

In the above, we have added  "[Tile] - " before the CASTV expression.

 

In addition, we will use search and replace to change the name of the created components, so instead of using Tiles Average they will use TPI.   The modified query appears below:

 

 

The full text of the modified query is:

 

-- $manifold$

 

-- prepare begin

 

CREATE TABLE [TPI Tiles] (

  [X] INT32,

  [Y] INT32,

  [mfd_id] INT64,

  [Tile] TILE,

  INDEX [mfd_id_x] BTREE ([mfd_id]),

  INDEX [X_Y_Tile_x] RTREE ([X], [Y], [Tile] TILESIZE (256, 16) TILETYPE FLOAT64),

  PROPERTY 'FieldCoordSystem.Tile' ComponentFieldCoordSystem([craterlake_1m Tiles], 'Tile'),

  PROPERTY 'FieldTileSize.Tile' '[ 256, 16 ]',

  PROPERTY 'FieldTileType.Tile' 'float64'

);

CREATE IMAGE [TPI] (

  PROPERTY 'Table' '[TPI Tiles]',

  PROPERTY 'FieldX' 'X',

  PROPERTY 'FieldY' 'Y',

  PROPERTY 'FieldTile' 'Tile',

  PROPERTY 'Rect' CAST(ComponentFieldBounds([craterlake_1m Tiles], 'Tile', FALSE) AS NVARCHAR)

);

 

-- prepare end

 

VALUE @image TABLE = CALL ComponentFieldImage([craterlake_1m Tiles], 'Tile');

VALUE @filterDef TILE = TileFilterDefSquare(3, 1);

DELETE FROM [TPI Tiles];

INSERT INTO [TPI Tiles] (

  [X], [Y],

  [Tile]

) SELECT

  [X], [Y],

  [TILE] - CASTV(TileRemoveBorder(TileFilter(

    TileCutBorder(@image, VectorMakeX2([X], [Y]), 3), 3, @filterDef), 3) AS FLOAT64)

FROM [craterlake_1m Tiles] THREADS SystemCpuCount();

 

TABLE CALL TileUpdateFieldPyramids([TPI Tiles], 'Tile');

ALTER IMAGE [TPI] (

  ADD PROPERTY 'Rect' Coalesce(

    CAST(ComponentFieldBounds([TPI Tiles], 'Tile', TRUE) AS NVARCHAR),

    ComponentProperty([TPI], 'Rect'))

);

 

To run it, we press the ! Run button on the main toolbar.  A new image called TPI and the image's table appear in the Project pane.

 

 

We drag and drop the image into our map, and Style it using the settings below in the Style pane.  Brighter spots show ridges, and darker spots show valleys.

 

 

In the Style pane we selected the three range rows for the channels, right-clicked, and chose 2.0 std autocontrast for a medium level of automatic contast.

 

 

We can combine the TPI image with the original image as seen above, positioning the original image layer as the top layer and setting it to 80% transparency in the Layers pane, thus allowing some of the brighter and darker pixels from the TPI layer below to modify the visual appearance of the combined map.  We can compare the effect to the original, shown below.

 

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.

Other Links

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

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.

 

Example: Enhance Terrain with Curvatures -  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.

 

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