# Plotting projected Albers grid with geographic data

**URL:** <https://forum.generic-mapping-tools.org/t/plotting-projected-albers-grid-with-geographic-data/2408>\
**Category:** PyGMT Q&A\
**Created:** [December 3, 2021, 5:32am UTC](https://forum.generic-mapping-tools.org/t/plotting-projected-albers-grid-with-geographic-data/2408 "2021-12-03T05:32:00Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![atrevisan21](https://yyz1.discourse-cdn.com/flex047/user_avatar/forum.generic-mapping-tools.org/atrevisan21/32/1857_2.png) [@atrevisan21](https://forum.generic-mapping-tools.org/u/atrevisan21)\
**Post date:** [December 3, 2021, 5:32am UTC](https://forum.generic-mapping-tools.org/t/plotting-projected-albers-grid-with-geographic-data/2408/1 "2021-12-03T05:32:00Z")

</div>

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.

```python
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](https://forum.generic-mapping-tools.org/uploads/short-url/Aoi9HVNTejO9d0E9BNrcTppQ8yP.zip) (831.8 KB)

---

<div class="post-metadata">

**Author:** ![KPinedaO](https://yyz1.discourse-cdn.com/flex047/user_avatar/forum.generic-mapping-tools.org/kpinedao/32/2402_2.png) [@KPinedaO](https://forum.generic-mapping-tools.org/u/KPinedaO)\
**Post date:** [January 6, 2023, 7:53pm UTC](https://forum.generic-mapping-tools.org/t/plotting-projected-albers-grid-with-geographic-data/2408/2 "2023-01-06T19:53:44Z")

</div>

Hi.  
Could you solve it?

---

<div class="post-metadata">

**Author:** ![Joaquim](https://yyz1.discourse-cdn.com/flex047/user_avatar/forum.generic-mapping-tools.org/joaquim/32/16_2.png) [@Joaquim](https://forum.generic-mapping-tools.org/u/Joaquim)\
**Post date:** [January 6, 2023, 8:32pm UTC](https://forum.generic-mapping-tools.org/t/plotting-projected-albers-grid-with-geographic-data/2408/3 "2023-01-06T20:32:30Z")

</div>

> [@atrevisan21](#):
>
> `@Kilauea.utm.nc`

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

---

<div class="post-metadata">

**Author:** ![EJFielding](https://yyz1.discourse-cdn.com/flex047/user_avatar/forum.generic-mapping-tools.org/ejfielding/32/82_2.png) [@EJFielding](https://forum.generic-mapping-tools.org/u/EJFielding)\
**Post date:** [January 15, 2023, 4:13am UTC](https://forum.generic-mapping-tools.org/t/plotting-projected-albers-grid-with-geographic-data/2408/4 "2023-01-15T04:13:24Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![Joaquim](https://yyz1.discourse-cdn.com/flex047/user_avatar/forum.generic-mapping-tools.org/joaquim/32/16_2.png) [@Joaquim](https://forum.generic-mapping-tools.org/u/Joaquim)\
**Post date:** [January 16, 2023, 4:18pm UTC](https://forum.generic-mapping-tools.org/t/plotting-projected-albers-grid-with-geographic-data/2408/5 "2023-01-16T16:18:26Z")

</div>

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](https://docs.generic-mapping-tools.org/latest/mapproject.html) and [grdproject](https://docs.generic-mapping-tools.org/latest/grdproject.html) 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.

---

<div class="post-metadata">

**Author:** ![EJFielding](https://yyz1.discourse-cdn.com/flex047/user_avatar/forum.generic-mapping-tools.org/ejfielding/32/82_2.png) [@EJFielding](https://forum.generic-mapping-tools.org/u/EJFielding)\
**Post date:** [January 16, 2023, 6:51pm UTC](https://forum.generic-mapping-tools.org/t/plotting-projected-albers-grid-with-geographic-data/2408/6 "2023-01-16T18:51:28Z")

</div>

@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.