Plotting projected Albers grid with geographic data

Q&A PyGMT Q&A
1

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)

2

Hi.
Could you solve it?

3

That grid is already project in UTM. One cannot plot it in another projection without … reprojecting it first.

4

To clarify the situation, the GMT and PyGMT projection system only knows how to project from geographic (long-lat) to UTM, Albers, or other projection, or the inverse. It does not know how to project from UTM to Albers. If the data is in UTM, then you need to reproject it to geographic before plotting in Albers.

5

This is (well …) actually not true but unfortunately also not that false for some transformations. Tried UTM → Alberts but it failed, possibly because the currentcode could not create a direct GDAL transform from one to the other. This is a part of GMT where coders help would be nice.

For mapproject and grdproject we can go directly from the referencing system A to B without the intermediate step of converting to geographic coordinates. That is obtained (like in cs2cs) by using the +to keyword. Example: -J EPSG:4326+to+proj=aeqd+ellps=WGS84+units=m. A much awaited bonus is also that we now do not need to set -R to do point coordinate conversions.

6

@Joaquim That is great news that this is now possible for some transformations in GMT. I use QGIS for some of my maps and if the source coordinate system and projection is defined, QGIS automatically calls GDAL with the right parameters to make a map in whatever projection you choose.