192 lines
6.8 KiB
C++
192 lines
6.8 KiB
C++
/********************************************************
|
|
* ██████╗ ██████╗████████╗██╗
|
|
* ██╔════╝ ██╔════╝╚══██╔══╝██║
|
|
* ██║ ███╗██║ ██║ ██║
|
|
* ██║ ██║██║ ██║ ██║
|
|
* ╚██████╔╝╚██████╗ ██║ ███████╗
|
|
* ╚═════╝ ╚═════╝ ╚═╝ ╚══════╝
|
|
* Geophysical Computational Tools & Library (GCTL)
|
|
*
|
|
* Copyright (c) 2022 Yi Zhang (yizhang-geo@zju.edu.cn)
|
|
*
|
|
* GCTL is distributed under a dual licensing scheme. You can redistribute
|
|
* it and/or modify it under the terms of the GNU Lesser General Public
|
|
* License as published by the Free Software Foundation, either version 2
|
|
* of the License, or (at your option) any later version. You should have
|
|
* received a copy of the GNU Lesser General Public License along with this
|
|
* program. If not, see <http://www.gnu.org/licenses/>.
|
|
*
|
|
* If the terms and conditions of the LGPL v.2. would prevent you from using
|
|
* the GCTL, please consider the option to obtain a commercial license for a
|
|
* fee. These licenses are offered by the GCTL's original author. As a rule,
|
|
* licenses are provided "as-is", unlimited in time for a one time fee. Please
|
|
* send corresponding requests to: yizhang-geo@zju.edu.cn. Please do not forget
|
|
* to include some description of your company and the realm of its activities.
|
|
* Also add information on how to contact you by electronic and paper mail.
|
|
******************************************************/
|
|
|
|
#include "gmt_JX_single.h"
|
|
|
|
#ifdef GCTL_GRAPHIC_GMT
|
|
|
|
gctl::gmt_JX_single::gmt_JX_single()
|
|
{
|
|
name.resize(5);
|
|
cmd.resize(5);
|
|
|
|
name[0] = "gmtset";
|
|
cmd[0] = "FONT_ANNOT_PRIMARY=10.5p,Times-Roman,black \
|
|
MAP_FRAME_PEN=thinnest,black \
|
|
MAP_GRID_PEN_PRIMARY=thinnest,black \
|
|
MAP_TICK_PEN_PRIMARY=thinnest,black \
|
|
MAP_TICK_LENGTH_PRIMARY=1p/0.5p \
|
|
MAP_TITLE_OFFSET=7.5p \
|
|
MAP_GRID_CROSS_SIZE_PRIMARY=2p \
|
|
FONT_LABEL=10.5p,Times-Roman,black \
|
|
MAP_LABEL_OFFSET=5p \
|
|
MAP_ANNOT_OFFSET_PRIMARY=2.5p";
|
|
name[1] = "grd2cpt";
|
|
cmd[1] = "-Crainbow -Z -D";
|
|
name[2] = "grdimage";
|
|
cmd[2] = "-Bxag+l\"x (m)\" -Byag+l\"y (m)\"";
|
|
name[3] = "psscale";
|
|
cmd[3] = "-Bxa -By+lm";
|
|
name[4] = "psconvert";
|
|
cmd[4] = "-A -TEG -E300";
|
|
}
|
|
|
|
void gctl::gmt_JX_single::set_command(std::string in_name, std::string in_cmd)
|
|
{
|
|
for (int i = 0; i < name.size(); ++i)
|
|
{
|
|
if (name[i] == in_name)
|
|
{
|
|
cmd[i] = in_cmd;
|
|
return;
|
|
}
|
|
}
|
|
return;
|
|
}
|
|
|
|
void gctl::gmt_JX_single::show_command(std::ostream &out)
|
|
{
|
|
for (int i = 0; i < name.size(); ++i)
|
|
{
|
|
out << "Module: " << name[i] << " Command: " << cmd[i] << std::endl;
|
|
}
|
|
return;
|
|
}
|
|
|
|
std::string gctl::gmt_JX_single::get_command(std::string in_name)
|
|
{
|
|
for (int i = 0; i < name.size(); ++i)
|
|
{
|
|
if (name[i] == in_name)
|
|
{
|
|
return cmd[i];
|
|
}
|
|
}
|
|
|
|
throw runtime_error("No command found for module: "+in_name+". From gmt_JX_single::get_command(...)");
|
|
}
|
|
|
|
void gctl::gmt_JX_single::plot(std::string plot_name, const array<double> &data, double xmin, double xmax, double ymin, double ymax, int xnum, int ynum)
|
|
{
|
|
std::string cpt_file = plot_name+".cpt";
|
|
std::string ps_file = plot_name+".ps";
|
|
|
|
// Initiate a new GMT session
|
|
// you need to destroy the session later
|
|
void *API = GMT_Create_Session("gctl_gmt_plot", 2U, 0, NULL);
|
|
|
|
// prepare header info
|
|
int m = xnum, n = ynum;
|
|
double wesn[4] = {xmin, xmax, ymin, ymax};
|
|
double inc[2] = {(xmax-xmin)/(xnum-1), (ymax-ymin)/(ynum-1)};
|
|
// create an empty data container of GMT_GRID type
|
|
// forcedly convert pointer type (C++ standard)
|
|
// you need to destroy the data later
|
|
struct GMT_GRID *G = reinterpret_cast<GMT_GRID*>(
|
|
GMT_Create_Data(API, GMT_IS_GRID, GMT_IS_SURFACE,
|
|
GMT_CONTAINER_AND_DATA, NULL, &wesn[0], &inc[0],
|
|
GMT_GRID_NODE_REG, -1, NULL));
|
|
|
|
// manipulate grid data
|
|
double *x_coord = GMT_Get_Coord(API, GMT_IS_GRID, GMT_X, G);
|
|
double *y_coord = GMT_Get_Coord(API, GMT_IS_GRID, GMT_Y, G);
|
|
|
|
int idx;
|
|
double minz = 1e+30, maxz = -1e+30;
|
|
for (int i = 0; i < G->header->n_rows; i++)
|
|
{
|
|
for (int j = 0; j < G->header->n_columns; j++)
|
|
{
|
|
idx = GMT_Get_Index(API, G->header, i, j);
|
|
G->data[idx] = data[j+(ynum-1-i)*xnum];
|
|
|
|
if (G->data[idx] < minz) minz = G->data[idx];
|
|
if (G->data[idx] > maxz) maxz = G->data[idx];
|
|
}
|
|
}
|
|
|
|
G->header->z_min = minz;
|
|
G->header->z_max = maxz;
|
|
|
|
// start the plotting
|
|
// 1. set GMT defaults
|
|
std::string args_defaults = get_command("gmtset");
|
|
// call gmtset
|
|
GMT_Call_Module (API, "gmtset", GMT_MODULE_CMD, (char*) args_defaults.c_str());
|
|
|
|
// load the grid to a virtual file for plotting
|
|
char grid_name[GMT_VF_LEN] = {""};
|
|
GMT_Open_VirtualFile (API, GMT_IS_GRID, GMT_IS_SURFACE, GMT_IN, G, grid_name);
|
|
// get string type of grid name
|
|
std::string grid_name_str = grid_name;
|
|
|
|
// prepare cpt file
|
|
std::string range_str = std::to_string(xmin)+"/"+std::to_string(xmax)+"/"+std::to_string(ymin)+"/"+std::to_string(ymax);
|
|
std::string args_cpt = grid_name_str + " " + get_command("grd2cpt") + " -R" + range_str + " ->" + cpt_file;
|
|
// call grd2cpt
|
|
GMT_Call_Module (API, "grd2cpt", GMT_MODULE_CMD, (char*) args_cpt.c_str());
|
|
|
|
// plot the image
|
|
std::string args_image = grid_name_str + " -R" + range_str + " -C" + cpt_file + " -JX1.5i/1.5i -X0.5i -Y0.5i --MAP_FRAME_AXES=WesNZ -K -P " + get_command("grdimage")
|
|
+ " ->" + ps_file;
|
|
// call grdimage
|
|
GMT_Call_Module (API, "grdimage", GMT_MODULE_CMD, (char*) args_image.c_str());
|
|
|
|
// plot color bar
|
|
std::string args_bar = get_command("psscale") + " -Dx0.1i/-0.2i+w1.3i/0.05i+h -C" + cpt_file + " -O -)" + ps_file;
|
|
// call psscale
|
|
GMT_Call_Module (API, "psscale", GMT_MODULE_CMD, (char*) args_bar.c_str());
|
|
|
|
// convert ps file to raster formats
|
|
std::string args_pic = ps_file + " " + get_command("psconvert");
|
|
// call psconvert
|
|
GMT_Call_Module (API, "psconvert", GMT_MODULE_CMD, (char*) args_pic.c_str());
|
|
|
|
// close virtual file
|
|
GMT_Close_VirtualFile (API, grid_name);
|
|
|
|
// destroy data container
|
|
GMT_Destroy_Data(API, G);
|
|
|
|
// end the GMT session
|
|
GMT_Destroy_Session(API);
|
|
|
|
// remove temporary files. including GMT defaults
|
|
remove(cpt_file.c_str());
|
|
remove(ps_file.c_str());
|
|
// if exist, delete
|
|
std::string hist_file = "gmt.history";
|
|
std::string conf_file = "gmt.conf";
|
|
if (!access(hist_file.c_str(), F_OK))
|
|
remove(hist_file.c_str());
|
|
if (!access(conf_file.c_str(), F_OK))
|
|
remove(conf_file.c_str());
|
|
|
|
return;
|
|
}
|
|
|
|
#endif // GCTL_GRAPHIC_GMT
|