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.  

 

Manifold's SQL functions for Kriging work in two steps: first, the TileInterpolateKriging or TileInterpolateKrigingPar function takes a drawing along with options and returns an interpolation option.  Second, the specific interpolation function desired takes the interpolation object and creates tiles.  

 

In this example we take a raster image that provides terrain elevation value for a region of the Alps and we create a vector drawing of contour lines.   This first part of the example simply makes it clear what vector data we are using.   Next, we write a query that takes the vector data in the form of contour lines and then rasterizes it using Kriging.      We do the example twice, first for a smaller data set and next for a larger data set, to show how the results table reports more tiles for the larger data set.

 

Example: Smaller Data

 

 

eg_sql_kriging01_01.png

 

 

We start with a single tile of SRTM terrain elevation data, seen above styled with a palette.   The data is one tile of the SRTM data used in the Example: Merge Images topic, showing terrain elevations in the Italian Alps near Bolzano, Italy.

 

eg_sql_kriging01_01a.png

 

Opening the image's table we see the data is a data type of int16, a typical type for terrain elevation data.

 

 

eg_sql_kriging01_02.png

 

 

Zooming further into the Bolzano region we can see some of the pixels in the data set are missing.  These are locations where the Space Shuttle instrumentation did not acquire an elevation.   

Create Contour Lines

The point of this example is to show how to create a raster image from a vector drawing using Kriging.  We will first create a vector drawing, contour lines, from our alps terrain elevation raster.   This will make it clear exactly where the vector drawing came from and then later, when we use Kriging to create a raster from that vector drawing, we can compare what the result of Kriging is with the original alps raster.

 

eg_sql_kriging01_03.png

 

With the focus on the alps layer we use the Transform panel and choose the Contour Lines template.  We create contour lines in steps of 250 meters form a minimum height of 0 to a maximum height of 3500.    See the Example: Contour Areas and Contour Lines topic for a step-by-step procedure.  

 

 

eg_sql_kriging01_04.png

 

 

As soon as we enter the minimum, maximum and step values, the Transform template previews what will be created using blue preview color.   We press Add Component to create the contour lines drawing, which we call alps_contours.

 

 

eg_sql_kriging01_05.png

 

 

Showing the alps_contours drawing in the map, with the lines styled in green color, against a black map background shows the very detailed contour lines that have been created.

 

eg_sql_kriging01_05a.png

 

Opening the alps_contours Table we see that each line object has a Value that gives the height of that contour line.  Even with this relatively small data set there are over 5320 contour lines.

 

Use Kriging to Create a Raster Image from Contour Lines in a Drawing

Now that we have a drawing with contour lines, we can write an SQL query that takes that drawing and uses Kriging to interpolate a raster image.    

 

To simplify the query, instead of creating a destination image from scratch in the query, we will first copy and paste the alps image and its table to create a destination image in the Project called alps_kriging and its table.   This saves us from having to write a paragraph or so of SQL that first creates a destination image, with all characteristics such as coordinate system correctly set up, ready to be filled with tiles created by the Kriging process.  Instead of doing that in SQL, we just manually use copy and paste to create a destination alps_kriging image and then we wipe it clean with a simple DELETE FROM statement.

 

 

eg_sql_kriging01_06.png

 

We run the following query, which produces 310 tiles as seen in the Result table above.

 

-- put height into geoms

UPDATE [alps_contours] SET [geom] = GeomSetZ([geom], [Value]);

-- drop old tiles

DELETE FROM [alps_kriging];

-- add new tiles

INSERT INTO [alps_kriging] ([x], [y], [tile])

SELECT [xy].[x], [xy].[y],

  CASTV(TileInterpolate(CALL TileInterpolateKriging([alps_contours], -1, -1, -1), [xy].[rect]) AS INT16)

FROM CALL ValueSequenceTileXY(

  CAST (ComponentProperty([alps_kriging], 'Rect') AS FLOAT64X4), -- rect in pixels

  CAST (ComponentProperty([alps_kriging Tiles], 'FieldTileSize.Tile') AS FLOAT64X2), -- tile size

  false -- tiles with any part inside of rect

) AS [xy];

-- update intermediate levels

EXECUTE CALL TileUpdatePyramids([alps_kriging]);

 

A few words on how the query works, step by step:

 

-- put height into geoms

UPDATE [alps_contours] SET [geom] = GeomSetZ([geom], [Value]);

 

Kriging works by extracting Z values that are part of the object's geometry.  However, contour lines are created with only X,Y coordinates in their geometry and the Z value is given by a Value attribute field as we saw in the table above.  This first part of the query simply takes whatever is the Z value in the Value field and stuffs it into the geometry to make the geom an X,Y,Z geom.  

 

-- drop old tiles

DELETE FROM [alps_kriging];

 

Rather than create a new image within the SQL query we simply copied and pasted the original raster image to provide a suitable target image.  This part of the query wipes that destination clean by removing any tiles it may contain.

 

-- add new tiles

INSERT INTO [alps_kriging] ([x], [y], [tile])

SELECT [xy].[x], [xy].[y],

  CASTV(TileInterpolate(CALL TileInterpolateKriging([alps_contours], -1, -1, -1), [xy].[rect]) AS INT16)

FROM CALL ValueSequenceTileXY(

  CAST (ComponentProperty([alps_kriging], 'Rect') AS FLOAT64X4), -- rect in pixels

  CAST (ComponentProperty([alps_kriging Tiles], 'FieldTileSize.Tile') AS FLOAT64X2), -- tile size

  false -- tiles with any part inside of rect

) AS [xy];

 

After cleaning out old tiles from the destination alps_kriging image, the INSERT repopulates what will be the interpolated image.   

 

The call to ValueSequenceTileXY in the FROM section of the INSERT produces the X, the Y, and the Rect values for tiles to interpolate using the Rect of the destination image and the tile size, both pulled out of the destination image's properties using calls to ComponentProperty.  

 

The tiles are produced using a call to TileInterpolate in the SELECT list. The interpolation object is created using TileInterpolateKriging, with all parameters set for automatic chocies, that is using Voronoi neighbors and automatic choice of model.  The produced tiles are converted to INT16 using CASTV, because that is the data type used for tiles by the destination image.

 

 

eg_sql_kriging01_07.png

 

 

Because we created the destination alps_kriging image by copying and pasting the original alps image, it still has within its properties a StylePixel property that contains the same styling the original alps image used.   When the query runs and populates the alps_kriging image with interpolated tiles, those tiles will be displayed using the same styling.   In the illustration above we see the interpolated alps_kriging image zoomed into the same view near Bolzano.

 

 

eg_sql_kriging01_08.png

 

 

Zoomed all the way out, we see that the alps_kriging image covers the same region as the original alps image.   When zoomed far out, it looks like the round-trip process of first vectorizing the alps image to create a drawing of contour lines and then rasterizing that drawing of contour lines using Kriging to create a raster image, produced a raster terrain elevation image remarkably close to the original.

 

eg_sql_kriging01_09.png eg_sql_kriging01_10.png

 

When we zoom closer and compare the alps_kriging image, at left, to the original alps data set, at right, we can see that the round trip process is, of course, lossy.  Vectorizing the relatively smooth alps terrain elevation raster into quantized contour lines at 250 meter intervals necessarily throws away the actual data that falls within those 250 meter contours.    Reconstructing the raster from quantized vector heights at 250 meters gives the somewhat terraced effect seen at left.  It is a tribute to the interpolative power of Kriging that the interpolated raster looks as good as it does.   Note also that it fairly convincingly fills in regions of missing data in the original.

A Larger Data Set

It can take dramatically longer to interpolate bigger data sets than just a single SRTM tile.    We can repeat the same process seen above using a merged raster created from four SRTM tiles.

 

eg_sql_kriging01_13.png

 

Starting with the terrain elevation seen above we create contour lines.

 

eg_sql_kriging01_14.png

 

We will now use the contour lines to create a raster, using the same query as before.

 

 

eg_sql_kriging01_11.png

 

Running the query we now get a result of 1170 tiles instead of only 310 tiles as before.

 

eg_sql_kriging01_12.png

 

The resulting alps_kriging image looks remarkably good when zoomed out as seen above.

 

eg_sql_kriging01_16.png

 

We can zoom in to see how it compares to the original.   The original alps raster is seen above.

 

eg_sql_kriging01_15.png

 

The contour lines for the zoomed in view are seen above.

 

eg_sql_kriging01_17.png

 

The alps_kriging raster created using Kriging from those contour lines is seen above.   The terraced effect from using vector data constrained to 250 meter intervals is clearly seen.

 

We can attempt to reduce the terracing effect by using the TileInterpolateKrigingMedianPolish function, applying median-polish as part of the interpolation.   The query would be virtually identical:

 

-- put height into geoms

UPDATE [alps_contours] SET [geom] = GeomSetZ([geom], [Value]);

-- drop old tiles

DELETE FROM [alps_kriging];

-- add new tiles

INSERT INTO [alps_kriging] ([x], [y], [tile])

SELECT [xy].[x], [xy].[y],

  CASTV(TileInterpolate(CALL TileInterpolateKrigingMedianPolish([alps_contours], -1, -1, -1), [xy].[rect]) AS INT16)

FROM CALL ValueSequenceTileXY(

  CAST (ComponentProperty([alps_kriging], 'Rect') AS FLOAT64X4), -- rect in pixels

  CAST (ComponentProperty([alps_kriging Tiles], 'FieldTileSize.Tile') AS FLOAT64X2), -- tile size

  false -- tiles with any part inside of rect

) AS [xy];

-- update intermediate levels

EXECUTE CALL TileUpdatePyramids([alps_kriging]);

 

eg_sql_kriging01_18.png  eg_sql_kriging01_09.png

 

At left is the median-polish result compared to Kriging without median-polish.   Whether it makes sense or not to use median-polish in a given situation depends on the data and the objectives.  At the cost of introducing spurious artifacts it can make the interpolation smoother.

 

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.  

 

Krige invented his method to help develop mines for gold-bearing ore bodies that were deep underground.  Drilling sampling boreholes was very costly, so estimates about the shape and location of such ore bodies had to be made on the basis of a sparse number of scattered boreholes.  Kriging, as his method became known, allowed Krige to make better estimates to interpolate the likely layout of ore bodies in between scattered boreholes.

 

The generalization of the technique allows interpolation of values from values at scattered, known locations.   For example, to create a raster surface from scattered point values or contour lines we do not know exactly what the terrain may be in between contour lines or in between specific, scattered point locations, but Kriging provides a rational way of generating a likely range of pixel heights for locations in between known values.   Kriging has thus become a very useful technique for creating raster surfaces from vector drawings.  

 

The term as generally used in GIS may or may not be capitalized.   It is generally pronounced in English as "Kree - ging" with a hard "g" as in the name of the tree, "ginko" or the word "geisha".

 

img_krige.png

 

Danie Krige (1919 - 2013)

 

A mining engineer and statistician, Danie Krige over the course of a long life revolutionized the spatial analysis of ore bodies while serving his country, South Africa, in the development of the nation's mineral wealth.  

 

Born in Bothaville, Free State, he graduated at age 19 from University of the Witwatersrand and in the 1930's embarked on a career of surveying and ore valuation of gold mines for Anglovaal, a major mining company in South Africa.   In the 1940's he moved into the department of the Government Mining Engineer, working closely with industry to help develop the country's mineral wealth.   In the 1930's and 1940's Krige developed his application of statistics to assessing the layout and valuation of ore bodies based on sparse data from boreholes, describing his method in a key paper in 1951 and in his Master's thesis.   He returned to Anglovaal in 1951 to apply Kriging for the first in the routine evaluation of large gold mines and reserves.  

 

In addition to improving the efficiency of economically risky development of gold ore bodies deep underground, Krige served South Africa in the 1940's by designing a pricing formula for uranium, to establish profitable trade in uranium between South Africa and the United States and the United Kingdom.  In the 1950's, Krige's analyses helped establish a solid foundation for South Africa's uranium industry.    

 

In 1981 Krige became a professor at the University of the Witwatersrand,  educating the next generation of South African experts in geostatistics and mining economics while continuing his research on the application of geostatistics to ore evaluation, publishing over 90 technical papers worldwide.    In addition to his doctorate from the University of the Witwatersrand, Krige received numerous awards and honors from institutions throughout the world, including both the Moscow State Mining University and the US National Academy of Engineering.   In English his last name is pronounced "Kreeg" with a hard "g" at the end.

 

img_matheron.png

 

Georges Matheron (1930 - 2000)

 

A French mathematician and geologist, George Matheron coined the term geostatistics to name the new science of applied statistics for estimating resources.   He also coined the term Kriging in honor of Danie Krige's seminal work upon which Matheron built after Krige's initial papers were translated into French in 1955. Matheron also invented the analytic approach he named mathematical morphology which is now widely used in image processing and other fields beyond geology.   Matheron founded  the Centre de Géostatistique et de Morphologie Mathématique at the Paris School of Mines in Fontainebleau, France, a major administrative achievement in addition to his four decades of inventive and productive work in applied mathematics.  It seems appropriate that such a first-rate mathematician as Matheron would have a family name that begins with "Math."

 

See Also

Images

 

Tables

 

Command Window

 

Data Types

 

SQL Functions

 

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.

 

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: Contour Areas and Contour Lines - In this example we use the Contour Areas transform template in the Transform panel 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.

 

Example: Vector to Raster using Kriging - Using the Interpolate, Kriging Transform panel template, 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.