Example: Transfer DEM Terrain Heights to Areas in a Drawing

Given a map with an image layer that shows terrain heights taken from a DEM, and a drawing layer that contains areas, using a small SQL query we transfer the average terrain height within each area to that area as a Height attribute for the area.   If our drawing had points instead of areas, this example would work exactly as written, transferring the height of the surface under each point into that point object's Height field.

 

An easier way to accomplish this example is to use the point and click Edit - Join dialog, which requires no SQL.   See the Edit - Join topic for step by step examples.  

 

To use the technique illustrated, the terrain elevation image from which heights are computed must have only a single channel, which is typical for raster images that convey data such as terrain elevations or other data.

 

The TileGeomAvg SQL function used in this example requires that both arguments, the raster image as well as the vector drawing, use the same coordinate system.  We can accomplish that either by using the same projection for both (it takes but a moment to reproject a component to ensure both use the same coordinate system), or by using the technique shown in this example, where we first create a coordinate converter object and then use CoordConvert to on the fly to convert the drawing's coordinate system into the same coordinate system as used by the image.  

Our Data

We start with the map seen below.  The 1475 PM layer is a raster image surface showing terrain elevation in the vicinity of Livermore, California.  It is the same the data set imported and styled in the Example: Import DDF SDTS DEM Raster File  topic and also used in the Style: Contouring using Colors topic.  

 

 

The Drawing has four areas in it.  The Drawing layer is shown with 50% opacity using the Layers pane so the terrain below the areas can partially show through.   The map also has a Height Labels layer, a labels layer which currently does not show anything since there are, as yet, no height values in the areas.

 

The terrain elevation image has been formatted using the following color and height intervals:

 

 

The drawing contains four areas, each with a Height field that is an INT64 type, an integer.

 

 

The Height for each area is currently empty.  

Our Task

Our task for each area is to find the height for each pixel of the terrain image that falls within the area, compute the average height for all those pixels, and to put the average height into the Height field for that area.

 

We can do that easily using the TileGeomAvg SQL function.

 

 

We launch View - New Command Window - SQL and enter the above query, copying and pasting it from the text below, if we like.

 

-- $manifold$

 

-- prepare to convert coordinates from drawing to image

 

VALUE @conv TABLE = CALL CoordConverterMake(

  ComponentCoordSystem([1475 PM]),

  ComponentCoordSystem([Drawing])

);

 

-- set 'height' field to average height taken from image

 

UPDATE (

  SELECT [mfd_id], [Height],

    TileGeomAvg([1475 PM], 0, CoordConvert(@conv, [Geom]))

      AS [computed]

  FROM [Drawing]

  THREADS SystemCpuCount()

) SET [Height] = [computed];

 

 

The VALUE statement creates a coordinate converter object that converts the coordinate system of the drawing into the coordinate system of the image, and then declares that object as a global value that can be used in the rest of the query.

 

The UPDATE statement uses the TileGeomAvg function to compute the average value of pixels under each Geom in the drawing, with the coordinate system used by the Geom converted on the fly into the coordinate system used by the image with CoordConvert and the coordinate converter global value @conv.

 

We press the ! Run button in the main toolbar to run the query.  Instantly, it computes the average height under each area and places the result into the Height field for each area:

 

 

The computation results in a floating point number.  Since we want simple labels, we have created the Height field in the table as an integer data type.  When floating point values are stored into that attribute they are automatically rounded into an integer value.

 

Now that there are Height values in the drawing's table, the Height Labels layer swings into action and shows those heights in the map.

 

 

Comparing the numbers we see for average height in each area to the height intervals used for coloring the terrain elevation, we can see the heights have been calculated accurately.

 

 

The blue regions in the surface start at 144 feet and go up to 275 feet, so an average Height of 193 for the leftmost area makes sense. The green regions go from 275 to 405, so an average height of 294 for the area in the center also makes sense, since it is right on the edge of the green region and partially overlaps canyons in blue.   The pale yellow regions of the surface go from a height of 405 to 536, so the average heights for the other two areas also seem reasonable.

 

Notes

When projections are the same -  If we know for sure that the coordinate systems are exactly the same for Drawing and for the 1475 PM image, we could have used a simpler query, which would not do coordinate system conversion on the fly:

 

UPDATE (

  SELECT [mfd_id], [Height],

    TileGeomAvg([1475 PM], 0, [geom]) AS [computed]

  FROM [Drawing]

  THREADS SystemCpuCount()

) SET [Height] = [computed];

 

This is tricky business, however, since every detail of the coordinate systems, including the metrics, must be identical for the simpler query to work.   It is safer to use the approach shown in the example query in this topic, first creating a coordinate converter object, and then using that coordinate converter object to convert coordinates on the fly.  That always works, regardless of the coordinate systems used for the image and the drawing.

 

Adapting the query - We can use the above query in our own work by simply substituting the name of our image for 1475 PM, the name of our drawing for Drawing, and the name of the height field in our drawing for Height.  Ctrl-h for search and replace, a standard Windows shortcut, works within the Command Window.

 

If we prefer to compute a different statistic instead of the average height under each area, we can replace the TileGeomAvg function with a different TileGeom series SQL function, such as TileGeomMax for the maximum height under each area, TileGeomMin for the minimum height under each area, TileGeomSum for the sum of all pixel values within the area, or other functions.

 

A short guide to statistics from other TileGeom functions that can be substituted in the above query:

 

TileGeomAvg

The average value found in pixels that are located within the area.

TileGeomCount

The number of pixels that are located within the area.

TileGeomDiversity

The total number of different values found in pixels that are located within the area.

TileGeomDiversityIndex

The diversity index for the pixels that are located within the area.

 

A diversity index provides a measure of diversity, computed as 1 - sum(individualcount^2) / (totalcount^2).  A diversity index of 0 means that all values are equal.

TileGeomMajor

The most frequently occurring value found in pixels that are located within the area.

TileGeomMax

The maximum value found in pixels that are located within the area.

TileGeomMedian

The median value found in pixels that are located within the area.

TileGeomMin

The minimum value found in pixels that are located within the area.

TileGeomSum

The sum of values found in pixels that are located within the area.

 

When might we use TileGeomSum?  Suppose the value for each pixel in a raster image is not height, but indicates the population at that pixel.  If we wanted to know the total population enclosed within an area in a drawing we would use TileGeomSum to get the sum.

See Also

Maps

 

Drawings

 

Tables

 

Images

 

Labels

 

Data Types

 

Queries

 

Style: Drawings

 

Style: Images

 

Style: Labels

 

Style: Images

 

Style: Contouring using Colors

 

Layers Pane

 

Command Window

 

Example: Import DDF SDTS DLG Vector File 

 

SQL Example: Acquire Pixel Values within Areas - Given a map with two layers, an image layer that provides terrain elevations and a drawing with area objects, we use the TileGeomToValues function to find the three highest pixels in the image within each area object in the drawing.

 

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.