What is a grid?

A grid is a form of spatial data that represents information (such as a gravity intensity, a magnetic reading, or a colour) at points organized as a 2-dimensional array on a right-handed cartesian plane:

Grid data on a plane

The plane on which the grid is located can be oriented in three dimensions relative to a Coordinate System on the Earth. The most common grids are located on a horizontal surface relative to the Coordinate System. For example, common surfaces might be sea-level, or the ground surface, or a constant elevation above or below the ground surface, or a constant elevation. Vertical cross-sections through the Earth are oriented to be orthogonal to the surface of the Earth.

A common term used with grids is the concept of a 'grid cell'. In Geosoft's usage, grids are an array of points at a point location, and a 'grid cell' is the rectangular area that extends half-way to the neighboring grid points.

Spatial reference angles throughout the geosoft.gxpy module will consistently use an angle in degrees azimuth, which is a clockwise-positive angle relative to a coordinate system frame (North, or positive Y). The geosoft.gxapi.GXIMG class specifies the rotation angle in degrees counterclockwise-positive, and other references within the geosoft.gxapi may differ.

Grids are stored as files on the file system, and there are many common grid file formats in existence. When working with grid files in GX Developer you define the grid file format with the use of a decorator string appended to the grid file name. Geosoft supports the 11 formats (and their many derivatives) described in the Grid File Name Decorations section of the GX Developer documentation. For example:

grid file string	grid file type
`'c:/project/mag.grd(GRD)'`	Geosoft format grid.
`'c:/project/mag.tif(TIF)'`	GeoTIF
`'c:/project/image.jpg(IMG;T=5)'`	jpeg image file
`'c:/project/mag.grd(GRD;TYPE=COLOR)'`	Geosoft colour grid

Convert a grid from one format to another

Problem: You have a grid in a Geosoft-supported format, and you need the grid in some other format to use in a different application.

Grid: elevation_surfer.grd, which is a Surfer v7 format grid file.

Approach:

Open the surfer grid with decoration '(SRF;VER=V7)', as documented in Grid File Name Decorations.
Use the gxgrid.Grid.copy class method to create an ER Mapper grid, which will have decoration '(ERM)'.

grid_convert_format.py

import geosoft.gxpy.gx as gx
import geosoft.gxpy.grid as gxgrid

# create context
gxc = gx.GXpy()

# open grid
grid_surfer = gxgrid.Grid.open('elevation_surfer.grd(SRF;VER=V7)')

# copy the grid to an ER Mapper format grid file
grid_erm = gxgrid.Grid.copy(grid_surfer, 'elevation.ers(ERM)', overwrite=True)

exit()

If you run this script with a debugger, break on the exit() (line 13) and inspect the objects. You will see the two grid instances, and you can open these to see the content and properties:

On the file system you will find new ER Mapper format grid files elevation and elevation.ers. You will also see files elevation.ers.gi and elevation.ers.xml, which are Geosoft-only files. The .gi file stores relevant grid information that is not supported natively by the grid file format, and the .xml file stores grid metadata.

Define the coordinate system

In Geosoft all spatial data should have a defined coordinate system which allows data to be located on the Earth. This also takes advantage of Geosoft's ability to reproject data as required. However, in this example the Surfer grid does not store the coordinate system information, but we know that the grid uses projection 'UTM zone 54S' on datum 'GDA94'. Let's modify this script to set the coordinate system, which will be saved as part of the ER Mapper grid, which does have the ability to store the coordinate system description.

In Geosoft, well-known coordinate systems like this can be described using the form 'GDA94 / UTM zone 54S', which conforms to the SEG Grid Exchange Format standard for describing coordinate systems. You only need to set the coordinate_system property of the grid_surfer instance:

grid_convert_format_with_coordinate_system.py

import geosoft.gxpy.gx as gx
import geosoft.gxpy.grid as gxgrid

# create context
gxc = gx.GXpy()

# open grid
grid_surfer = gxgrid.Grid.open('elevation_surfer.grd(SRF;VER=V7)', mode=gxgrid.FILE_READWRITE)

# define the coordinate system
grid_surfer.coordinate_system = 'GDA94 / UTM zone 54S'

# copy the grid to an ER Mapper format grid file
grid_erm = gxgrid.Grid.copy(grid_surfer, 'elevation.ers(ERM)', overwrite=True)

print(str(grid_erm.coordinate_system)) # prints 'GDA94 / UTM zone 54S'

exit()

Note

In this example we have also opened the Surfer grid using mode FILE_READWRITE. This will cause the coordinate system to be saved with the Surfer grid so that the next time the Surfer grid file is opened by Geosoft the coordinate system will already be defined. Geosoft stores the coordinate system information together with other Geosoft specific information in the companion .gi file. This demonstrates how Geosoft uses this extra file to extend the capabilities of limited file formats.

Run this script in a debugger and break on line 16. Inspect the coordinate system property of the grid_erm instance and you will see it contains complete coordinate system details, which have been pulled from Geosoft coordinate system engine. There is another exercise that provides more information on working with coordinate systems in Geosoft:

Calculate basic statistics about data contained in a grid

In this exercise we will work with the data stored in a grid. One common need is to determine some basic statistical information about the grid data, such as the minimum, maximum, mean and standard deviation. This exercise will work with the grid data a number of ways that demonstrate some useful patterns.

Using numpy

The smallest code and most efficient approach is to read the grid into a numpy array and then use the optimized numpy methods to determine statistics. This has the benefit of speed and simplicity at the expense memory, which may be a concern for very large grids, though on modern 64-bit computers with most grids this would be the approach of choice.

grid_statistics_numpy.py

import numpy as np
import geosoft.gxpy.gx as gx
import geosoft.gxpy.grid as gxgrid

# create context
gxc = gx.GXpy()

# open the grid
with gxgrid.Grid.open('elevation_surfer.grd(SRF;VER=V7)') as grid:

    # get the data in a numpy array
    data_values = grid.xyzv()[:, :, 3]

# print statistical properties
print('minimum: ', np.nanmin(data_values))
print('maximum: ', np.nanmax(data_values))
print('mean: ', np.nanmean(data_values))
print('standard deviation:', np.nanstd(data_values))

# prints:
# minimum:  88.061668396
# maximum:  339.345184326
# mean:  167.952990336
# standard deviation: 61.8165048181

line 9	This shows the best way to open a grid file, using the Python with... construct. This ensures that the grid file is closed and resources are freed upon exit from the with...: block. Most geosoft.gxpy classes support this use, and this is supported for all classes that depend on a file on the file system.
line 12	the Grid method xyzv() returns the grid data in a numpy array with shape (nx, ny, 4), where the 4 items for each grid point are the (x, y, z, grid_value). For the statistics we only need the grid_value, and we slice the numpy array to this item only. Grid Null Values Each point in a grid will either have a data value, or will be null. A null is normally stored within a grid file using a special value called the dummy value. Grid methods will handle null values in reasonable ways, converting to Geosoft type-equivalent dummy values within Geosoft classes, to None in a method that returns a grid value, or to numpy.nan when working with numpy float data.
line 16-19	We use the nan-aware numpy statistical methods to calculate the various statistics.

Using Geosoft VVs

Many Geosoft methods will work with a geosoft.gxpy.vv.GXvv, which wraps the geosoft.gxapi.GXVV class that deals with very long single-value vectors. The Geosoft GXVV methods works with Geosoft data types and, like numpy, is optimized to take advantage of multi-core processors to improve performance. The pattern in this exercise reads a grid one grid row at a time, returning a GXvv instance and accumulate statistics in an instance of the geosoft.gxapi.GXST class.

grid_statistics_gxvv.py

import geosoft.gxapi as gxapi
import geosoft.gxpy.gx as gx
import geosoft.gxpy.grid as gxgrid

# create context
gxc = gx.GXpy()

# create a gxapi.GXST instance to accumulate statistics
stats = gxapi.GXST.create()

# open the grid
with gxgrid.Grid.open('elevation_surfer.grd(SRF;VER=V7)') as grid:

    # add data from each row to the stats instance
    for row in range(grid.ny):
        stats.data_vv(grid.read_row(row).gxvv)

# print statistical properties
print('minimum: ', stats.get_info(gxapi.ST_MIN))
print('maximum: ', stats.get_info(gxapi.ST_MAX))
print('mean: ', stats.get_info(gxapi.ST_MEAN))
print('standard deviation:', stats.get_info(gxapi.ST_STDDEV))

# prints:
# minimum:  88.0616683959961
# maximum:  339.3451843261719
# mean:  167.95299033603075
# standard deviation: 61.818035062225405

line 1	Include the geosoft.gxapi module to access the GXST class at line 12.
line 14	We loop through all grid rows.
line 15	The read_row() method returns a geosoft.gxpy.GXvv instance. The gxvv property of a GXvv instance the gxapi.GXVV handle that can be passed into gxapi methods, which in this case is passed to the gxapi.GXST.data_vv() method. All geosoft.gxpy classes that wrap geosoft.gxapi classes will, by convention, have a property named after the gxapi class name but lower-case. In this example the geosoft.gxpy.GXvv instance has property gxvv that is the geosoft.gxapi.GXVV instance handle.

Grid Iterator

A grid instance also behaves as an iterator that works through the grid data points by row, then by column, each iteration returning the (x, y, z, grid_value). In this example we will iterate through all points in the grid and accumulate the statistics a point at a time. This is the least-efficient way to work through a grid, but the pattern can be useful to deal with a very simple need. For example, any Geosoft supported grid can be easily converted to an ASCII file that has lists the x, y, z, grid_value for all points in a grid, which we show at the end of this section.

Here is the code that calculates statistics using the grid iterator:

grid_statistics_iterator.py

import geosoft.gxapi as gxapi
import geosoft.gxpy.gx as gx
import geosoft.gxpy.grid as gxgrid
import geosoft.gxpy.utility as gxu

# this example requires version 9.2.1, which adds iteration support
gxu.check_version('9.2.1')

# create context
gxc = gx.GXpy()

# create a gxapi.GXST instance to accumulate statistics
stats = gxapi.GXST.create()

# add each data to stats point-by-point (slow, better to use numpy or vector approach)
number_of_dummies = 0
with gxgrid.Grid.open('elevation_surfer.grd(SRF;VER=V7)') as grid:
    for x, y, z, v in grid:
        if v is None:
            number_of_dummies += 1
        else:
            stats.data(v)
	total_points = grid.nx * grid.ny

# print statistical properties
print('minimum: ', stats.get_info(gxapi.ST_MIN))
print('maximum: ', stats.get_info(gxapi.ST_MAX))
print('mean: ', stats.get_info(gxapi.ST_MEAN))
print('standard deviation:', stats.get_info(gxapi.ST_STDDEV))
print('number of dummies: ', number_of_dummies)
print('number of valid data points: ', total_points - number_of_dummies)

# prints:
# minimum:  88.0616683959961
# maximum:  339.3451843261719
# mean:  167.95299033603075
# standard deviation: 61.818035062225405
# number of dummies:  4875
# number of valid data points:  20199

line 7	Check that the minimum version required for this code to work. This example depends on the Grid class being iterable, which was introduced in version 9.2.1.
line 18	Iterate through all values in the grid, adding non-dummy items to the stats instance.

As another very simple example, here we use the grid iterator pattern to print the coordinate system followed by the (x, y, z, grid_value) for every non-dummy value in a grid:

import geosoft.gxpy as gxpy
gxc = gxpy.gx.GXpy()
with gxpy.grid.Grid.open('elevation_surfer.grd(SRF;VER=V7)') as grid:
    print('coordinate_system: ', grid.coordinate_system)
    for x, y, z, v in grid:
        if v is not None:
            print(x, y, z, v)