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.

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.

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 first part is a series of comment lines that report the settings which were used in the Transform pane when the Edit Query button was pressed. In future illustrations in this topic we will drop those, to reduce the size of illustrations and the length of the topic.
The second part is between the -- prepare begin and -- prepare end comment lines, and creates the destination tile table and image.
The third part uses the TileFilterDefSquare function with a radius of 3 and a center weight of 1, which results in an average filter, to do the average computation. In the listing above we have added some extra whitespace for greater legibility.
Finally, the query updates pyramids and the Rect property for the image that is created.

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.

SQL Example: Create Topographic Position Index TPI Displays

Videos

Other Links

See Also