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.
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 use the TileGeomToValues function to find the three highest pixels in the image within each area object in the drawing. That is very quick. Most of this example topic shows how to use the result in a way that teaches useful methods.
We first illustrate the query in action to create a results table in the Command Window. We next alter the query slightly to write the result into a new table, and then we create a drawing from that table. We add the drawing to the map so we can see where the three highest pixels are within each area.
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 TileGeomToValues 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 used by the image.
Important: Results values created by the TileGeomToValues function include X and Y coordinates in whatever coordinate system is used by the image. For everything to line up, if the image's coordinate system is applied to a new drawing created from those results, the applied coordinate system must override metrics. That is easy to do using simple steps shown in the Favorite Coordinate Systems topic. In this example, we show additional steps using the Favorites dialog to create a favorite that overrides metrics.
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 terrain elevation image has been formatted using the following color and height intervals:
Popping open the drawing's table, we can see it contains only areas.
Our task is to find the three highest pixels in the 1475 PM image that fall within each area. We can do that using the TileGeomToValues 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$
-- take all pixels, return n highest
FUNCTION highestCoords(@t TABLE, @n INT32) TABLE AS (
SELECT * FROM @t ORDER BY [value] DESC FETCH @n
) END;
-- prepare to convert coordinates from drawing to image
VALUE @conv TABLE = CALL CoordConverterMake(
ComponentCoordSystem([1475 PM]), -- to
ComponentCoordSystem([Drawing]) -- from
);
-- select 3 highest pixels under each object in the drawing.
SELECT [mfd_id], SPLIT CALL highestCoords(
CALL TileGeomToValues([1475 PM], CoordConvert(@conv, [Geom])), 3
) FROM [Drawing];
The FUNCTION statement defines a highestCoords function for the use of this query. The function takes a table of pixels t and an integer number n, and it returns a table giving the n highest pixel values in the table.
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 SELECT statement first uses the TileGeomToValues function to return a table of all pixels values in the image that fall within each Geom in the drawing. The coordinate system used by the Geom is converted on the fly into the coordinate system used by the image with CoordConvert and the coordinate converter global value @conv. Next, the table returned by the TileGeomToValues function is processed by the highestCoords function with an integer argument of 3 to return a table of the three highest coordinates. The SPLIT clause works with the SELECT to generate permutations so we end up with a table of all geoms with the three highest pixel values for each.
We press the ! Run button in the main toolbar to run the query.
Instantly, the query generates a results table that lists the X, Y coordinates and pixel value for the three highest pixels within each area in the drawing. The image uses Pseudo-Mercator projection, so the X, Y coordinates are in meters, as adjusted by the metrics in use by the image.
The above query and results table shows
how the function works, the primary goal of this topic. We
can now continue the topic to show how the results of the function might
be used in real life, to illustrate how various other dialogs and features
are used.
While the above result is wonderful, in this example we would prefer to capture the results within a new table. We can then create a drawing based on that table, to show the points where the highest pixels are located.
We modify our query slightly, to convert the SELECT into a SELECT ... INTO statement, creating a new table called Points Table.
-- $manifold$
-- take all pixels, return n highest
FUNCTION highestCoords(@t TABLE, @n INT32) TABLE AS (
SELECT * FROM @t ORDER BY [value] DESC FETCH @n
) END;
-- prepare to convert coordinates from drawing to image
VALUE @conv TABLE = CALL CoordConverterMake(
ComponentCoordSystem([1475 PM]), -- to
ComponentCoordSystem([Drawing]) -- from
);
-- select 3 highest pixels under each object in the drawing.
SELECT SPLIT CALL highestCoords(
CALL TileGeomToValues([1475 PM], CoordConvert(@conv, [Geom])), 3
) INTO [Points Table] FROM [Drawing];
In addition, we drop the mfd_id field from the SELECT list, since we cannot create multiple records with the same mfd_id field in a table.
We press the ! Run button in the main toolbar to run the query.
The query creates a new Points Table that contains the results of the query. The table has no mfd_id identity field and no index, so it is read-only and displays fields in gray background color used for read-only fields.
The Points Table created above is a geocoded table, in which each record has an X, Y location, but not yet a Geom field that captures the geometry as a point. We will add a Geom geometry field to the table using the process shown in the Example: Create a Drawing from a Geocoded Table topic. But first, we make the table writable by adding an identity field and index, using the procedure shown in the Add an Index to a Table topic.
With the focus on the opened table, we choose Edit - Schema from the main menu, or press Ctrl-E as the keyboard shortcut, to launch the Schema dialog.
We
click the Add Identity button
in the toolbar.
A new mfd_id field of type int64 and a new mfd_id_x index are provisionally added to the schema. The mfd_id_x index is a btree index on the new mfd_id field. The new field and index are shown with provisional, bluish background color to indicate the addition of the field and index have not yet been committed to the schema.
We press Save Changes to update the schema, to add the new field and index to the table.
An mfd_id field appears in the table, and the table switches to white background color to show it is writable. In the illustration above we have used the Layers pane to move the new mfd_id field to the left of the table.
We choose Edit - Schema again, or press Ctrl-E, to launch the Schema dialog so we can create a Geom field.
Press the Add command
button, and choose Field in
the drop down menu.
We enter the name Geom and choose type geom. We press OK.
There is a shortcut
we could have done in the above dialog, to make a slight adjustment in
the coordinate system. In this example, instead of taking the shortcut
in the above dialog we will make the adjustment later on, to allow a more
explicit discussion of what the adjustment accomplishes.
Back in the Schema dialog we see the new field in provisional, bluish color. To commit changes to the table, adding the field, we press Save Changes.
A new Geom field appears in the table, as yet unpopulated.
Using the process shown in the Example: Create a Drawing from a Geocoded Table topic we use the Transform pane to populate the Geom field.
We choose the Geom field as the Target and then we click on Compose Point to choose that template. In the X and Y boxes that appear we choose the X and Y fields respectively.
Press Update Field to commit the action.
The Compose Point template populates the Geom field with geoms that contain point objects, each of which is located at the X, Y coordinate specified by the X and Y fields.
We can now right-click onto the Points Table in the Project pane and choose Create - New Drawing to create a drawing.
We enter a name of Points for the new drawing, leave defaults as is, and then press Create Drawing.
There is a shortcut
we could have done in the above dialog, to make a slight adjustment in
the coordinate system. In this example, instead of taking the shortcut
in the above dialog we will make the adjustment later on, to allow a more
explicit discussion of what the adjustment accomplishes.
A new drawing, called Points, appears in the Project pane.
We drag and drop the Points drawing into the map, but no points appear. The reason no points appear is that the projection used by the new drawing was created using Manifold's default Pseudo-Mercator coordinate system, which does not override metrics. However, the X, Y coordinates used to populate the geoms for the points are values in the coordinate system of the 1475 PM image. The coordinate system of the 1475 PM image uses metrics that are different from defaults, the usual case with images. In coordinate systems for images, metrics are a key part of how pixel sizes are defined.
We need to adjust the coordinate system used by the Points drawing to use identically the same metrics as used by the 1475 PM image. We could do that by launching the Coordinate System Metrics dialog for the 1475 PM image, writing down the metrics, and then launching the Coordinate System Metrics dialog for the Points drawing and entering the same values for the drawing's metrics.
An easier, less error-prone way to ensure the Points drawing and the 1475 PM image both use the same metrics, is to save the coordinate system used by the 1475 PM image as a favorite that overrides metrics, and to then apply that favorite to the Points drawing with one click. Much easier!
We will use the procedure given in the Favorite Coordinate Systems topic to add a new favorite that overrides metrics.
We click the 1475 PM tab in the map to move the focus onto that tab.
In the Info pane,
we click the coordinate system picker button and choose Repair
Initial Coordinate System. and then in the resulting
sub-menu we choose Favorites.
The Favorites dialog opens up with the coordinate system used by the 1475 PM image already loaded in the lower pane as a potential new favorite. Ctrl-click that coordinate system to select it.
We press Add to Favorites.
That adds the new favorite, with 2 appended to the name to differentiate it from the default Pseudo-Mercator coordinate system already in the list.
We double-click
into the metrics box and choose override
metrics. This is a key step, because it will apply
the metrics the image uses whenever this favorite is applied.
Press OK. That creates a new version of the default Pseudo-Mercator coordinate system, a version that applies the metrics values used by the 1475 PM image instead of whatever metrics were previously used.
In the map, we click on the Points layer to make it the active layer.
In the Component pane,
we click the coordinate system picker button and choose Repair
Coordinate System. We then choose our new, override
metrics, favorite. In addition to the slightly different name, it
is marked with a hash # character
to indicate it overrides metrics.
Instantly, the points appear in the Points layer, since the drawing now uses a correct coordinate system. We have not styled the drawing, so the points appear in default, gray color and round point symbol.
We use Style to format the points to a brighter, more visible color and symbol. To take a closer look, we right-click and drag to zoom box around the upper area in the display.
Not surprisingly, the three pixels within the area with the highest elevation values are all near each other, in this case lined up on a ridge.
We can create labels using the height value saved to the points in the Points drawing's table. In the illustration above, the middle point's label does not appear because it is automatically clipped by the too-close proximity of the other labels.
Zoomed out, the labels appear for the points. Note that in cases where areas overlap terrain that slopes upward, the highest pixels are likely to be found on the edge of the area. Only in cases where the area surrounds a mountain peak or similar local maximum will the highest pixels appear within the area.
How did we get labels without many digits after the decimal point? - Looking at the last illustrations above, the height values do not show many digits after a decimal point, but are round integer values. We accomplished that by using the Round SQL function in the text pattern for the label. If our label text simply used [Value] the result would have many digits after the decimal point. If the label text uses Round( [Value] ) then the height value will be rounded to a whole integer number. See the Label Text topic.
using the Schema dialog to add an integer field to the table, and then using the Transform pane to copy the floating point Value numbers into the integer field, which automatically rounds them to the nearest integer value. We then created labels from the integer field.
Adapting the query - We can use the above query in our own work by simply substituting the name of our image for 1475 PM, and the name of our drawing for Drawing. Ctrl-h for search and replace, a standard Windows shortcut, works within the Command Window.
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:
-- take all pixels, return n highest
FUNCTION highestCoords(@t TABLE, @n INT32) TABLE AS (
SELECT * FROM @t ORDER BY [value] DESC FETCH @n
) END;
-- select 3 highest pixels under each object in the drawing.
SELECT SPLIT CALL highestCoords(
CALL TileGeomToValues([1475 PM], [Geom]), 3
) INTO [Points2 Table] FROM [Drawing];
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.
Style: Contouring using Colors
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. Easy!
Example: Import DDF SDTS DLG Vector File
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: 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.