I see the figure 28 example on GMT with UTM coordinates. This makes sense because the method for way the UTM projection is specified in GMT. I can’t seem to get it to work with Albers projections because it assumes lat/long.
I’ve tried supplying the filename as the region, manually specifying the projected coordinates and converting the region coordinates to degrees (which runs but the plot doesn’t line up because the coastal plot is being plotted on the geographic axis) but I still get a map projection error. Is it even possible to plot a projected albers grid with geographic units in GMT? Here’s what I’ve tried.
import pygmt
import subprocess
def figure28():
subprocess.call(["gmt", "makecpt", "-Ccopper", "-T0/1500"])
fig = pygmt.Figure()
filename = "@Kilauea.utm.nc"
fig.grdimage(filename, projection="x1:160000", shading="+d")
fig.coast(region=filename, projection="u5Q/1:160000", water="lightblue", frame=["5mg5m", "NEsw"])
fig.savefig("figure28.pdf")
def fails():
dem = "dem.dem"
reg = "1535157.3089179904200137,1641757.3089179904200137,1810605.1008632136508822,1937405.1008632136508822"
fig = pygmt.Figure()
fig.grdimage(grid=dem, projection="x1:290000")
fig.coast(region=reg + "+umeter", projection="b-96/23/29.5/45.5/1:290000", rivers="DODGERBLUE2")
fig.savefig("test.pdf")
if __name__ == "__main__":
figure28()
fails()
I attached an example with a python script and a test dem.
example.zip (831.8 KB)