next up previous contents index
Next: 7.29 Gridding spherical surface Up: 7. Creating GMT Graphics Previous: 7.27 Plotting Sandwell/Smith Mercator   Contents   Index


7.28 Mixing UTM and geographic data sets

Next, we present a similar case: We wish to plot a data set given in UTM coordinates 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 there are a few things you must consider:

  1. You need to know the lower left and upper right UTM coordinates of your map. Given the UTM zone you can use mapproject to recover the lon/lat of those two points. Conversely, if you instead know the lon/lat corners then you need to convert those to UTM coordinates. You now have the ability to specify two domains with the -R setting: The linear UTM meter domain when plotting UTM data and the geographic domain (remember to use the rectangular variant of -R that ends with the modifier r) when plotting lon/lat data.
  2. Make sure you use the same scale (and not width) with both the linear and UTM projection.




#!/bin/sh
#               GMT EXAMPLE 28
#
# Purpose:      Illustrates how to mix UTM data and UTM projection
# GMT progs:    makecpt, grdgradient, grdimage, grdinfo, pscoast, pstext, mapproject
# Unix progs:   rm, cut, grep, $AWK
#
ps=example_28.ps

# Get intensity grid and set up a color table
grdgradient Kilauea.utm.nc -Nt1 -A45 -GKilauea.utm_i.nc
makecpt -Ccopper -T0/1500/100 -Z > Kilauea.cpt
# Save min/max UTM coordinates with enough precision
grdinfo Kilauea.utm.nc --D_FORMAT=%.10g -C > tmp.txt
# Use inverse UTM projection to determine the lon/lat of the lower left and upper right corners
LL=`cut -f2,4 tmp.txt | mapproject -Ju5Q/1:1 -F -C -I --OUTPUT_DEGREE_FORMAT=ddd:mm:ss.x | \
        $AWK '{printf "%s/%s\n", $1, $2}'`
UR=`cut -f3,5 tmp.txt | mapproject -Ju5Q/1:1 -F -C -I --OUTPUT_DEGREE_FORMAT=ddd:mm:ss.x | \
        $AWK '{printf "%s/%s\n", $1, $2}'`
# Lay down the UTM topo grid using a 1:17,000 scale
grdimage Kilauea.utm.nc -IKilauea.utm_i.nc -CKilauea.cpt -Jx1:170000 -P -K -B5000g5000WSne \
        -U"Example 28 in Cookbook" --D_FORMAT=%.10g --ANNOT_FONT_SIZE_PRIMARY=9 \
        --GRID_CROSS_SIZE_PRIMARY=0.1i > $ps
# Overlay geographic data and coregister by using correct region and projection with the same scale
pscoast -R$LL/${UR}r -Ju5Q/1:170000 -O -K -Df+ -Slightblue -W0.5p -B5mg5mNE \
        --ANNOT_FONT_SIZE_PRIMARY=12 --PLOT_DEGREE_FORMAT=ddd:mmF >> $ps
psbasemap -R -J -O -K --ANNOT_FONT_SIZE_PRIMARY=9 -Lf155:07:30W/19:15:40N/19:23N/5k+l1:17,000+u \
        --LABEL_FONT_SIZE=10 >> $ps
echo 155:16:20W 19:26:20N 12 0 1 CB KILAUEA | pstext -R -J -O >> $ps
# Clean up

rm -f Kilauea.utm_i.nc Kilauea.cpt tmp.txt .gmt*


Our script illustrates how we would plot a UTM grid of elevations near Kilauea volcano on the Big Island of Hawaii. Given we are in UTM zone 5Q, the script determines the geographic coordinates of the lower left and upper right corner of the UTM grid, then uses that region when overlaying the coastline and light blue ocean. We place a scale bar and label Kilauea crater to complete the figure.

Figure 7.28: Mixing UTM and geographic data sets requires knowledge of the map region domain in both UTM and lon/lat coordinates and consistent use of the same map scale.
\includegraphics[scale=0.5]{scripts/example_28}


next up previous contents index
Next: 7.29 Gridding spherical surface Up: 7. Creating GMT Graphics Previous: 7.27 Plotting Sandwell/Smith Mercator   Contents   Index
Paul Wessel 2011-02-27