Mercator map centered on the reference point

Q&A
1

Hi all;

I’ve plot a map on a Mercator projection using earthquake data. What i want to achieve is to get the map centered on single earthquake coordinate. I’ve try to calculate the bounding box using GMT info but it just not getting the map centered. Below is the code that I’ve used to plot the map. Any help would be appreciated.

#!/bin/sh

#
# GMT set up
#
export GMT_DATADIR=`pwd`
export GMT_CACHEDIR=`pwd`

gmt set FORMAT_GEO_MAP=dddF \
    PS_MEDIA A4 \
    MAP_FRAME_PEN=dimgray \
    MAP_FRAME_WIDTH=0.1c \
    MAP_TITLE_OFFSET=1c \
    MAP_ANNOT_OFFSET=0.1c \
    MAP_TICK_PEN_PRIMARY=thinner,dimgray \
    MAP_GRID_PEN_PRIMARY=thin,white \
    MAP_GRID_PEN_SECONDARY=thinnest,white \
    FONT_TITLE=12p,Palatino-Roman,black \
    FONT_ANNOT_PRIMARY=7p,0,dimgray \
    FONT_LABEL=7p,0,dimgray

varcurrentDate=`date +"%Y-%m-%d %H:%M:%S %p %Z"`
varMapPSFile="map-test-admis.ps"
#
# Earthquake Data
#   
myEQOriginTimeUTC="2022-11-18 7:15:45 UTC"
myEQLat=6.047
myEQLong=116.590
myEQLocation="Ranau, Sabah"
myEMagnitude="6.0"
myEQDepth="9"

#
# Map Setting
#
varMapSize="7.0";            # Map Size in Inch
varMapBoundZoom="5/5/5/5";
varMapTick="1";
#
# Get Map Boundaries from Location Coordinates 
#
varInfo="-C -I${varMapBoundZoom} -rp"

minMax=`echo $myEQLong $myEQLat | gmt info ${varInfo}`

xMin=`echo $minMax | awk '{ print $1 }'`
xMax=`echo $minMax | awk '{ print $2 }'`
yMin=`echo $minMax | awk '{ print $3 }'`
yMax=`echo $minMax | awk '{ print $4 }'`

varMapBound="$xMin/$xMax/$yMin/$yMax"

midlon=`gmt gmtmath -Q $xMin $xMax ADD 2 DIV =`
midlat=`gmt gmtmath -Q $yMin $yMax ADD 2 DIV =`

varMapBase="-Bpx${varMapTick}f${varMapTick} -Bpy${varMapTick}g${varMapTick}f${varMapTick} -Bsxg${varMapTick} -Bsy "

varMapProj="M${myEQLong}/${myEQLat}/${varMapSize}i"
echo "$minMax | $xMin $xMax $yMin $yMax | Lat:$myEQLat Lon:$myEQLong | Center Point: $midlat $midlon | Map Proj: $varMapProj | Map Bound: $varMapBound Base:$varMapBase" 
#
# Start plot map
#
gmt pscoast -R$varMapBound -J$varMapProj -P -N1/thin,olivedrab1 -Na -Df -Slightblue -W0.01 -K > $varMapPSFile 

gmt psbasemap -R -J \
    --MAP_FRAME_AXES=WEsN \
    --FORMAT_GEO_MAP=ddd:mm:ssF \
     ${varMapBase}  \
    --MAP_TITLE_OFFSET=0.5i \
    --FONT_ANNOT_PRIMARY=10p,0,black \
    --FONT_LABEL=12p,25,black \
    --FONT_TITLE=18p,25,black \
    -O -K >> $varMapPSFile 
#
# Plot Earthquake Coordinate
#
gmt psxy -R -J -Wfaint -Sa0.3i -W0.5p -O -Gred << END >> $varMapPSFile
$myEQLong $myEQLat
END
#
# Output image file to jpeg
#
gmt psconvert $varMapPSFile -A1.7c -Tj
2

varMapBound=$(echo $myEQLong $myEQLat | gmt info -I5 -D)

use $varMapBound replace -R$varMapBound

3

I did not read through your script, but here are some nice options when you want to define a map center, and draw a map around it.

  • $ gmt grdimage @earth_relief_02m -pdf test -R-500/500/-500/500+uk -JM10/60/15c -B

    (10,60) is the center and draw 500km to the left, up, right, down.

  • $ gmt pscoast -pdf square -Bafg -Ggray -JS0/60/10c -R400+uk

    (0,60) is center, and a circle of r = 400k, -Rradius+uunit

  • $ gmt pscoast -pdf rectangular -Bafg -Ggray -JS0/60/10c -R400/800+uk

    Same as above but -Rhalfwidth/halfheight+uunit.

Bad explanations - try them out and you’ll understand.

4

Thanks @Alpiner for the suggestion it works. Before your suggestion i try to use this crude implementation to get min and max by coordinates:

xMin="$(bc <<<"$myEQLong - $varIncrement")"
xMax="$(bc <<<"$myEQLong + $varIncrement")"
yMin="$(bc <<<"$myEQLat - $varIncrement")"
yMax="$(bc <<<"$myEQLat + $varIncrement")"