Joseph Barraud
This notebook presents a few examples of how to use interpies.
At the core of interpies is the grid class that is used to load gridded data into memory. A large number of grid formats can be read thanks to rasterio and the underlying GDAL library. Once it has been created, the grid can be transformed, filtered and displayed alongside other grids or objects.
Let's start by importing the interpies module.
In [ ]:
import interpies
In [1]:
# show the plots in the notebook
% matplotlib inline
Thanks to rasterio, one of the required libraries of interpies, raster data can be easily loaded into a grid object, which is basically a numpy 2D array with attributes such as cell size, coordinate system and nodata value attached to it.
As an example, let's read some aeromagnetic data collected by the USGS over the town of Blanca in Colorado (Bankey and Grauch, 2004). The survey lines, the gridded data and the report are available from this page.
I have downloaded the Reduced-to-the-pole (RTP) magnetic data that is provided in GXF format. The grid can be loaded with the function open.
In [4]:
inFile = r'..\data\brtpgrd.gxf'
grid1 = interpies.open(inFile)
grid objects have a few attributes such as cellsize and ncols:
In [5]:
print('The cell size is {} m.'.format(grid1.cellsize))
print('The number of columns is {}.'.format(grid1.ncols))
This kind of metadata can be conveniently listed with the info method:
In [6]:
grid1.info()
In [7]:
grid1.data
Out[7]:
In [8]:
grid1.show()
Out[8]:
You may have noticed that the show function returns matplotlib axes. This is useful to keep a handle on the figure and display additional objects.
The show function has a lot of options and parameters. As you can see, it displays the data by default with a hillshading effect, the colormap is the classic clra color table by Geosoft, and histogram equalization is applied to rebalance the distribution of colours that would otherwise be dominated by extreme negative and positive anomalies.
It is possible to remove the hillshade, to use different colormaps and to apply linear normalisation in this way:
In [10]:
ax = grid1.show(hs=False, cmap='coolwarm', cmap_norm='none', title='a grid in coolwarm and no hillshade')
Zooming on a specific area can be done by clipping. Let's first look at the extent of the grid - these are the outer bounds of the grid, i.e. they include half a cell size in all directions. Additionally, the convention here is that the geographic location of the grid is defined by the centre of the cells.
In [11]:
# west, south, east, north
grid1.extent
Out[11]:
The clip method returns a new grid.
In [13]:
grid2 = grid1.clip(xmin=450000, xmax=460000, ymin=4135000, ymax=4150000)
grid2.show();
In [14]:
grid1.clip(xmin=456000, xmax=466000, ymin=4120000, ymax=4130000).show();
In [16]:
# Horizontal derivative in the x (west-east) direction
# The vertical exaggeration (zf) has to be increased for the hillshade to look good.
ax = grid1.dx().show(zf=1000)
In [17]:
# tilt angle
ax = grid1.tilt().show()
To conclude this first notebook, let's calculate the horizontal gradient magnitude of the pseudogravity. This particular transform is described in the report of Bankey and Grauch (2004). The pseudogravity is basically a vertical integral of the RTP magnetic anomalies. The horizontal gradient highlights the steepest slopes.
The calculation can be done in one line using a chain of functions. I am starting with a linear detrend to improve the quality of the vertical integration.
In [19]:
ax = grid1.detrend().vi().hgm().show(title='Horizontal Gradient Magnitude of the Pseudogravity')