Dev Direct Solution Center

How To Use GRID to Enhance Your GIS Development with Map Suite

Introduction

This solution shows how to use GRID, a really neat way of taking point data readings and interpolating them to produce an attractive visual representation of that data on a map. Utilizes the Map Suite Desktop .NET mapping control.

Detail

GRID can come handy in numerous situations where you need to create custom map displays, perform interpolation or do spatial analysis. We hope that by reading this article you will see more clearly how to apply GRID to your own industry and reap the benefits it provides.

For the purposes of this article we have used the Map Suite Desktop Edition to show the functionality of GRID, but please note that GRID functionality is also available in the Map Suite Web and Engine editions.

Note: This article is meant to be used as a tutorial for using GRID as a GIS technique, not as a definitive tool for soil analysis. The soil quality examples included here are fictitious, and the soil quality calculations based on pH and Magnesium levels do not come from an authoritative source in agriculture.

GRID Explained

A GRID is a raster format that defines a geographic space as an array of equally sized squares (cells) arranged in rows and columns. Each cell stores a numeric value that represents a geographic attribute (such as elevation, surface slope, soil pH, etc.) for that unit of space. Each grid cell is referenced by its x, y coordinate location.

File Format

The format of the GRID file is ASCII and contains two main parts, the header and the body. The header consists of the first few lines for references of the grid such as the number of columns, the number of rows, the x and y coordinate of the lower left of the grid, the cell size and the no data value (which will be used in the body of the GRID file to represent a null value). The body contains the cell values listed in the order they would naturally appear from left to right and from top to bottom.

Example of ASCII GRID Format:

ncols         4
nrows         6
xllcorner     0.0
yllcorner     0.0
cellsize      50.0
NODATA_value  -9999
-9999 -9999 5 2
-9999 20 100 36
3 8 35 10
32 42 50 6
88 75 27 9
13 5 1 -9999
etc...

Getting Started

Now that we have covered the theory of the GRID, we are going to use a very concrete example to show how useful GRID can be in practical GIS applications.

Let's say that we have a point-based shapefile representing some sample points where the pH (acidity level) of the soil was measured on a field.

Fig 1: Point-based shapefile of sample pH levels.

Fig 2: Map legend.

Now with those few sample points, we want the entire field to be covered with a pH value. This is when GRID comes in handy. Using an interpolation method, IDW (Inverse Distance Weighing), we will create a GRID and then display it on the map to show clearly the different pH values throughout the field.

Creating the GRID

Before creating the GRID, we need to decide the extent of the grid, the size of the cells and the value for no data. In our case, we are going to take a slightly larger extent of the field and we are going to have a cell size of 0.0001 (we are in decimal degrees, so 0.0001 degrees are about 9.6 meters). The no data value is -9999. This means that any cell that is outside the field (where the pH value is irrelevant), we will assign it the -9999 value.

In the Click event of a Button, we will write the logic to create the GRID.

Private Sub btnIDW_Click(ByVal sender As System.Object, ByVal e As System.EventArgs) Handles btnIDW.Click
        Dim gridExtent As New StraightRectangle(Map1.CurrentExtent)
        Dim cellSize As Double = 0.0001
        Dim NoDataValue As Double = -9999

        AddHandler GridLayer.CalculateCellValue, AddressOf CalculateCellValueHandlerIDW
        Dim Result As Boolean = GridLayer.CreateGridFile("..\PHIDW_2_200_A.grd", gridExtent, cellSize, NoDataValue)

        If Result = True Then
            MessageBox.Show("Grid File saved successfully.", "Information", MessageBoxButtons.OK, MessageBoxIcon.Information)
        Else
            MessageBox.Show("Saving Grid File Failed.", "Information", MessageBoxButtons.OK, MessageBoxIcon.Information)
        End I
    End Sub

At the creation of the GRID, the event CalculateCellValue is going to be raised for each cell to assign a value. So we are going to write the interpolation algorithm in the handler CalculateCellValueHandlerIDW. In our example, we are using IDW interpolation, or Inverse Distance Weighing. This is a common method for interpolation that assigns values to unknown locations based on some scattered set of known points – which is exactly the case in our example. This interpolation has two variables: The number of points to take into account (based on distance) and the power parameter. Playing with those two variables is going to give a result with smoother or sharper peaks.

Notes: We are using IDW interpolation in this example, but you are free to use any other type of interpolation or algorithm to calculate the cell value. For more information on interpolation, check the following site: http://en.wikipedia.org/wiki/Multivariate_interpolation

Private Sub CalculateCellValueHandlerIDW(ByVal CellExtent As StraightRectangle, ByRef CellValue As Double)
        Dim CellULx As Double = CellExtent.UpperLeftPoint.X
        Dim CellULy As Double = CellExtent.UpperLeftPoint.Y
        Dim CellLRx As Double = CellExtent.LowerRightPoint.X
        Dim CellLRy As Double = CellExtent.LowerRightPoint.Y

        Dim xp As Double = CellExtent.Center.X
        Dim yp As Double = CellExtent.Center.Y
        Dim CellCenterPointShape As New PointShape(xp, yp)
        Dim CenterRecords() As Integer = Map1.Layers(0).QueryByDistance(CellExtent, 15, MapLengthUnits.metres, Map1.MapUnit)

        If CenterRecords.GetLength(0) = 1 Then

            Dim NearestRecords() As Integer = Map1.Layers(1).QueryByDistance(CellCenterPointShape, 200, MapLengthUnits.metres, Map1.MapUnit)
            If NearestRecords.GetLength(0) > 0 Then
                Dim NearestPointShapes() As PointShape = CType(Map1.Layers(1).GetShapes(NearestRecords), PointShape())
                Dim NearestValues(NearestRecords.GetLength(0) - 1) As String
                For i As Integer = 0 To NearestRecords.GetLength(0) - 1
                    NearestValues(i) = Map1.Layers(1).DataQuery(NearestRecords(i), "pH")
                Next i
                'This is the algorithm for IDW interpolation
                Dim Power As Double = 2
                Dim Zi, Zx1, Zx2, Zx, Di As Double
                For i As Integer = 0 To NearestValues.GetLength(0) - 1
                    Zi = CDbl(NearestValues(i))
                    Di = CellCenterPointShape.DistanceTo(NearestPointShapes(i), Map1.MapUnit, MapLengthUnits.metres, DistanceTypes.Shortest)
                    Zx1 = Zx1 + (Math.Pow((1 / Di), Power) * (Zi))
                    Zx2 = Zx2 + ((Math.Pow((1 / Di), Power)))
                Next
                Zx = Zx1 / Zx2
                CellValue = Zx
            Else
                CellValue = -9999
            End If
        Else
            CellValue = -9999
        End If
    End Sub

In our example, we are using a distance of 200 meters to get the sample points within that distance and we are using a power of 2. By specifying a higher power, more emphasis is placed on the nearest points, giving a smoother result. By specifying a lower power, more influence is placed on the points that are farther away, giving a sharper result.

Displaying the GRID

Now that we have the GRID created, we need to display it. This article will demonstrate two different ways of displaying the GRID. For the first way, we will apply a specific color to any cell belonging to a certain class. The second method uses a gradient to go from one color to another. The purpose is to show different ways one GRID can be displayed.

Using ClassBreaks

We are going to display the cells of the GRID in the same way as we did for the original shapefile of sample points. Note that you might want to remove the two shapefiles from the Map so that they do not appear on top of the GRID.

Dim STR_GridFile As String = "..\SoilQuality.grd"
        Dim gridlayer As GridLayer = New GridLayer(STR_GridFile)
        gridlayer.LowerThreshold = 0
        gridlayer.UpperThreshold = Double.MaxValue
        gridlayer.ThresholdUnit = ThresholdUnits.miles
        gridlayer.IsNonValueTransparent = True
        Dim gridranges As RangeCollection = New RangeCollection()
        gridranges.Add(New Range(7.54, 7.9, Color.FromArgb(0, 255, 0)))
        gridranges.Add(New Range(7.21, 7.54, Color.FromArgb(128, 255, 128)))
        gridranges.Add(New Range(7.15, 7.21, Color.FromArgb(224, 251, 132)))
        gridranges.Add(New Range(7.08, 7.15, Color.FromArgb(225, 255, 0)))
        gridranges.Add(New Range(7.0, 7.08, Color.FromArgb(245, 210, 10)))
        gridranges.Add(New Range(6.83, 7.0, Color.FromArgb(255, 128, 0)))
        gridranges.Add(New Range(6.2, 6.83, Color.FromArgb(255, 0, 0)))
        gridlayer.Ranges = gridranges

        Map1.GridLayers.Add(gridlayer)

Fig 3: GRID resulting from sample points interpolation displayed with Class Breaks.

Using Gradient

We are going to display the cells of the GRID with a gradient going progressively from red (for low pH) to green (for high pH). Note that in addition to Map, Layer (the vector shapefile layer) has a GRIDLayers collection so that GRID files can be displayed in the desired order. Here we are displaying the GRID file between the Field Layer and the sample points Layer.

Dim STR_GridFile As String = "..\SoilQuality.grd"
        Dim gridlayer As GridLayer = New GridLayer(STR_GridFile)
        gridlayer.LowerThreshold = 0
        gridlayer.UpperThreshold = Double.MaxValue
        gridlayer.ThresholdUnit = ThresholdUnits.miles
        gridlayer.IsNonValueTransparent = True
        Dim gridranges As RangeCollection = New RangeCollection()

        Dim Alpha As Integer = 255
        gridranges.Add(New GradientRange(6.2, Color.FromArgb(Alpha, Color.FromArgb(255, 0, 0)), 7.9, Color.FromArgb(Alpha, Color.FromArgb(0, 255, 0))))
        gridlayer.Ranges = gridranges
     Map1.Layers(0).GridLayers.Add(gridlayer)

Fig 4: GRID resulting from sample points interpolation displayed with gradient. The Layer with the sample points appears on top.

Multi GRID operations

One of the big advantages of GRID is to be able to create a new GRID by analyzing and performing calculations on the cell values of an existing GRID, or by combining two existing GRIDs to create a third one. For our next example, we have a second GRID that represents Magnesium levels throughout the field from 40 to 95 ppm. We will combine the pH value GRID and the Magnesium level GRID to create a third, composite GRID.

Displaying the GRID for Magnesium

We are going to display the GRID Magnesium_Level using a gradient. (This previously-created GRID is supplied to us for the purpose of this example.)

Dim STR_GridFile As String = "..\Magnesium_Level.grd"
        Dim gridlayer As GridLayer = New GridLayer(STR_GridFile)
        gridlayer.LowerThreshold = 0
        gridlayer.UpperThreshold = Double.MaxValue
        gridlayer.ThresholdUnit = ThresholdUnits.miles
        gridlayer.IsNonValueTransparent = True
        Dim gridranges As RangeCollection = New RangeCollection()
        gridranges.Add(New GradientRange(40, Color.FromArgb(254, 215, 186), 90,    Color.FromArgb(105, 46, 1)))
        gridlayer.Ranges = gridranges
        Map1.Layers(0).GridLayers.Add(gridlayer)

Fig 5: GRID of Magnesium level in PPM from 40 to 95. The darker color the higher the value.

Creating a third GRID from pH GRID and Magnesium GRID

For our example crop, the combination of a pH value between 7.0 and 7.5 and a high magnesium level represents the best soil. The resulting GRID will represent soil quality for that crop by combining the two GRIDs. The cell values will be calculated according to a common formula used by agricultural engineers. We attribute a factor for each pH value and Magnesium level, then we combine the two factors to get the soil quality for that crop.

To determine the soil quality, we will combine the pH factor with the Mg factor. For example, if the pH is 6.67 and the Mg level is 55, the soil quality index will be 28 (7 * 4). The soil quality index goes from a minimum of 5 to a maximum of 100.

Private Sub btnSoilQuality_Click(ByVal sender As System.Object, ByVal e As System.EventArgs) Handles btnSoilQuality.Click
        Dim gridExtent As New StraightRectangle(Map1.Layers(0).GridLayers(0).Extent)
        Dim cellSize As Double = 0.0001
        Dim NoDataValue As Double = -9999

        AddHandler GridLayer.CalculateCellValue, AddressOf CalculateCellValueSoilQuality

        Dim Result As Boolean = GridLayer.CreateGridFile("..\SoilQuality.grd", gridExtent, cellSize, NoDataValue)

        If Result = True Then
            MessageBox.Show("Grid File saved successfully.", "Information", MessageBoxButtons.OK, MessageBoxIcon.Information)
        Else
            MessageBox.Show("Saving Grid File Failed.", "Information", MessageBoxButtons.OK, MessageBoxIcon.Information)
        End If
    End Sub

Private Sub CalculateCellValueSoilQuality(ByVal CellExtent As StraightRectangle, ByRef CellValue As Double)
        Dim CellULx As Double = CellExtent.UpperLeftPoint.X
        Dim CellULy As Double = CellExtent.UpperLeftPoint.Y
        Dim CellLRx As Double = CellExtent.LowerRightPoint.X
        Dim CellLRy As Double = CellExtent.LowerRightPoint.Y

        Dim xp As Double = CellExtent.Center.X
        Dim yp As Double = CellExtent.Center.Y
        Dim CellCenterPointShape As New PointShape(xp, yp)

        Dim pHValue As Double = Map1.Layers(0).GridLayers(0).GetCellValue(xp, yp)
        Dim MgValue As Double = Map1.Layers(0).GridLayers(1).GetCellValue(xp, yp)

        If pHValue <> -9999 And MgValue <> -9999 Then
            CellValue = GetSoilQualityIndex(pHValue, MgValue)
        Else
            CellValue = -9999
        End If
    End Sub

Private Function GetSoilQualityIndex(ByVal pHValue As Double, ByVal MgValue As Double) As Double
        Dim pHFactor, MgFactor As Integer
      
        If pHValue < 6.4 Then
            pHFactor = 5
        ElseIf pHValue < 6.6 Then
            pHFactor = 6
        ElseIf pHValue < 6.8 Then
            pHFactor = 7
        ElseIf pHValue < 7.0 Then
            pHFactor = 8
        ElseIf pHValue < 7.2 Then
            pHFactor = 9
        ElseIf pHValue < 7.4 Then
            pHFactor = 10
        ElseIf pHValue < 7.6 Then
            pHFactor = 9
        ElseIf pHValue < 7.8 Then
            pHFactor = 8
        Else
            pHFactor = 7
        End If

        If MgValue < 45 Then
            MgFactor = 1
        ElseIf MgValue < 50 Then
            MgFactor = 2
        ElseIf MgValue < 55 Then
            MgFactor = 3
        ElseIf MgValue < 60 Then
            MgFactor = 4
        ElseIf MgValue < 65 Then
            MgFactor = 5
        ElseIf MgValue < 70 Then
            MgFactor = 6
        ElseIf MgValue < 75 Then
            MgFactor = 7
        ElseIf MgValue < 80 Then
            MgFactor = 8
        ElseIf MgValue < 85 Then
            MgFactor = 9
        Else
            MgFactor = 10
        End If

        Return pHFactor * MgFactor

    End Function

Now that we have created the GRID for Soil Quality, we are ready to display it.

Dim STR_GridFile As String = "..\SoilQuality.grd"
        Dim gridlayer As GridLayer = New GridLayer(STR_GridFile)
        gridlayer.LowerThreshold = 0
        gridlayer.UpperThreshold = Double.MaxValue
        gridlayer.ThresholdUnit = ThresholdUnits.miles
        gridlayer.IsNonValueTransparent = True
        Dim gridranges As RangeCollection = New RangeCollection()
        gridranges.Add(New Range(0, 25, Color.FromArgb(216, 253, 193)))
        gridranges.Add(New Range(25, 40, Color.FromArgb(182, 251, 140)))
        gridranges.Add(New Range(40, 55, Color.FromArgb(137, 248, 71)))
        gridranges.Add(New Range(55, 70, Color.FromArgb(96, 239, 10)))
        gridranges.Add(New Range(70, 85, Color.FromArgb(75, 188, 7)))
        gridranges.Add(New Range(85, 100, Color.FromArgb(52, 131, 5)))
        gridlayer.Ranges = gridranges
        Map1.Layers(0).GridLayers.Add(gridlayer)

Fig 6: GRID for soil quality.

Fig 7: Soil quality map legend.

Conclusion

We hope that the simple examples in this article have shown you more clearly how GRID can help you perform and enhance various tasks in your industry. Keep in mind that the examples shown are very simplistic and that GRID is capable for far more sophisticated uses. In the future, we plan to add more white papers related to GRID to show how this technology can be used for a variety of additional tasks, such as surface analysis, diffusion modeling, zonal functions and more types of interpolations.

Downloads

• Download GRID Demo Application
• Note: You will need to download the Map Suite Desktop .NET GIS Control Free Evaluation in order to run the example in this article. (Registration Required)

Have Questions?

Do you have any questions or comments about this article? If so, feel free post your questions on the Map Suite Discussion Forums or leave feedback at the ThinkGeo Blog.

Additional Resources

• Map Suite QuickStart Guides
• Map Suite API Documentation
• Map Suite FAQs

About the Author

Val Guillou is a GIS analyst and developer at ThinkGeo, a software company specializing in geospatial software with an emphasis on software development tools and GPS tracking solutions.