|
| 1 | +""" |
| 2 | +Calculating grid gradient and radiance |
| 3 | +-------------------------------------- |
| 4 | +The :meth:`pygmt.grdgradient` method calculates the gradient of a grid file. |
| 5 | +In the example shown below we will see how to calculate a hillshade map based |
| 6 | +on a Data Elevation Model (DEM). As input :meth:`pygmt.grdgradient` gets |
| 7 | +a :class:`xarray.DataArray` object or a path string to a grid file, calculates |
| 8 | +the respective gradient and returns it as an :class:`xarray.DataArray` object. |
| 9 | +We will use the ``radiance`` parameter in order to set the illumination source |
| 10 | +direction and altitude. |
| 11 | +""" |
| 12 | + |
| 13 | +import pygmt |
| 14 | + |
| 15 | +# Define region of interest around Yosemite valley |
| 16 | +region = [-119.825, -119.4, 37.6, 37.825] |
| 17 | + |
| 18 | +# Load sample grid (3 arc second global relief) in target area |
| 19 | +grid = pygmt.datasets.load_earth_relief(resolution="03s", region=region) |
| 20 | + |
| 21 | +# calculate the reflection of a light source projecting from west to east |
| 22 | +# (azimuth of 270 degrees) and at a latitude of 30 degrees from the horizon |
| 23 | +dgrid = pygmt.grdgradient(grid=grid, radiance=[270, 30]) |
| 24 | + |
| 25 | +fig = pygmt.Figure() |
| 26 | +# define figure configuration |
| 27 | +pygmt.config(FORMAT_GEO_MAP="ddd.x", MAP_FRAME_TYPE="plain") |
| 28 | + |
| 29 | +# --------------- plotting the original Data Elevation Model ----------- |
| 30 | + |
| 31 | +pygmt.makecpt(cmap="gray", series=[200, 4000, 10]) |
| 32 | +fig.grdimage( |
| 33 | + grid=grid, |
| 34 | + projection="M12c", |
| 35 | + frame=['WSrt+t"Original Data Elevation Model"', "xa0.1", "ya0.1"], |
| 36 | + cmap=True, |
| 37 | +) |
| 38 | + |
| 39 | +fig.colorbar(position="JML+o1.4c/0c+w7c/0.5c", frame=["xa1000f500+lElevation", "y+lm"]) |
| 40 | + |
| 41 | +# --------------- plotting the hillshade map ----------- |
| 42 | + |
| 43 | +# Shift plot origin of the second map by 12.5 cm in x direction |
| 44 | +fig.shift_origin(xshift="12.5c") |
| 45 | + |
| 46 | +pygmt.makecpt(cmap="gray", series=[-1.5, 0.3, 0.01]) |
| 47 | +fig.grdimage( |
| 48 | + grid=dgrid, |
| 49 | + projection="M12c", |
| 50 | + frame=['lSEt+t"Hillshade Map"', "xa0.1", "ya0.1"], |
| 51 | + cmap=True, |
| 52 | +) |
| 53 | + |
| 54 | +fig.show() |
0 commit comments