# Using oblique mercator with earth relief bathymetry

**URL:** <https://forum.generic-mapping-tools.org/t/using-oblique-mercator-with-earth-relief-bathymetry/5474>\
**Category:** PyGMT Q&A\
**Created:** [October 21, 2024, 4:15pm UTC](https://forum.generic-mapping-tools.org/t/using-oblique-mercator-with-earth-relief-bathymetry/5474 "2024-10-21T16:15:34Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![ledeczi](https://avatars.discourse-cdn.com/v4/letter/l/df788c/32.png) [@ledeczi](https://forum.generic-mapping-tools.org/u/ledeczi)\
**Post date:** [October 21, 2024, 4:15pm UTC](https://forum.generic-mapping-tools.org/t/using-oblique-mercator-with-earth-relief-bathymetry/5474/1 "2024-10-21T16:15:34Z")

</div>

I am creating maps in oblique mercator projections for a publication.  
I have succeeded in creating these maps with the following code:

```auto
fig = pygmt.Figure()
fig.coast(projection="OC-3/146.0/145.7/-10.7/6c",
    region="143/-7/160/0+r",
    frame="afg",
    land="olivedrab",
    shorelines="1/thin",
   water="aliceblue",
)
fig.show(width=1000, dpi=500)

```

However, for some projections (like the one above) the coastlines look very jagged. Is this just a result of the projection? If I try the same with a different projection (OC-3/146.0/143.8/36.5/6c) the coastlines look more normal. I’m hoping to make these panels in the same figure and I’m wondering why they look so different.

 ![Screenshot 2024-10-21 at 9.10.03 AM](https://canada1.discourse-cdn.com/flex047/uploads/gmt/original/2X/9/9ef63d765055f4546f1e3e6781a5c93edb565b78.jpeg)

Ultimately, I need to add bathymetry to the background. I am trying to do this by:

```auto
grid = pygmt.datasets.load_earth_relief(resolution="30s", region=[143.0, 154.0, -6.0, -2.0])
fig = pygmt.Figure()
fig.grdimage(grid=grid, projection="OC-3/146.0/145.7/-10.7/6c", frame="a", cmap="geo")
fig.colorbar(frame=["a1000", "x+lElevation", "y+lm"])
fig.show()

```

I am getting the following error:  
Module ‘psconvert’ failed with status code 79:  
psconvert [ERROR]: System call [gs -q -dNOSAFER -dNOPAUSE -dBATCH -sDEVICE=bbox -DPSL\_no\_pagefill -dMaxBitmap=2147483647 -dUseFastColor=true ‘/Users/annaledeczi/.gmt/sessions/gmt\_session.29035/gmt\_27.ps-’ 2\> ‘/Users/annaledeczi/.gmt/sessions/gmt\_session.29035/psconvert\_29035c.bb’] returned error 256.

Is it possible to use an oblique mercator proejction on earth relief imagery?  
Thanks!

---

<div class="post-metadata">

**Author:** ![mkononets](https://yyz1.discourse-cdn.com/flex047/user_avatar/forum.generic-mapping-tools.org/mkononets/32/4029_2.png) [@mkononets](https://forum.generic-mapping-tools.org/u/mkononets)\
**Post date:** [October 22, 2024, 12:28am UTC](https://forum.generic-mapping-tools.org/t/using-oblique-mercator-with-earth-relief-bathymetry/5474/2 "2024-10-22T00:28:44Z")

</div>

> [@ledeczi](#):
>
> for some projections (like the one above) the coastlines look very jagged. Is this just a result of the projection? If I try the same with a different projection (OC-3/146.0/143.8/36.5/6c) the coastlines look more normal. I’m hoping to make these panels in the same figure and I’m wondering why they look so different.

the bottom part of your region size is actually very small, so you need to specify the use of high-resolution coastline, see [Coastlines and borders — PyGMT](https://www.pygmt.org/latest/tutorials/basics/coastlines.html#resolutions)  
your “different projection” actually defines a much larger region so the default coastline resolution looks “more normal”.

> [@ledeczi](#):
>
> Is it possible to use an oblique mercator proejction on earth relief imagery?

yes it is possible.  
Thing is you have to extract a rectangular earth relief dataset that covers your whole oblique region of `region="143/-7/160/0+r"`

in command-line interface the rectangular region can be autodetected using `mapproject`’s `-Wr` or `-WR` option, see this [illustration to `mapproject`’s `-W` option](https://docs.generic-mapping-tools.org/latest/mapproject.html#id1)

as far as I can understand, `mapproject` is not directly available in `pyGMT`. You can call it using pygmt.clib.Session(), like below (or find out the extent of the rectangular region extent visually by plotting grid with a small step size):

```auto
import pygmt
from pygmt.clib import Session
from pygmt.helpers import GMTTempFile

projection="OC-3/146.0/145.7/-10.7/6c"
region = "143/-7/160/0+r"

with Session() as s:
    with GMTTempFile() as fout:
        s.call_module('mapproject', args=f'-R{region} -J{projection} -Wr > {fout.name}')
        R = fout.read().strip()

W,E,S,N = R.split()
print(f'The rectangular region is [{W}, {E}, {S}, {N}]')

grid = pygmt.datasets.load_earth_relief(resolution="30s", region=[W, E, S, N])
fig = pygmt.Figure()
fig.grdimage(grid=grid, projection=projection, region=region, frame="afg1", cmap="geo")
fig.colorbar(frame=["a1000", "x+lElevation", "y+lm"])
fig.savefig("obr.png", show=True)

```

NB `frame="afg1"` shows gridlines plotted every 1 degree, visually illustrating the extent of the rectangular region that needs to be extracted using `pygmt.datasets.load_earth_relief()`, approx` [135, 160, -8, 8]`

---

<div class="post-metadata">

**Author:** ![ledeczi](https://avatars.discourse-cdn.com/v4/letter/l/df788c/32.png) [@ledeczi](https://forum.generic-mapping-tools.org/u/ledeczi)\
**Post date:** [October 22, 2024, 5:34pm UTC](https://forum.generic-mapping-tools.org/t/using-oblique-mercator-with-earth-relief-bathymetry/5474/3 "2024-10-22T17:34:06Z")

</div>

Thank you so much! I’ve been puzzling over this for weeks.