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:

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.

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:

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:

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

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:

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()

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. 

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

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.

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

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:

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

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)