(27) Plotting Sandwell/Smith Mercator img grids — GMT 5.1.1 documentation

### Navigation

- index
- next
- previous
- Project Home »
- GMT 5.1.1 documentation »
- The GMT historical collection »
- Search:

# (27) Plotting Sandwell/Smith Mercator img grids

Next, we show how to plot a data grid that is distributed in projected form. The gravity and predicted bathymetry grids produced by David Sandwell and Walter H. F. Smith are not geographical grids but instead given on a spherical Mercator grid. The GMT supplement imgsrc has tools to extract subsets of these large grids. If you need to make a non-Mercator map then you must extract a geographic grid using img2grd and then plot it using your desired map projection. However, if you want to make a Mercator map then you can save time and preserve data quality by avoiding to re-project the data set twice since it is already in a Mercator projection. This example shows how this is accomplished. We use the -M option in img2grd [1] to pull out the grid in Mercator units (i.e., do not invert the Mercator projection) and then simply plot the grid using a linear projection with a suitable scale (here 0.25 inches per degrees of longitude). To overlay basemaps and features that has geographic longitude/latitude coordinates we must remember two key issues:

1. This is a spherical Mercator grid so we must use –PROJ_ELLIPSOID=Sphere with all commands that involve projections (or use gmtset to change the setting).
2. Select Mercator projection and use the same scale that was used with the linear projection.

This map of the Tasman Sea shows the marine gravity anomalies with land painted black. A color scale bar was then added to complete the illustration.

    #!/bin/bash
    #		GMT EXAMPLE 27
    #
    # Purpose:	Illustrates how to plot Mercator img grids
    # GMT progs:	makecpt, mapproject, grdgradient, grdimage, grdinfo, pscoast
    # GMT supplement: img2grd (to read Sandwell/Smith img files)
    # Unix progs:	rm, grep, $AWK
    #
    ps=example_27.ps

    # Gravity in tasman_grav.nc is in 0.1 mGal increments and the grid
    # is already in projected Mercator x/y units.
    # First get gradients.

    gmt grdgradient tasman_grav.nc -Nt1 -A45 -Gtasman_grav_i.nc

    # Make a suitable cpt file for mGal

    gmt makecpt -T-120/120/240 -Z -Crainbow > grav.cpt

    # Since this is a Mercator grid we use a linear projection

    gmt grdimage tasman_grav.nc=ns/0.1 -Itasman_grav_i.nc -Jx0.25i -Cgrav.cpt -P -K > $ps

    # Then use gmt pscoast to plot land; get original -R from grid remark
    # and use Mercator gmt projection with same scale as above on a spherical Earth

    R=`gmt grdinfo tasman_grav.nc | grep Remark | $AWK '{print $NF}'`

    gmt pscoast $R -Jm0.25i -Ba10f5 -BWSne -O -K -Gblack --PROJ_ELLIPSOID=Sphere \
    	-Cwhite -Dh+ --FORMAT_GEO_MAP=dddF >> $ps

    # Put a color legend on top of the land mask justified with 147E,31S

    pos=`echo 147E 31S | gmt mapproject -R -J -Di --PROJ_ELLIPSOID=Sphere | \
    	$AWK '{printf "%si/%si\n", $1, $2}'`
    gmt psscale -D$pos/2i/0.15i -Cgrav.cpt -Bx50f10 -By+lmGal -I -O -T+gwhite+p1p >> $ps

    # Clean up

    rm -f grav.cpt *_i.nc

Plotting Sandwell/Smith Mercator img grids.

[1] | You could also use img2grd directly - your only option under DOS

### Navigation

- index
- next
- previous
- Project Home »
- GMT 5.1.1 documentation »
- The GMT historical collection »
- Search: