How to get a dip grid

Q&A PyGMT Q&A
1

With PyGMT, how do I get a grid of maximum dip (in degrees) of a topography grid?
I see how to get azimuth …

azimuth_grid = pygmt.grdgradient(grid=topo_grid, direction=‘a’)

… but I don’t see an option for dip.

Thanks!

2

Hi Trevor

Check the -S option. (I am not sure if it was implemented in pygmt).

3

Thanks. I think -S is implemented with the “slope_file=” option, but I’m still struggling to understand the right syntax and what it actually does … the help pages are still under development I think.

4

Here’s a work-around … write a temporary netCDF file to disk and read it back in.

[1] Remove old version of temporary netCDF file

try:
os.remove(‘tmp_slope.nc’)
except OSError:
pass

[2] Use PyGMT to get a dataarray for azimuth, and also write a temporary netCDF of slope to file

da_az = pygmt.grdgradient(da_topo, direction=‘a’, slope_file=‘tmp_slope.nc’)

[3] Read slope netCDF back into a dataarray

da_slope = xr.open_dataarray(‘tmp_slope.nc’)

[4] Calculate dip in degrees from slope dataarray using Numpy

da_dip = np.arctan(da_slope) * (180/np.pi)

3 Likes