Mixing UTM and geographic data sets

Q&A PyGMT Q&A
1

Hi, everyone. I hope someone can help me with this issue.
Last few days. I´ve been struggling with this issue.

I want to plot a UTM elevation grid and plot above data in geographic coordinates.
I saw this information (https://docs.generic-mapping-tools.org/dev/gallery/ex28.html); however, I´m not able to make it work.

So far, I can plot the UTM grid, but I want to zoom to a specific region and I´m not sure how to do that. This is my code:

import pygmt
fig = pygmt.Figure()
filename = "grid_UTM.tif"
fig.grdimage(filename, projection="x1:160000", shading="+d")
fig.basemap(region=filename, projection="u16P/1:160000", frame=["WSne", "af"])
fig.show()

So in basemap, should I add the region coordinates in lat, lon or in UTM format?
Thanks for the help.

2

I have been combining UTM data and geographic coordinates data for many years in regular GMT, now called classic mode. I always include the region and projection information in every GMT call, and I have a separate “projUTM” and “projGeog”, with corresponding regions. I am not sure exactly how PyGMT handles the projection and region information, but I would guess you can do the same thing in PyGMT by adding the relevant UTM or geographic projection and region to each call.

3

One other thing about combining the UTM and geographic projections is that you need to ensure that the regions for the two projections are the same. In classic GMT, you can use the mapproject program to project corner points from geographic to UTM. I have lines like this in many of my old GMT scripts:

mapproject Acapulco-corn-dd.xy -R-102/-98/10/20 -Ju14/1:1 -F -C > Acapulco-corn-UTM14.xy

In Python, you can do this much more easily than my old-school method of hand-editing the Acapulco-corn-dd.xy input file and getting the numbers from the output file.

4

If you want to display geographic coordinates
You can use grdproject to transform grid_UTM.tif from UTM coordinates to geographic coordinates. Usually, you’d plot the basemap before the grdimage, so that the grdimage can take the projection and regional extent defined in the basemap. The projection of the basemap can be Mercator cylindrical, Equidistant cylindrical, etc.

There are different projections available at https://docs.generic-mapping-tools.org/dev/cookbook/options.html#coordinate-transformations-and-map-projections-the-j-option

If you want to display UTM coordinates
You can use grdinfo to find the extent of grid_UTM.tif. Plot the grdimage using the UTM projection and extent obtained from grdinfo. You can forget about the basemap; however, if the basemap is needed, it’s better to plot it before grdimage.

5

Could you show me an example using pygmt?

6

Reading the manuals normally helps.