Topography in UTM [m] projection

earthistics · October 28, 2025, 11:44am

Hello,
I’m trying to plot a topographic map using Transverse Mercator projection (with thicks like 420000 → 470000 Easting and 4510000->4535000 North). I’m giving a try with this:

import pygmt
import pyproj
import xarray as xr

region_utm = [420000, 470000, 4510000, 4535000]  # [xmin, xmax, ymin, ymax]
utm_proj = pyproj.Proj(proj="utm", zone=33, ellps="WGS84", north=True)

lon_min, lat_min = utm_proj(region_utm[0], region_utm[2], inverse=True)
lon_max, lat_max = utm_proj(region_utm[1], region_utm[3], inverse=True)
region_lonlat = [lon_min, lon_max, lat_min, lat_max]

grid_geo = pygmt.datasets.load_earth_relief(resolution="03s", region=region_lonlat)

grid_utm = pygmt.grdproject(grid=grid_geo, projection="EPSG:32633")

grid_utm = pygmt.grdcut(grid=grid_utm, region=region_utm)

fig = pygmt.Figure()

fig.grdimage(
    grid=grid_utm,
    region=region_utm,
    projection="X15c/10c",
    shading=True,
    cmap="gray",
    frame=["xaf", "yaf"]
)

# coastlines (same UTM region)
fig.coast(
    region=region_utm,
    shorelines=True,
    water="lightblue",
    resolution="f",
    frame=["xaf", "yaf"]
)

fig.show()

However, it is not working. I get this:

pygmt.exceptions.GMTCLibError: Module 'coast' failed with status code 74:
coast [ERROR]: Map region exceeds 360 degrees
coast [ERROR]: General map projection error

So, the coast module does not work in this case (if I comment the coast module lines, the plot works!).
I would appreciate your support with this issue.

Thanks

EJFielding · October 28, 2025, 1:14pm

You projected your topographic data into UTM coordinates, and then you set up a linear projection for your map in UTM coordinates. The coastline database is still in lat-long, so you would need to tell PyGMT how to project the coastline into the UTM projection for your area by using the geographic coordinates of the corners that you calculated as the region and the UTM projection as the projection.
This is an inconvenience in GMT and PyGMT that you have to keep track yourself of what data is in what projection and provide the correct projection in the plotting statements when mixing data from different projections.

mkononets · October 28, 2025, 4:57pm

yeah right and guess the map scale for the user and what not.

The topic starter has already done all the corner coords calculations, using his pyproj.Proj calls. The two things left was 1) to specify the region correctly so coast recognizes it’s corners, and 2) use scale, not size (height/width) for both Cartesian and geographical projections.

The corrected script is below. Note region and projection syntax.

import pygmt
import pyproj
import xarray as xr

region_utm = [420000, 470000, 4510000, 4535000]  # [xmin, xmax, ymin, ymax]
utm_proj = pyproj.Proj(proj="utm", zone=33, ellps="WGS84", north=True)

lon_min, lat_min = utm_proj(region_utm[0], region_utm[2], inverse=True)
lon_max, lat_max = utm_proj(region_utm[1], region_utm[3], inverse=True)
region_lonlat = [lon_min, lon_max, lat_min, lat_max]

grid_geo = pygmt.datasets.load_earth_relief(resolution="03s", region=region_lonlat)

grid_utm = pygmt.grdproject(grid=grid_geo, projection="EPSG:32633")

grid_utm = pygmt.grdcut(grid=grid_utm, region=region_utm)

fig = pygmt.Figure()

scale="1:500000"
# note one must use the same scale for Cartesian and UTM projection specification
# one cannot use X15c/10c with arbitrary map width and height and hope
# the following coast call resolves it to lon/lat
# NB map frame specification: cartesian tics bottom (South) and left (West) axes
fig.grdimage(
    grid=grid_utm,
    region=region_utm,
    projection=f"x{scale}",
    shading=True,
    cmap="gray",
    frame=["SW", "xaf", "yaf"]
)

# coastlines (same UTM region)
# notes
# 1) region synthax in form of {lon_min}/{lat_min}/{lon_max}/{lat_max}+r
# +r tells its low left and upper right corner coords, needed
# to produce a square map frame matching utm
# the ...+r syntax is needed as lon/lat grid is not rectilinear in UTM coords
# 2) one must use the same scale for Cartesian and UTM projection specification
# I don't think pygmt documents this, neither gmt does anymore.
# the (in)famous GMT Exapmle 28
# https://docs.generic-mapping-tools.org/latest/gallery/ex28.html#example-28
# uses (not so) modern syntactic sugar to retrieve region by specifying
# the grid file. Here this is not possible as there's no grid file.
#
# NB map frame specification: geographical tics top (North) and right (East) axes
fig.coast(
    region=f"{lon_min}/{lat_min}/{lon_max}/{lat_max}+r",
    projection=f"EPSG:32633/{scale}",
    shorelines=True,
    water="lightblue",
    resolution="f",
    frame=["NE", "xaf", "yaf"]
)

fig.savefig("Topography_in_UTM.png", show=True)

earthistics · October 30, 2025, 7:04pm

Thanks a lot @mkononets and @EJFielding for replying! I appreciate your solution, @mkononets ! However, it seems that pyGMT is very “reluctant” when working with UTM projection. I’ve being struggling with the map scale and the very last option I’ve found is a manual solution, but I’m not convinced. Before adopting a manual solution, I’ve been trying with these three options:

##Option 1
scale_x = region_utm[0] 
scale_y = region_utm[2] 
fig.basemap(map_scale=f"j{scale_x}/{scale_y}+w5k+u") 

##Option 2
fig.basemap(map_scale="jBL+w5k+lkm")

##Option 3
utm_e_center = 0.5 * (region_utm[0] + region_utm[1])
utm_n_center = 0.5 * (region_utm[2] + region_utm[3])
lon_center, lat_center = utm_proj(utm_e_center, utm_n_center, inverse=True)
fig.basemap(
    region=region_utm,
    projection=f"x{scale}",
    frame=True,
    map_scale=f"jBR+c{lon_center}/{lat_center}+w5k+u+l5\\ km+o0.2c"

None of these options have worked. Do you maybe have any idea to fix this?

Thanks!

mkononets · October 30, 2025, 7:14pm

what do you want to do? explain and present some complete code that can be executed and the resulting plot.

earthistics · October 30, 2025, 7:43pm

Continuing with the previous code (the one you improved), from the UTM map, I’m trying to add the map scale bar and tested the three different options I mentioned in my previous post (in the code option1, option2, option3). However, none of these options worked: the scale bar doesn’t appear.

import pygmt
import pyproj
import xarray as xr

region_utm = [420000, 470000, 4510000, 4535000]  # [xmin, xmax, ymin, ymax]
utm_proj = pyproj.Proj(proj="utm", zone=33, ellps="WGS84", north=True)

lon_min, lat_min = utm_proj(region_utm[0], region_utm[2], inverse=True)
lon_max, lat_max = utm_proj(region_utm[1], region_utm[3], inverse=True)
region_lonlat = [lon_min, lon_max, lat_min, lat_max]

grid_geo = pygmt.datasets.load_earth_relief(resolution="03s", region=region_lonlat)
grid_utm = pygmt.grdproject(grid=grid_geo, projection="EPSG:32633")
grid_utm = pygmt.grdcut(grid=grid_utm, region=region_utm)

fig = pygmt.Figure()
scale="1:500000"

fig.grdimage(
    grid=grid_utm,
    region=region_utm,
    projection=f"x{scale}",
    shading=True,
    cmap="gray",
    frame=["SW", "xaf", "yaf"]
)

#Map-scale bar
##Option 1
scale_x = region_utm[0] 
scale_y = region_utm[2] 
fig.basemap(map_scale=f"x{scale_x}/{scale_y}+w5k+u") 

##Option 2
fig.basemap(map_scale="jBL+w5k+lkm")

##Option 3
utm_e_center = 0.5 * (region_utm[0] + region_utm[1])
utm_n_center = 0.5 * (region_utm[2] + region_utm[3])
lon_center, lat_center = utm_proj(utm_e_center, utm_n_center, inverse=True)
fig.basemap(
    region=region_utm,
    projection=f"x{scale}",
    frame=True,
    map_scale=f"jBR+c{lon_center}/{lat_center}+w5k+u+l5\\ km+o0.2c")

fig.coast(
    region=f"{lon_min}/{lat_min}/{lon_max}/{lat_max}+r",
    projection=f"EPSG:32633/{scale}",
    shorelines=True,
    water="lightblue",
    resolution="f",
    frame=["NE", "xaf", "yaf"]
)

fig.savefig("Topography_in_UTM.png", show=True)

mkononets · October 30, 2025, 7:51pm

now run your code and show the plot here, and explain what you want to achieve and how it is different from what you get.

earthistics · October 30, 2025, 7:58pm

The scale bar doesn’t appear. Only the label (bottom-right corner).

But something like this is expected:

mkononets · October 30, 2025, 9:22pm

First, you got two axes in projected meters. That is 420000 is exactly 10 km away from 430000 on x axis. You don’t really need a scale, right?

Second, you’ve apparently looked at Example 28. Why didn’t you just copy their map scale specification? You need to call basemap to plot your map scale after coast to retain geographical projection and region otherwise it’ll get water plotted over it:

...
fig.coast(
...
)

lon, lat = utm_proj(region_utm[0], region_utm[2], inverse=True)
fig.basemap(map_scale=f'jBL+c{lat}+f+w10k+l{scale}+u+o0.5c')

fig.savefig("Topography_in_UTM.png", show=True)

UPD there’s no need to separately call basemap as you can apparently specify map scale right in the coast call:

lon, lat = utm_proj(region_utm[0], region_utm[2], inverse=True)
fig.coast(
    region=f'{lon_min}/{lat_min}/{lon_max}/{lat_max}+r',
    projection=f'EPSG:32633/{scale}',
    shorelines=True,
    water='lightblue',
    resolution='f',
    frame=['NElb', 'xaf', 'yaf'],
    map_scale=f'jBL+c{lat}+f+w10k+l{scale}+u+o0.5c'
)

UPD2 Example 28 used a separate basemap call to plot scale bar because they wanted different FONT_ANNOT_PRIMARY for their coast and basemap calls. Apparently FONT_ANNOT_PRIMARY is used both for map frame and for scale bar tick annotations, which may be not desirable.