[API] Output from grdmath cannot be read when GMT_IS_REFERENCE is used

Q&A
1

Hello,

As far as I know, and to the extend of the reference doc, all modules accept that their output be stored in a virtual file opened with the GMT_IS_REFERENCE flag. However, using this flag with grdmath makes it impossible to read the ouput. Following is a simple test program:


#include 
#include 

int main () {
  void *API;
  struct GMT_GRID *Gin = NULL;
  struct GMT_GRID *Gout = NULL;
  char input[GMT_VF_LEN] = {""};
  char output[GMT_VF_LEN] = {""};
  char command[2*GMT_VF_LEN + 10];
  double wesn[4] = {-10., 0., 40., 50.}, inc[2] = {1./60., 1./60.};
  
  API = GMT_Create_Session("GMT_grdmath", 2, 0, NULL);
  Gin = GMT_Create_Data(API, GMT_IS_GRID, GMT_IS_SURFACE, GMT_CONTAINER_AND_DATA | GMT_GRID_XY, NULL, wesn, inc, GMT_GRID_NODE_REG, -1, NULL);
  GMT_Open_VirtualFile(API, GMT_IS_GRID, GMT_IS_SURFACE, GMT_IN|GMT_IS_REFERENCE, Gin, input);
  GMT_Open_VirtualFile(API, GMT_IS_GRID, GMT_IS_SURFACE, GMT_OUT|GMT_IS_REFERENCE, NULL, output);
  //Without GMT_IS_REFERENCE on output, test passes
  //GMT_Open_VirtualFile(API, GMT_IS_GRID, GMT_IS_SURFACE, GMT_OUT, NULL, output);
  sprintf(command, "%s SQR = %s", input, output);
  GMT_Call_Module(API, "grdmath", GMT_MODULE_CMD, command);
  printf("close virtual input file\n");
  GMT_Close_VirtualFile(API, input);
  printf("destroy input grid\n");
  GMT_Destroy_Data(API, &Gin);
  printf("read output grid\n");
  Gout = GMT_Read_VirtualFile(API, output);
  printf("close virtual output file\n");
  GMT_Close_VirtualFile(API, output);
  printf("destroy output grid\n");
  GMT_Destroy_Data(API, &Gout);
  GMT_Destroy_Session(API);
  exit(0);
  
}

producing this ouput:


close virtual input file
destroy input grid
read output grid
[Session GMT_grdmath (0)]: Error returned from GMT API: GMT_OBJECT_NOT_FOUND (60)
close virtual output file
destroy output grid

When removing the GMT_IS_REFERENCE flag on the ouput virtual file opening, the output grid can be read without error. This is probably not strictly a bug, but perhaps the doc misses the point? I use the API in many programs, and it is disturbing not to be able to use the same syntax everywhere without understanding the reasons behind.

Thanks for your help,

Olivier