One would have thought that Gauss should have dealt a death blow to new map projections when he proved a flat map of a sphere could not be distortion free. Because after all, with all the map projections known about up to the time of Gauss, surely one of them would be OK for whatever task people wanted to do.
But on the contrary, there seems to be no end of new projections to try and display data for a certain area better than before. Take this example from Snyder:
Almost perfect annotations can be added to arbitrary map projections if one computes the grid lines in the arbitrary projection coordinates, and then simply plots the grid lines as quoted lines, with latitude or longitude labels placed where the grid lines intersect a box near the border of the map. Quoted line functionality seems to have been devised for exactly this scenario.
Here’s all that’s needed in shell script (you’ll need to adjust the range of longitude and latitudes for the grid lines to suit the map you’re making). It would be trivial to extend this function to parse a GMT style -R argument and work out what lines to compute from this.
function grid(){
for lat in $(seq -80 10 80)
do
echo "> ${lat}"
gmt math -T70/330/1 ${lat} = | gmt mapproject -J"${proj}"
done
for lon in $(seq 70 10 330)
do
echo "> ${lon}"
gmt math -o1,0 -T-80/80/1 ${lon} = | gmt mapproject -J"${proj}"
done
}
And then you use this function like this (where ${la}, ${ba} etc are the coordinates of a box just inside the frame of the map (in projected coordinates). In theory you can use two letter codes (LB, LT etc) for the corners too, but I ran into problems with this, so put my box 1 cm inside the frame instead of right on it.
grid | gmt plot -Jx${scale} -R${proj_region} -Wthinner,gray60 \
-Sql${la}/${ba}/${la}/${ta},${la}/${ta}/${ra}/${ta},${ra}/${ta}/${ra}/${ba}:+f12,Helvetica,gray60+Lh
Here is an example (I added a bit more complexity by converting angles like -60 to 60S etc).
There are situations this method doesn’t handle, such as when the frame of the plot is not rectangular (I can envisage how non-rectangular frames could be handled). And issues with labels overlapping each other can occur. But it’s a pretty good start … at least I reckon so.
Edit: overlapping labels and non-rectangular frames has been “solved” with the -Sqx... feature. Here is a current draft version of my map where I made a series of lines around a non rectangular frame for the -Sqx... feature to label.
It’s being able to hack GMT into doing things like this which makes it so powerful.
I think this need to go to showcase with a full script 
I’ll tidy up the script as much as I can. It’s bash … mainly because I know how to do GMT with the shell. And learning the newer ways is a gradual process. The key “breakthrough” for me (I’m sure others were aware of this though) was finding that quoted lines solved many of the annotation problems quite easily - in fact, quoted lines solved problems I’d never even imagined could be done short of working on the C code.
Here is a script which plots a coastline map for an arbitrary map projection supported by proj, and adds “nice” annotations around the borders of the map.
#!/usr/bin/env bash
proj="+proj=robin +ellps=GRS80 +lon_0=-165"
scale="1:50000000"
left=-8000000
right=8000000
bottom=-8000000
top=8000000
proj_region="${left}/${right}/${bottom}/${top}"
function annot_frame(){
# Create an interior frame inside the actual map frame where annotations
# will be placed. If an argument is provided, this is the spacing between
# the border and the annotations (default 1 cm).
local spc=${1:-0.01}
local la=$(bc -l -e "${left}+${spc}*${scale:2}")
local ra=$(bc -l -e "${right}-${spc}*${scale:2}")
local ba=$(bc -l -e "${bottom}+${spc}*${scale:2}")
local ta=$(bc -l -e "${top}-${spc}*${scale:2}")
cat <<- END
>
${la} ${ba}
${la} ${ta}
>
${la} ${ta}
${ra} ${ta}
>
${ra} ${ta}
${ra} ${ba}
>
${ra} ${ba}
${la} ${ba}
END
}
function grid(){
# Calculate gridlines for the whole globe projected onto the map
# projection. Optionally pass the gridline spacing as an argument
# (default 10 degrees).
spc=${1:-10}
for lat in $(seq -90 ${spc} 90)
do
awk -v l=${lat} \
'BEGIN{printf("> %s\n",(l>0)?l"\260N":(l<0)?-l"\260S":l"\260");}'
gmt math -T70/330/0.1 ${lat} = | gmt mapproject -J"${proj}"
done
for lon in $(seq 0 ${spc} 360)
do
awk -v l=${lon} \
'BEGIN{printf("> %s\n",(l>180)?(360-l)"\260W":(l<180)?l"\260E":l"\260");}'
gmt math -o1,0 -T-80/80/0.1 ${lon} = | gmt mapproject -J"${proj}"
done
}
function split_polygon(){
# Split polygons which contain sequential points more than a specified
# distance apart (default: 1000 km) into separate polygons.
gawk -v tol=${1:-1000000} \
'(NR==1){
xp = $1;
yp = $2;
}
(NR>1){
if ((($1-xp)^2 + ($2-yp)^2) > tol^2){
printf("> Split\n%s\n",$0);
} else {
printf("%s\n",$0);
}
xp = $1;
yp = $2;
}'
}
gmt begin testmap pdf
gmt coast -Di -Rg -M -W | gmt mapproject -J"${proj}" | split_polygon | \
gmt plot -Jx${scale} -R${proj_region} -Wthinnest,black -BWSEN
annot_frame 0.005 > frame.txt
grid 10 | gmt plot -Jx${scale} -R${proj_region} -Wthinner,gray10 \
-Sqxframe.txt+r1c:+f8,Helvetica,gray25+Lh
gmt end
How it works:
The map projection is set in the ${proj} variable. In this example it is set to Robinson’s map projection centred at 165 West. The rectangular boundaries of the map are specified by the varialbles ${left}, ${right}, ${top} and ${bottom}. The units of these boundaries are projected metres. The map scale is set by ${scale}.
Firstly we dump the coastlines from the in-built GSHHG coastline dataset, and project it to the chosen map projection. GSHHG needn’t be used - more modern datasets such as OpenStreetMap could be used at this point.
After we have the projected coastlines, we need to split those coastline polygons which wrap around from the right hand edge of our map to the left hand edge. This is important. If we don’t do this, the coastline looks horrible with lines connecting portions of the coast on the right hand edge of the map with portions on the left hand edge. The way the splitting is done is to see if any sequential points in the projected coastline exceed a critical distance (default 1000 km). If this happens, it’s almost certainly points wrapping around the edge of the map. I wrote an awk script to do this splitting (split_polygons), but I wouldn’t be surprised if a GMT command to do this exists - I just don’t know what it is called.
The coastlines are then plotted with gmt plot.
To get nice grid lines and annotations, we draw the meridians and parallels as annotated lines (-Sq). One of the features of the annotated lines is that annotations can be placed where the line intersects other lines (-Sqxfile.txt option) specified in a file. This is the key to nice annotations. We compute a series of lines around the edge of our map where we want our annotations to lie. The example script has the “annotation frame” a small distance inside the edge of the map. However, arbitrary “frames” can be made - the limit is your imagination. The +r modified for annotated lines ensures no annotations overlap each other.
The grid lines, which are annotated are computed using the grid shell function. The annotation “frame” is computed with the annot_frame function.
And that’s all there is to it. There are gotchas with certain map projections, but this is the basic procedure for making nicely annotated and gridded maps on any projection supported by proj.
See “coast -gD”
Magic. I knew there’d be a solution, just didn’t know what it was.
Using -gd (I used it in the plot call rather than coast, because of the way I’m specifying the region in projected coordinates rather than geographical coordinates), the shell script is even simpler.
#!/usr/bin/env bash
proj="+proj=robin +ellps=GRS80 +lon_0=-165"
scale="1:50000000"
left=-8000000
right=8000000
bottom=-8000000
top=8000000
proj_region="${left}/${right}/${bottom}/${top}"
function annot_frame(){
# Create an interior frame inside the actual map frame where annotations
# will be placed. If an argument is provided, this is the spacing between
# the border and the annotations (default 1 cm).
local spc=${1:-0.01}
local la=$(bc -l -e "${left}+${spc}*${scale:2}")
local ra=$(bc -l -e "${right}-${spc}*${scale:2}")
local ba=$(bc -l -e "${bottom}+${spc}*${scale:2}")
local ta=$(bc -l -e "${top}-${spc}*${scale:2}")
cat <<- END
>
${la} ${ba}
${la} ${ta}
>
${la} ${ta}
${ra} ${ta}
>
${ra} ${ta}
${ra} ${ba}
>
${ra} ${ba}
${la} ${ba}
END
}
function grid(){
# Calculate gridlines for the whole globe projected onto the map
# projection. Optionally pass the gridline spacing as an argument
# (default 10 degrees).
spc=${1:-10}
for lat in $(seq -90 ${spc} 90)
do
awk -v l=${lat} \
'BEGIN{printf("> %s\n",(l>0)?l"\260N":(l<0)?-l"\260S":l"\260");}'
gmt math -T70/330/0.1 ${lat} = | gmt mapproject -J"${proj}"
done
for lon in $(seq 0 ${spc} 360)
do
awk -v l=${lon} \
'BEGIN{printf("> %s\n",(l>180)?(360-l)"\260W":(l<180)?l"\260E":l"\260");}'
gmt math -o1,0 -T-80/80/0.1 ${lon} = | gmt mapproject -J"${proj}"
done
}
gmt begin testmap pdf,png
gmt coast -Di -Rg -M -W | gmt mapproject -J"${proj}" | \
gmt plot -Jx${scale} -R${proj_region} -Wthinnest,black -BWSEN -gd1000000
annot_frame 0.005 > frame.txt
grid 10 | gmt plot -Jx${scale} -R${proj_region} -Wthinner,gray10 \
-Sqxframe.txt+r1c:+f8,Helvetica,gray25+Lh
gmt end
Does anyone else have a problem with the bc command?
bc: invalid option -- 'e'
Could be the version of bc on my Mac is different than on Linux (for example).
Anyhow, all I’m doing is computing a rectangle (in projected coordinates) which lies just inside the edge of the map. awk would able to do the job too. This is a good example of the problems one encounters with shell scripts. If this was in Julia or Python it wouldn’t be a problem.
Could be the version of bc on my Mac is different than on Linux (for example).
Probably. My bc doesn’t have an -e option, and yours presumably do.
I’m starting to learn zsh (it’s been the default of the Mac for some time, but I’m still using bash). zsh can do floating point arithmetic which eliminates the need for bc
This zsh script seems to work on my Mac:
#!/usr/bin/env zsh
proj="+proj=robin +ellps=GRS80 +lon_0=-165"
scale="1:50000000"
left=-8000000
right=8000000
bottom=-8000000
top=8000000
proj_region="${left}/${right}/${bottom}/${top}"
function annot_frame(){
# Create an interior frame inside the actual map frame where annotations
# will be placed. If an argument is provided, this is the spacing between
# the border and the annotations (default 1 cm).
local spc=${1:-0.01}
cat <<- END
>
$((${left}+${spc}*${scale:2})) $((${bottom}+${spc}*${scale:2}))
$((${left}+${spc}*${scale:2})) $((${top}-${spc}*${scale:2}))
>
$((${left}+${spc}*${scale:2})) $((${top}-${spc}*${scale:2}))
$((${right}-${spc}*${scale:2})) $((${top}-${spc}*${scale:2}))
>
$((${right}-${spc}*${scale:2})) $((${top}-${spc}*${scale:2}))
$((${right}-${spc}*${scale:2})) $((${bottom}+${spc}*${scale:2}))
>
$((${right}-${spc}*${scale:2})) $((${bottom}+${spc}*${scale:2}))
$((${left}+${spc}*${scale:2})) $((${bottom}+${spc}*${scale:2}))
END
}
function grid(){
# Calculate gridlines for the whole globe projected onto the map
# projection. Optionally pass the gridline spacing as an argument
# (default 10 degrees).
local spc=${1:-10}
for ((lat=-90; lat<=90; lat+=${spc}))
do
awk -v l=${lat} \
'BEGIN{printf("> %s\n",(l>0)?l"\260N":(l<0)?-l"\260S":l"\260");}'
gmt math -T0/360/0.1 ${lat} = | gmt mapproject -J"${proj}"
done
for ((lon=0; lon<360; lon+=${spc}))
do
awk -v l=${lon} \
'BEGIN{printf("> %s\n",(l>180)?(360-l)"\260W":(l<180)?l"\260E":l"\260");}'
gmt math -o1,0 -T-90/90/0.1 ${lon} = | gmt mapproject -J"${proj}"
done
}
gmt begin testmap pdf,png
gmt coast -Di -Rg -M -W | gmt mapproject -J"${proj}" | \
gmt plot -Jx${scale} -R${proj_region} -Wthinnest,black -BWSEN -gd1000000
annot_frame 0.005 > frame.txt
grid 10 | gmt plot -Jx${scale} -R${proj_region} -Wthinner,gray10 \
-Sqxframe.txt+r1c:+f8,Helvetica,gray25+Lh
gmt end
I was able to run the zsh-script. Looks great!
Below is zsh code with no dependencies other than GMT.
#!/usr/bin/env zsh
proj="+proj=robin +ellps=GRS80 +lon_0=-165"
scale="1:50000000"
left=-8000000
right=8000000
bottom=-8000000
top=8000000
proj_region="${left}/${right}/${bottom}/${top}"
function annot_frame(){
# Create an interior frame inside the actual map frame where annotations
# will be placed. If an argument is provided, this is the spacing between
# the border and the annotations (default 1 cm).
local spc=${1:-0.01}
cat <<- END
>
$((${left}+${spc}*${scale:2})) $((${bottom}+${spc}*${scale:2}))
$((${left}+${spc}*${scale:2})) $((${top}-${spc}*${scale:2}))
>
$((${left}+${spc}*${scale:2})) $((${top}-${spc}*${scale:2}))
$((${right}-${spc}*${scale:2})) $((${top}-${spc}*${scale:2}))
>
$((${right}-${spc}*${scale:2})) $((${top}-${spc}*${scale:2}))
$((${right}-${spc}*${scale:2})) $((${bottom}+${spc}*${scale:2}))
>
$((${right}-${spc}*${scale:2})) $((${bottom}+${spc}*${scale:2}))
$((${left}+${spc}*${scale:2})) $((${bottom}+${spc}*${scale:2}))
END
}
function grid(){
# Calculate gridlines for the whole globe projected onto the map
# projection. Optionally pass the gridline spacing as an argument
# (default 10 degrees).
local spc=${1:-10}
for (( lat=-90.0; lat<=90.0; lat+=${spc} ))
do
if (( ${lat} > 0 )); then
printf "> %.f\260N\n", ${lat}
elif (( ${lat} < 0 )); then
printf "> %.f\260S\n", $(( -1*${lat} ))
else
printf "> %.f\260\n" ${lat}
fi
gmt math -T0/360/0.1 ${lat} = | gmt mapproject -J"${proj}"
done
for (( lon=0.0; lon<360.0; lon+=${spc} ))
do
if (( ${lon} < 180 )); then
printf "> %.f\260E\n", ${lon}
elif (( ${lon} > 180 )); then
printf "> %.f\260W\n", $(( 360.0-${lon} ))
else
printf "> %.f\260\n" ${lon}
fi
gmt math -o1,0 -T-90/90/0.1 ${lon} = | gmt mapproject -J"${proj}"
done
}
gmt begin testmap pdf,png
gmt coast -Di -Rg -M -W | gmt mapproject -J"${proj}" | \
gmt plot -Jx"${scale}" -R"${proj_region}" -Wthinnest,black -BWSEN -gd1000000
annot_frame 0.005 > frame.txt
grid 10 | gmt plot -Jx"${scale}" -R"${proj_region}" -Wthinner,gray10 \
-Sqxframe.txt+r1c:+f8,Helvetica,gray25+Lh
gmt end
Tim, back to this. You can make your script slightly friendlier by using the new mapproject -Ws and thus not require that user pass in a map scale (you could also use limits in geographical coordinates). And does it work for the Lee projection? I had several headaches with the 0 and 20N parallels that reenter in the plot domain. Other shit that annoyed and made loose lots of time was on how to workaround a gdalwarp +proj=lee_os bug.
But I can now do
G, cl = leepacific("@earth_relief_20m_p");
grdimage(G, shade=true, plot=(data=cl,), cmap=:geo, B=:none)
plotgrid!(G, show=true)
and, with GMT master and using your quoted lines trick
grdimage(G, shade=true, plot=(data=cl,), cmap=:geo, B=:bare)
plotgrid!(G, inner=true, show=true)
The frame issues (grid lines, annotations etc. ) with this projection looks fixed in your upper figure - is this possible in normal (ie. non-Julia) GMT as well?
If you do the kind of programming Tim shows in his script, yes. What I do is to compute the intersections between the frame and the grid lines using gmtspatial and then use them to plot custom axis.
So the magic that you do, happens under the hood with Julia?
(All this with 2-3 lines of Julia?!)
Yes, and with the added bonus that you can also control fig size, line colors, styles, etc… All of that with a human readable/understandable syntax.