img2mercgrd
img2mercgrd - Extract region of img, preserving Mercator,
save as grd
SYNOPSIS
img2mercgrd imgfile -Ggrdfile -Rwest/east/south/north
-Ttype [ -Nnavg ] [ -Sscale ] [ -V ] [ -mminutes ] [
-xmaxlon ] [ -yminlat/maxlat ]
DESCRIPTION
img2mercgrd reads an img format file and creates a grd
file. The Spherical Mercator projection of the img file is
preserved, so that the region -R set by the user is modi
fied slightly; the modified region corresponds to the
edges of pixels [or groups of navg pixels]. The grdfile
header is set so that the x and y axis lengths represent
distance from the west and south edges of the image, mea
sured in user default units, with -Jm1 and the adjusted
-R. By setting the gmtdefaults ELLIPSOID = Sphere, the
user can make overlays with the adjusted -R so that they
match. See examples below. The adjusted -R is also writ
ten in the grdheader remark, so it can be found later. The
-Ttype selects all data or only data at constrained pix
els, and can be used to create a grid of 1s and 0s indi
cating constraint locations. The output grd file is pixel
registered; it inherits this from the img file.
imgfile
An img format file such as the marine gravity or
seafloor topography fields estimated from satellite
altimeter data by Sandwell and Smith. If the user
has set an environment variable GMT_IMGDIR, then
img2mercgrd will try to find imgfile in
$GMT_IMGDIR; else it will try to open imgfile
directly.
-G grdfile is the name of the output grdfile.
-R west, east, south, and north specify the Region of
interest. To specify boundaries in degrees and min
utes, use the dd:mm format.
-T type handles the encoding of constraint informa
tion. type = 0 indicates that no such information
is encoded in the img file (used for pre-1995 ver
sions of the gravity data; all more recent files do
not support this choice) and gets all data. type >
0 indicates that constraint information is encoded
(1995 and later (current) versions of the img
files) so that one may produce a grd file as fol
lows: -T1 gets data values at all points, -T2 gets
data values at constrained points and NaN at inter
polated points; -T3 gets 1 at constrained points
-N Average the values in the input img pixels into
navg by navg squares, and create one output pixel
for each such square. If used with -T3 it will
report an average constraint between 0 and 1. If
used with -T2 the output will be average data value
or NaN according to whether average constraint is >
0.5. navg must evenly divide into the dimensions
of the imgfile in pixels. [Default 1 does no aver
aging].
-S Multiply the img file values by scale before stor
ing in grd file. [Default is 1.0]. (img topo files
are stored in (corrected) meters; gravity files in
mGal*10; vertical deflection files in microradi
ans*10, vertical gravity gradient files in
Eotvos*10. Use -S0.1 for those files)
-V Selects verbose mode, which will send progress
reports to stderr [Default runs "silently"]. Par
ticularly recommended here, as it is helpful to see
how the coordinates are adjusted.
-m Indicate minutes as the width of an input img pixel
in minutes of longitude. [Default is 2.0]
-x Indicate maxlon as the maximum longitude extent of
the input img file. Versions since 1995 have had
maxlon = 360.0, while some earlier files had maxlon
= 390.0. [Default is 360.0]
-y Indicate minlat/maxlat as the latitude extent of
the input img file. All versions to date have used
-72.006/72.006. [Default is -72.006/72.006]
EXAMPLES
To extract data in the region -R-40/40/-70/-30 from
world_grav.img.7.2 try
img2mercgrd world_grav.img.7.2 -Gmerc_grav.grd
-R-40/40/-70/-30 -T1 -V
Note that the -V option tells us that the range was
adjusted to -R-40/40/-70.0004681551/-29.9945810754 We can
also use grdinfo to find that the grd file header shows
its region to be -R0/80/0/67.9666667 This is the range of
x,y we will get from a Spherical Mercator projection using
-R-40/40/-70.0004681551/-29.9945810754 and -Jm1. Thus, to
take ship.lonlatgrav and use it to sample the
merc_grav.grd, we can do this:
gmtset ELLIPSOID Sphere
mapproject -R-40/40/-70.0004681551/-29.9945810754 -Jm1
latgravsat
It is recommended to use the above method of projecting
and unprojecting the data in such an application, because
then there is only one interpolation step (in grdtrack).
If one first tries to convert the grd file to lon,lat and
then sample it, there are two interpolation steps (in con
version and in sampling).
To make a lon,lat grid from the above grid we can use
grdproject merc_grav.grd
-R-40/40/-70.0004681551/-29.9945810754 -Jm1 -I -F -D2m
-Ggrav.grd
In some cases this won't be easy as the -R in the two
coordinate systems won't align well. When this happens,
we can also use (in fact, it may be always better to use)
grd2xyz merc_grav.grd | mapproject
-R-40/40/-70.0004681551/-29.994581075 -Jm1 -I | surface
-R-40/40/-70/70 -I2m -Ggrav.grd
To make a Mercator map of the above region, suppose our
gmtdefault MEASURE_UNIT is inch. Then since the above
merc_grav.grd file is projected with -Jm1 it is 80 inches
wide. We can make a map 8 inches wide by using -Jx0.1 on
any map programs applied to this grid (e.g., grdcontour,
grdimage, grdview), and then for overlays which work in
lon,lat (e.g., psxy, pscoast) we can use the above
adjusted -R and -Jm0.1 to get the two systems to match up.
However, we can be smarter than this. Realizing that the
input img file had pixels 2.0 minutes wide (or checking
the nx and ny with grdinfo merc_grav.grd) we realize that
merc_grav.grd used the full resolution of the img file and
it has 2400 by 2039 pixels, and at 8 inches wide this is
300 pixels per inch. We decide we don't need that many and
we will be satisfied with 100 ' pixels per inch, so we
want to average the data into 3 by 3 squares. (If we want
a contour plot we will probably choose to average the data
much more (e.g. 6 by 6) to get smooth contours.) Since
2039 isn't divisible by 3 we will get a different adjusted
-R this time:'
img2mercgrd world_grav.img.7.2 -Gmerc_grav_2.grd
-R-40/40/-70/-30 -T1 -N3 -V
This time we find the adjusted region is
-R-40/40/-70.023256525/-29.9368261101 and the output is
800 by 601 pixels, a better size for us. Now we can create
an artificial illumination file for this using grdgradi
and if we also have a cpt file called "grav.cpt" we can
create a color shaded relief map like this:
grdimage merc_grav_2.grd -Iillum.grd -Cgrav.cpt -Jx0.1 -K
> map.ps
psbasemap -R-40/40/-70.023256525/-29.9368261101 -Jm0.1
-Ba10 -O >> map.ps
Suppose you want to obtain only the constrained data val
ues from an img file, in lat/lon coordinates. Then run
img2mercgrd with the -T2 option, use grd2xyz to dump the
values, pipe through grep -v NaN to eliminate NaNs, and
pipe through mapproject with the inverse projection as
above.
SEE ALSO
gmt(l), grdproject(l), mapproject(l)
Man(1) output converted with
man2html