Mercator map centered on the reference point

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.


# 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_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 \

varcurrentDate=`date +"%Y-%m-%d %H:%M:%S %p %Z"`
# Earthquake Data
myEQOriginTimeUTC="2022-11-18 7:15:45 UTC"
myEQLocation="Ranau, Sabah"

# Map Setting
varMapSize="7.0";            # Map Size in Inch
# 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 }'`


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 "

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 \
    --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
# Output image file to jpeg
gmt psconvert $varMapPSFile -A1.7c -Tj

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

use $varMapBound replace -R$varMapBound

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.

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")"