This article will introduce the Map Suite developer community to the GRID format and help demonstrate how useful GRID can be. 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.
Introduction
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.
| pH Level |
Factor |
| 6.20 to 6.40 |
5 |
| 6.40 to 6.60 |
6 |
| 6.60 to 6.80 |
7 |
| 6.80 to 7.00 |
8 |
| 7.00 to 7.20 |
9 |
| 7.20 to 7.40 |
10 |
| 7.40 to 7.60 |
9 |
| 7.60 to 7.80 |
8 |
| 7.80 to 8.00 |
7 |
Table 1: pH factors
| Mg Level |
Factor |
| 40 to 45 |
1 |
| 45 to 50 |
2 |
| 50 to 55 |
3 |
| 55 to 60 |
4 |
| 60 to 65 |
5 |
| 65 to 70 |
6 |
| 70 to 75 |
7 |
| 75 to 80 |
8 |
| 80 to 85 |
9 |
| 85 to 90 |
10 |
Table 2: Mg factors
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
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
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.