Topography in UTM [m] projection

[mkononets]

@EJFielding

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.

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

@EJFielding

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.

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)