Re: MAPM & GMT Package design?

[Nir Krakauer]

Dear Clark,

This would be excellent. Some time ago I wrote an Octave wrapper
(below) to produce and run a GMT script for basic mapping. Feel free
to make use of it in your project.

Nir

_____
## Copyright (C) 2014 Nir Krakauer
##
## This program is free software; you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
## the Free Software Foundation; either version 3 of the License, or
## (at your option) any later version.
##
## This program is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
## GNU General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with this program; If not, see <http://www.gnu.org/licenses/>.

%this function uses GMT to make a global map representing a gridded quantity
%input parameters (only array, lons, lats are required):
%array: (lats * lons) to be gridded
%lons, lats: (degrees E and N) coordinates of array
%range_option: a switch - if zero, the map domain is global; if
nonzero, the map domain is determined to fit the range of the lats and
lons arrays; default: 0
%name: name of the output map file; default: 'example'
%title_use: string to use for the map title; default: 'example'
%colorbar_label: string to use for labeling the color bar; default: none
%file_coords: if a filename is given, plot the coordinates in the file
on the map
%scale: a 2-element vector with the minimum and maximum values on the
color bar - otherwise the minimum and maximum values in array will be
used

function []=map_gmt(array, lons, lats, range_option, name, title_use,
colorbar_label, file_coords, scale)

if ~exist('name', 'var') || isempty(name)
name = 'example';
end

if ~exist('title_use', 'var') || isempty(title_use)
title_use = 'example';
end

if ~exist('colorbar_label', 'var')
colorbar_label = '';
end

if ~exist('range_option', 'var') || ~range_option
range_string = '0/360/-90/90'; %global domain
else
if max(lons) >= 0
range_string = [num2str(min(lons)) '/' num2str(max(lons)) '/'
num2str(min(lats)) '/' num2str(max(lats))];
else %for some reason pscoast can't handle a negative eastern-boundary longitude
range_string = [num2str(min(lons)+360) '/' num2str(max(lons)+360) '/'
num2str(min(lats)) '/' num2str(max(lats))];
end
%range_string = [num2str(min(lons)) '/' num2str(min(lats)) '/'
num2str(max(lons)) '/' num2str(max(lats)) 'r'];
end

n_lons = numel(lons);
n_lats = numel(lats);

if size(array, 1) ~= n_lats
array = array'; %flip so that latitude is the first dimension
end

if ~exist('scale', 'var') || isempty(scale)
minval = nanmin(nanmin(array));
maxval = nanmax(nanmax(array));
else
minval = scale(1);
maxval = scale(2);
end
n_colors = 25;
color_step = (maxval-minval)/n_colors;
label_step = (maxval-minval)/5;

%first create a netcdf file with the information

filename = ['temp.nc'];

nc = netcdf(filename, 'c');

% Define dimensions.

nc('Latitude') = n_lats;
nc('Longitude') = n_lons;

% Define variables.
nc{'Longitude'} = ncdouble('Longitude');
nc{'Latitude'} = ncdouble('Latitude');
nc{'var'} = ncdouble('Latitude', 'Longitude');

%write values
nc{'Longitude'}(:) = lons;
nc{'Latitude'}(:) = lats;
nc{'var'}(:) = array;

nc{'Latitude'}.units = 'degrees_north';                    % Attributes.
nc{'Longitude'}.units = 'degrees_east';

ncclose(nc)

%now write a GMT script to read the file

scriptname = ['gmt_script.txt'];

system(['echo "#!/bin/bash" > ' scriptname]);
system(['echo "export PROJ=\"-JX6i/4i\"" >> ' scriptname]);
system(['echo "export RANGE=\"-R' range_string '\"" >> ' scriptname]);
system(['echo "gmtset GRID_CROSS_SIZE 0 ANNOT_FONT_SIZE_PRIMARY 12
LABEL_FONT_SIZE 12 HEADER_FONT_SIZE 12" >> ' scriptname]);
system(['echo "gmtset PAPER_MEDIA A4" >> ' scriptname]);
system(['echo "gmtset DOTS_PR_INCH  600" >> ' scriptname]);


system(['echo "makecpt -Cno_green -T\"'  num2str(minval) '/'
num2str(maxval) '/' num2str(color_step) '\" > mycpt.cpt" >> '
scriptname]);


system(['echo "export OUTFILE="' name '.ps"" >> ' scriptname]);
system(['echo "export PDFOUTFILE="' name '.pdf"" >> ' scriptname]);

system(['echo "export FILENAME="' filename '"" >> ' scriptname]);


system(["echo \'psbasemap $RANGE $PROJ -Y1.4i -X0.9i -Bf20:.\""
title_use "\": -K  > $OUTFILE\' >> " scriptname]);

system(["echo \'grdimage $RANGE $PROJ $FILENAME -Cmycpt.cpt -Bg20 -O
-K >> $OUTFILE\' >> " scriptname]);

system(["echo \'psscale -Cmycpt.cpt -D2.3i/-0.2i/4.5i/0.4ih -O -K -B"
num2str(label_step) ":\"\":/:\"" colorbar_label "\": >> $OUTFILE\' >>
" scriptname]);

if exist('file_coords', 'var') && exist(file_coords, 'file')
system(["echo \'psxy " file_coords " -fig $RANGE $PROJ -O -K -Sx0.4c
-Wblue >> $OUTFILE\' >> " scriptname]);
end

system(["echo \'pscoast $RANGE  $PROJ -O -W -Na -Dc >> $OUTFILE\' >> "
scriptname]);

system(["echo \'ps2pdf $OUTFILE $PDFOUTFILE\' >> " scriptname]); %also
create a PDF copy, to facilitate Mac viewing

%system(["echo \'rm $OUTFILE\' >> " scriptname]); %remove the PostScript version


system(['chmod 744 ' scriptname]);
system([scriptname]);
return



endfunction