Example: Vector to Raster using Kriging

Using the Interpolate : Kriging Transform pane operation we take a vector drawing of contour lines where each line has a Height attribute and we create a raster image that is a terrain elevation surface.   This example, uses a point and click Transform template to accomplish the same task as shown using an SQL query in the SQL Example: Kriging topic.   The data set and the project contents are the same as in that topic as well.

 

A typical task in GIS is to take a drawing with objects, such as a grid of points, a set of contour lines, or a set of areas, each of which has a height value associated with it, and to create a smooth image that represents an elevation surface, interpolating in between locations where we have precise height data for a vector object in the drawing, to estimate what the heights of pixels should be in between locations where we have data.  

 

One such method for interpolating is Kriging, a method which in many applications is especially useful for interpolating data in between sparsely distributed "known" locations.   Kriging variations supported by Manifold include the Interpolate: Kriging with median polish and Interpolate: Kriging with regression transforms, as documented in the Transform - Geometry: Interpolate topic.

 

In this example, we use the Interpolate: Kriging operation to create a surface that is smooth and continuous, even though the data from which we start is a set of contour lines with no data in between the lines.   Creating a raster using Kriging in GIS is often done based on data sets of points, but the Manifold Kriging transform can use lines, such as contour lines, as well as points.

 

To fit into this documentation, illustrations show a small Manifold desktop, with only a few panes, docked to the right side.  In real life we use a much larger Manifold desktop, and more panes would be turned on, with panes docked to the left or to the right, or undocked, as we prefer.   Right-click a pane's tab to change where it is docked.  Manifold will remember that new arrangement for our next session.

 

 

 

 

To show exactly what the contour lines represent, we create contour lines from a raster image that represents a terrain elevation surface.  The sample image is space shuttle SRTM data that shows terrain elevations in the Italian Alps, with a zoomed-in view around Bolzano shown above.   The original image is called alps and is seen above, styled with a useful palette and hill shaded to provide a sense of depth.   Black regions indicate missing values, where there is no data.   In such places the pixels are transparent and allow the black background to show through.

 

 

 

 

Following the procedure in the Example: Contour Areas and Contour Lines topic, using the alps image we have created a drawing containing contour lines from 0 to 3500 meters with a step of 250 meters.    The image above shows the contour drawing lines overlaid in a map above the alps image.   We Style the contour lines using orange-yellow color for line stroke color, with a right-side accessory line in dark gray that is .3 point to the right of the main line.  That provides a nice pseudo-depth effect to the lines.

 

 

 

 

In the illustration above we have turned off the alps image to reduce visual clutter, showing the yellow contour lines in the alps_contour drawing clearly against a black background.  We have opened the drawing's table, to see how each line object has a Height attribute that gives the height of that line.  The default contour line creation template creates a height field called Value.  For clarity, we have used the Schema dialog to rename that field Height.  Each line provides precise data: exactly on the line segment the height of the terrain is the Height value.  In between the contour lines there is no data.   The first and last records in the table have NULL values, since there are no contours at 0 height and at 3500 height.

 

This is not as sparse a vector data set as if there were just a scattering of points with a height for each point and no data in between.   All the same we will use Kriging interpolation since Kriging is quite good at creating smooth interpolations both from sparse data and data that is not sparse.

 

With the focus on the map window, in the Transform pane we choose the alps_contours drawing layer and the Geom field within that layer.  We double-click the Interpolate template to launch it.   

 

The Transform pane is shown below undocked, so we can resize it to a taller size to provide ample room for all of the many controls that are used for Kriging.

 

 

In the Interpolate template, from the pull down menu for Interpolation we choose Kriging as the interpolation method.  

 

 

For the  Z parameter, from the pull down menu we choose the Height field.

 

The alps image from which the contours were created was in Latitude / Longitude coordinate system, so the alps_contours drawing created from it is also in Latitude / Longitude and thus the default units are  degrees.   For the Resolution, we enter 0.000278.   That many degrees is approximately 30 meters at the latitude in view, approximately 46.5 degrees North.   If we wanted a less sloppy approach, we would have re-projected the alps image into some reasonable, meter-based projection, like Lambert Conformal Conic or even Pseudo-Mercator.  We could then work with linear units like meters directly, instead of messing about with estimates of fractions of a degree.

 

The Result destination for Interpolation operations is always New Table, meaning to create a new table and a new image to visualize the tiles in the new table.   We enter Kriged alps as the name of the New image, and Kriged alps Tiles as the name of the New table.   We can use whatever names we want, but it makes sense to use names that will later remind us what those components are supposed to be.

 

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.

 

To see a preview, press the Preview button.   Previews for Kriging take as long as actually doing the transform, so many people will skip previews for such interpolations, but previews can still be useful in such cases to check workflow before proceeding.  

 

 

 

 

A preview appears in blue preview colors.  The preview layer is a virtual layer that appears above all other layers and contents in the map.  A blue preview caption bar appears at the top of the window, giving the name of the template being previewed.  We can right-click the caption bar for a menu of choices to control the preview display.

 

 

 

 

For example, if we would like to compare the preview to a view of layers below the preview, we can right-click the caption bar and choose 75% opacity to render the preview with partial opacity, so what is below can be seen.  

 

 

 

 

We can right-click the caption bar and choose Left or Right to show the preview only on the left or right side of the window, dragging the thin blue vertical separator line left or right to adjust the width of the preview.

 

If we like what we see, to apply the template we press Transform.  

 

Manifold thinks a bit, and a new Kriged alps image and table appear in the Project pane.

 

 

 

 

The result is that an image called Kriged alps is created, along with the image's table.   We drag and drop the new image into the map, as a layer below the alps_contours drawing layer.  The image appears using default grayscale formatting, since it has not yet been Styled.  

 

 

 

 

So we can better compare this new image we have created by Kriging from the contour lines drawing to the original alps image from which those contour lines were created, we would like to style the new image using exactly the same palette and shading settings as the alps image uses.  That is easy to do by copying the StylePixel property from the alps image to the new image.   We right-click onto the alps image and choose Properties.

 

 

In the Properties dialog for the alps image we right-click onto the value cell for the StylePixel property and we choose Copy.  That copies the long JSON text string that specifies the styling of the alps image to the Windows clipboard.

 

 

Next, we right-click onto the Kriged alps image and choose Properties.  In the properties dialog for the Kriged alps image we click into the first * cell and create a new StylePixel property, and then we can right-click into the value cell for that property and Paste the copied JSON string.   Click OK and we have added a StylePixel property to the Kriged alps image which is identical to that used in the alps image.

 

 

 

 

Like magic, the Kriged alps image immediately changes from default, grayscale formatting to using the same palette and shading as the alps image.

 

 

 

 

We can turn off the alps_contours layer for a clearer look at the new image we have created by interpolating from objects in the contour lines drawing using Kriging.

 

 

Comparing the Kriged alps image at left with the original alps image at right we can see how missing data locations in the alps image have been filled in with interpolated values in the Kriged alps image.   We can also see that the Kriged alps image is lower resolution, that is, less detailed.   That makes sense, since it was created from height data in lines that were 250 meters  apart, a distance that in most locations is further apart than more detailed SRTM elevation data.

 

The Kriged alps image also shows a terraced effect, also expected since each contour line represents the data as if it were made up of a series of stepped plateaus or terraces.   If the contour lines were at closer intervals the resulting interpolated surface would be smoother and have better detail.   

The DEST Algorithm

Another way to reduce terracing effects is to use the Interpolate : triangulation with segments transform with the Remove flat areas (DEST) option checked.   The DEST (Determination of Earth Surface Structures) algorithm helps to remove or reduce terracing effects when contouring from contour lines or other regular data.  DEST is a triangulation interpolation method that uses a modified Delaunay triangulation.

 

 

We can choose triangulation with segments as the Interpolation method in the Transform pane, and then do the transform once with Remove flat areas (DEST) box unchecked, and then again with the box checked, saving the results to suitably named new components.

 

We then copy the StylePixel property from the alps drawing to the resulting images, so they are formatted the same for easy comparison.

 

 

In the illustration at left above we have used the Interpolate : triangulation with segments template, without the Remove flat areas (DEST) option being checked.   The terrain elevation surface image was created from the alps_contours drawing of relatively coarse contour lines at 250 meter elevation intervals, from which the template creates an unsmooth surface that shows terracing effects despite the use of triangulation interpolation.      

 

The illustration at right above shows the result of exactly the same triangulation with the Remove flat areas (DEST) option being checked.  Applying the DEST algorithm provides some variation even in regions where relatively large flat surfaces appear without DEST.

Notes

Kriging - Kriging is a geostatistical process invented by Danie Krige, a South African statistician and mining engineer, that uses data from locations where good data is available to make estimates, by interpolation, for locations where data is not available.  See details in the Notes to the SQL Example: Kriging topic.

 

Resolution - The Resolution parameter sets the siize of resulting pixels, expressed as a multiplier of drawing coordinate system units.  Do not use a Resolution of 1 when drawings are in Latitude / Longitude coordinate system. That creates pixels which are 1 degree in size, many kilometers in size in most parts of the world.     If drawings are in a coordinate system that uses meters, a Resolution of 30 creates pixels that are 30 meters by 30 meters in size.  To keep units straight, it is best to do interpolations using drawings that have been projected into coordinate systems which use linear units, such as meters or feet, and not  Latitude / Longitude, which uses degrees as a unit.

Videos

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!

See Also

Drawings

 

Images

 

Tables

 

Transform Pane

 

Example: How Images use Tiles from Tables - An example showing how an image is made up from data stored in a table in tiles.

 

Example: Create Two Images From One Table - More than one image can show data from the same table, including from the same tile field.

 

Example: An Image using Computed Fields in a Table - How an image can be created from tiles where the data for the tiles is taken from a field that is computed on the fly.

 

SQL Example: Re-tile an Image using a Different Tile Size - Starting with an image that uses a tile size of 128 x 128 pixels this SQL example creates a copy of the image using 500 x 500 pixel tiles.

 

SQL Example: Kriging - We use SQL functions to create a raster terrain elevation image from vector contour lines in a drawing, using SQL functions for Kriging interpolation.

 

Example: Merge Images - A step-by-step example using the Merge Images command showing how to merge dozens of images showing SRTM terrain elevation data into one image, with various tricks for faster workflow as an experienced Manifold user would do the job.  After creating the new image we style it with a palette and use hill shading to better show terrain elevation.

 

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 hillshading.  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: Contour Areas and Contour Lines - In this example we use the Contour Areas transform template in the Transform pane for images to create a drawing with vector areas showing height contours at desired altitude steps.   We color the areas using the attribute fields automatically created by the template.  Next, we apply a similar procedure using the Contour Lines transform template to create a drawing with vector lines showing height contours at the desired intervals.