(28) Mixing UTM and geographic data sets¶
Next, we present a similar case: We wish to plot a data set given in UTM coordinates (meter) and want it to be properly registered with overlying geographic data, such as coastlines or data points. The mistake many GMT rookies make is to specify the UTM projection with their UTM data. However, that data have already been projected and is now in linear meters. The only sensible way to plot such data is with a linear projection, yielding a UTM map. In this step one can choose to annotate or tick the map in UTM meters as well. To plot geographic (lon/lat) data on the same map you simply have to specify the region using the UTM meters but supply the actual UTM projection parameters. Make sure you use the same scale with both the linear and UTM projection.
Our script illustrates how we would plot a UTM grid (with coordinates in meters) of elevations near Kilauea volcano on the Big Island of Hawaii and overlay geographical information (with longitude, latitude coordinates). We first lay down the UTM grid using the linear projection. Then, given we are in UTM zone 5Q, we use the UTM domain and the UTM projection when overlaying the coastline and light blue ocean. We do some trickery by converting the UTM domain to km so that we can add custom annotations to the map. Finally, we place a scale bar and label Kilauea crater to complete the figure.
#!/bin/bash
# GMT EXAMPLE 28
# $Id$
#
# Purpose: Illustrates how to mix UTM data and UTM gmt projection
# GMT modules: makecpt, grdimage, grdmath, pscoast, pstext,
# Unix progs: rm, echo
#
ps=example_28.ps
# Set up a color table
gmt makecpt -Ccopper -T0/1500 > Kilauea.cpt
# Lay down the UTM topo grid using a 1:16,000 scale
gmt grdimage Kilauea.utm.nc -I+a45+nt1 -CKilauea.cpt -Jx1:160000 -P -K \
--FORMAT_FLOAT_OUT=%.10g --FONT_ANNOT_PRIMARY=9p > $ps
# Overlay geographic data and coregister by using correct region and gmt projection with the same scale
gmt pscoast -RKilauea.utm.nc -Ju5Q/1:160000 -O -K -Df+ -Slightblue -W0.5p -B5mg5m -BNE \
--FONT_ANNOT_PRIMARY=12p --FORMAT_GEO_MAP=ddd:mmF >> $ps
echo 155:16:20W 19:26:20N KILAUEA | gmt pstext -R -J -O -K -F+f12p,Helvetica-Bold+jCB >> $ps
gmt psbasemap -R -J -O -K --FONT_ANNOT_PRIMARY=9p -LjRB+c19:23N+f+w5k+l1:16,000+u+o0.2i \
--FONT_LABEL=10p >> $ps
# Annotate in km but append ,000m to annotations to get customized meter labels
gmt psbasemap -RKilauea.utm.nc+Uk -Jx1:160 -B5g5+u"@:8:000m@::" -BWSne -O --FONT_ANNOT_PRIMARY=10p \
--MAP_GRID_CROSS_SIZE_PRIMARY=0.1i --FONT_LABEL=10p >> $ps
# Clean up
rm -f Kilauea.cpt tmp.txt