Address issues raised in #8438 and more. (#8440)

src/gmt_init.c

@@ -15794,9 +15794,14 @@ struct GMT_CTRL *gmt_init_module (struct GMTAPI_CTRL *API, const char *lib_name,
  GMT->current.ps.active = is_PS; /* true if module will produce PS */
 
  /* Check if there is an input remote grid or memory file so we may set geographic to be true now before we must decide on auto-J */
- if (options && (opt = GMT_Find_Option (API, GMT_OPT_INFILE, *options))) {
- if (gmt_remote_dataset_id (API, opt->arg) != GMT_NOTSET) /* All remote data grids/images are geographic */
- gmt_set_geographic (GMT, GMT_IN);
+ if (options && (opt = GMT_Find_Option(API, GMT_OPT_INFILE, *options))) {
+ if (gmt_remote_dataset_id(API, opt->arg) != GMT_NOTSET) { /* All remote data grids/images are geographic */
+ gmt_set_geographic(GMT, GMT_IN);
+ /* Above is not enough for tilled grids. That case looses this info when creating the final grid but a not
+ consulting the units in the header of any of the tiles. Se we reinforce it with a fake -fg option.
+ */
+ API->GMT->common.f.is_geo[0] = true;
+ }
  else if (gmtlib_data_is_geographic (API, opt->arg)) /* This dataset, grid, image, matrix, or vector is geographic */
  gmt_set_geographic (GMT, GMT_IN);
  }
@@ -18238,7 +18243,7 @@ unsigned int gmt_parse_inc_option (struct GMT_CTRL *GMT, char option, char *item
 
 GMT_LOCAL int gmtinit_parse_proj4 (struct GMT_CTRL *GMT, char *item, char *dest) {
  /* Deal with proj.4 or EPSGs passed in -J option */
- char *item_t1 = NULL, *item_t2 = NULL, wktext[32] = {""}, *pch;
+ char *item_t1 = NULL, *item_t2 = NULL, epsg2proj[GMT_LEN256] = {""}, wktext[32] = {""}, *pch;
  bool do_free = false;
  int error = 0, scale_pos;
  size_t k, len;
@@ -18270,7 +18275,7 @@ GMT_LOCAL int gmtinit_parse_proj4 (struct GMT_CTRL *GMT, char *item, char *dest)
  /* Don't remember anymore if we still need to call gmt_importproj4(). We don't for simple
  mapproject usage but maybe we still need for mapping purposes. To-be-rediscovered.
  */
- item_t2 = gmt_importproj4 (GMT, item_t1, &scale_pos); /* This is GMT -J proj string */
+ item_t2 = gmt_importproj4 (GMT, item_t1, &scale_pos, epsg2proj); /* This is the GMT -J proj string */
  if (item_t2 && !GMT->current.ps.active && !strcmp(item_t2, "/1:1")) {
  /* Even though it failed to do the mapping we can still use it in mapproject */
  GMT->current.proj.projection_GMT = GMT_NO_PROJ;
@@ -18290,6 +18295,10 @@ GMT_LOCAL int gmtinit_parse_proj4 (struct GMT_CTRL *GMT, char *item, char *dest)
  }
  else /* Not particularly useful yet because mapproject will fail anyway when scale != 1:1 */
  GMT->current.proj.projection_GMT = GMT_NO_PROJ;
+
+ /* If input was a EPSG code, save the epsg -> proj conversion string. Save also if input was a +proj string */
+ if (epsg2proj[0] != '\0') snprintf(GMT->common.J.proj4string, GMT_LEN256-1, "%s", epsg2proj);
+ else if (item_t1[0] == '+') snprintf(GMT->common.J.proj4string, GMT_LEN256-1, "%s", item_t1);
 
  /* Check if the scale is 1 or 1:1, and don't get fooled with, for example, 1:10 */
  pch = &item_t2[scale_pos];

src/gmt_plot.c

@@ -8014,9 +8014,11 @@ GMT_LOCAL void gmtplot_wipe_substr(char *str1, char *str2) {
 }
 #endif
 
-char *gmt_importproj4 (struct GMT_CTRL *GMT, char *pStr, int *scale_pos) {
+char *gmt_importproj4 (struct GMT_CTRL *GMT, char *pStr, int *scale_pos, char *epsg2proj) {
  /* Take a PROJ.4 projection string or EPSG code and try to find the equivalent -J syntax
  scale_pos is position on the return string where starts the scale sub-string.
+ The pStr pointer is changed when in input it contains an EPSG code. On return it will
+ have the proj4 string corresponding to the EPSG code. That pointer can later be freed with free()
  */
  unsigned int pos = 0;
  //bool got_lonlat = false;
@@ -8031,7 +8033,7 @@ char *gmt_importproj4 (struct GMT_CTRL *GMT, char *pStr, int *scale_pos) {
  pch[0] = '\0';
  }
  else {
- GMT->current.ps.active ? sprintf(scale_c, "14c") : sprintf(scale_c, "1:1");
+ GMT->current.ps.active ? sprintf(scale_c, "15c") : sprintf(scale_c, "1:1");
  }
 
  if (isdigit(szProj4[1])) { /* A EPSG code. By looking at 2nd char instead of 1st both +epsg and epsg work */
@@ -8066,6 +8068,7 @@ char *gmt_importproj4 (struct GMT_CTRL *GMT, char *pStr, int *scale_pos) {
  return (pStrOut);
  }
  snprintf(szProj4, GMT_LEN256-1, "%s", pszResult);
+ snprintf(epsg2proj, GMT_LEN256-1, "%s", pszResult); /* This one now has the proj4 conversion of the EPSG code. */
  if (scale_c[0] != '\0') strcat (szProj4, scale_c); /* Add the width/scale found above */
  CPLFree(pszResult);
  OSRDestroySpatialReference(hSRS);

src/gmt_prototypes.h

@@ -280,8 +280,8 @@ EXTERN_MSC void gmt_map_text (struct GMT_CTRL *GMT, double x, double y, struct G
 EXTERN_MSC void gmt_map_title (struct GMT_CTRL *GMT, double x, double y);
 EXTERN_MSC void gmt_linearx_grid (struct GMT_CTRL *GMT, struct PSL_CTRL *P, double w, double e, double s, double n, double dval);
 EXTERN_MSC int gmt_ps_append (struct GMT_CTRL *GMT, char *source, unsigned int mode, FILE *dest);
-EXTERN_MSC char * gmt_export2proj4 (struct GMT_CTRL *GMT);
-EXTERN_MSC char * gmt_importproj4 (struct GMT_CTRL *GMT, char *szProj4, int *scale_pos);
+EXTERN_MSC char *gmt_export2proj4 (struct GMT_CTRL *GMT);
+EXTERN_MSC char *gmt_importproj4 (struct GMT_CTRL *GMT, char *szProj4, int *scale_pos, char *epsg2proj);
 EXTERN_MSC int gmt_strip_layer (struct GMTAPI_CTRL *API, int nlayers);
 EXTERN_MSC void gmt_textpath_init (struct GMT_CTRL *GMT, struct GMT_PEN *BP, double Brgb[]);
 EXTERN_MSC void gmt_draw_map_rose (struct GMT_CTRL *GMT, struct GMT_MAP_ROSE *mr);

src/grdblend.c

@@ -1084,8 +1084,8 @@ EXTERN_MSC int GMT_grdblend (void *V_API, int mode, void *args) {
  n_fill++; /* One more cell filled */
  if (z[col] < Grid->header->z_min) Grid->header->z_min = z[col]; /* Update the extrema for output grid */
  if (z[col] > Grid->header->z_max) Grid->header->z_max = z[col];
- if (gmt_M_is_zero (z[col]))
- m = 1;
+ if (gmt_M_is_zero (z[col]))
+ m = 1;
  }
  else /* No grids covered this node, defaults to the no_data value */
  z[col] = no_data_f;
@@ -1162,5 +1162,11 @@ EXTERN_MSC int GMT_grdblend (void *V_API, int mode, void *args) {
  GMT_Report (GMT->parent, GMT_MSG_ERROR, "Failed to delete file %s\n", outfile);
  }
 
+ /* These two were not copied back into GMT and would be lost in gmt_end_module. Needed to set the correct
+ @PROJ parameters in PS files. See issue #8438
+ */
+ GMT_cpy->current.io.col_type[0][0] = GMT->current.io.col_type[0][0];
+ GMT_cpy->current.io.col_type[0][1] = GMT->current.io.col_type[0][1];
+
  Return (GMT_NOERROR);
 }

src/grdimage.c

@@ -1511,11 +1511,16 @@ EXTERN_MSC int GMT_grdimage (void *V_API, int mode, void *args) {
 
  /* Determine if grid/image is to be projected */
  need_to_project = (gmt_M_is_nonlinear_graticule (GMT) || Ctrl->E.dpi > 0);
- if (need_to_project)
- GMT_Report (API, GMT_MSG_DEBUG, "Projected grid is non-orthogonal, nonlinear, or dpi was changed\n");
- else if (Ctrl->D.active) /* If not projecting no need for a pad */
- gmt_set_pad(GMT, 0);
+ if (need_to_project && GMT->current.proj.projection == GMT_PROJ4_PROJS) {
+ if (GMT->current.proj.is_proj4)
+ if (strstr(GMT->common.J.proj4string, "longlat") || strstr(GMT->common.J.proj4string, "lonlat") || strstr(GMT->common.J.proj4string, "latlon"))
+ need_to_project = false;
+ }
+
+ if (Ctrl->D.active) gmt_set_pad(GMT, 0); /* If not projecting no need for a pad */
+
  pad_mode = (need_to_project) ? GMT_GRID_NEEDS_PAD2 : 0;
+
  /* Read the illumination grid header right away so we can use its region to set that of an image (if requested) */
  if (use_intensity_grid) { /* Illumination grid desired */
  if (Ctrl->I.derive) { /* Illumination grid must be derived */
@@ -1689,16 +1694,10 @@ EXTERN_MSC int GMT_grdimage (void *V_API, int mode, void *args) {
  Return (API->error);
  }
  if (!GMT->common.J.active) { /* No map projection was supplied, set default */
- if ((Grid_orig->header->ProjRefWKT != NULL) || (Grid_orig->header->ProjRefPROJ4 != NULL)) {
- static char *Jarg[2] = {"X15c", "Q15c+"}; /* The two default projections for Cartesian and Geographic */
- unsigned int kind = gmt_M_is_geographic (GMT, GMT_IN);
- gmt_parse_common_options (GMT, "J", 'J', Jarg[kind]); /* No projection specified, use linear or equidistant */
- GMT->common.J.active = true; /* Now parsed */
- }
- else if (GMT->current.setting.run_mode == GMT_CLASSIC) {
- GMT_Report (API, GMT_MSG_ERROR, "Must specify a map projection with the -J option\n");
- Return (GMT_PARSE_ERROR);
- }
+ static char *Jarg[2] = {"X15c/0", "Q15c+"}; /* The two default projections for Cartesian and Geographic */
+ unsigned int kind = gmt_M_is_geographic (GMT, GMT_IN);
+ gmt_parse_common_options (GMT, "J", 'J', Jarg[kind]); /* No projection specified, use linear or equidistant */
+ GMT->common.J.active = true; /* Now parsed */
  }
  }
  if (!Ctrl->C.active) /* Set no specific CPT so we turn -C on to use current or default CPT */
@@ -1718,11 +1717,13 @@ EXTERN_MSC int GMT_grdimage (void *V_API, int mode, void *args) {
 
  if (Ctrl->E.dpi == 0) grdimage_adjust_R_consideration (GMT, header_work); /* Special check for global pixel-registered plots */
 
- /* Initialize the projection for the selected -R -J */
- if (gmt_map_setup (GMT, GMT->common.R.wesn)) Return (GMT_PROJECTION_ERROR);
+ /* Initialize the projection for the selected -R -J. It shouldn't be necessary when -A is in effect because it
+ doesn't produce a map, but the grid projection machinery uses variables set by it. Can go away only when we implement
+ image projections with GDAL instead. 
+ */
+ if (gmt_map_setup(GMT, GMT->common.R.wesn)) Return (GMT_PROJECTION_ERROR);
 
  /* Determine the wesn to be used to read the grid file; or bail if file is outside -R */
-
  if (!gmt_grd_setregion (GMT, header_work, wesn, need_to_project * GMT->common.n.interpolant))
  nothing_inside = true; /* Means we can return an empty plot, possibly with a basemap */
  else if (use_intensity_grid && !Ctrl->I.derive && !gmt_grd_setregion (GMT, Intens_orig->header, wesn, need_to_project * GMT->common.n.interpolant))
@@ -1997,12 +1998,16 @@ tr_image: GMT_Report (API, GMT_MSG_INFORMATION, "Project the input image\n");
  GMT_Report (API, GMT_MSG_INFORMATION, "Project the input grid\n");
  if ((Grid_proj = GMT_Duplicate_Data (API, GMT_IS_GRID, GMT_DUPLICATE_NONE, Grid_orig)) == NULL)
  Return (API->error); /* Just to get a header we can change */
+
  /* Determine the dimensions of the projected grid */
+ if (Ctrl->A.active && GMT->current.proj.rect[0] == GMT->current.proj.rect[1]) /* Somewhere behind we missed doing this */
+ gmt_M_memcpy(GMT->current.proj.rect, Grid_orig->header->wesn, 4, double);
  grdimage_set_proj_limits (GMT, Grid_proj->header, Grid_orig->header, need_to_project, mixed);
+
  if (grid_registration == GMT_GRID_NODE_REG) /* Force pixel if a dpi was specified, else keep as is */
  grid_registration = (Ctrl->E.dpi > 0) ? GMT_GRID_PIXEL_REG : Grid_orig->header->registration;
  if (gmt_M_err_fail (GMT, gmt_project_init (GMT, Grid_proj->header, inc, nx_proj, ny_proj, Ctrl->E.dpi, grid_registration),
- Ctrl->In.file)) Return (GMT_PROJECTION_ERROR);
+  Ctrl->In.file)) Return (GMT_PROJECTION_ERROR);
  gmt_set_grddim (GMT, Grid_proj->header); /* Recalculate projected grid dimensions */
  if (GMT_Create_Data (API, GMT_IS_GRID, GMT_IS_SURFACE, GMT_DATA_ONLY, NULL, NULL, NULL, 0, 0, Grid_proj) == NULL)
  Return (API->error); /* Failed to allocate memory for the projected grid */
@@ -2077,6 +2082,10 @@ tr_image: GMT_Report (API, GMT_MSG_INFORMATION, "Project the input image\n");
  else /* Use a different reference header for the image */
  header_work = Img_proj->header; /* Later when need to refer to the header, use this copy */
 
+ if (header_work->ProjRefPROJ4 == NULL && GMT->current.proj.is_proj4) header_work->ProjRefPROJ4 = strdup(GMT->common.J.proj4string);
+ if (!need_to_project && Ctrl->A.active && header_work->ProjRefPROJ4 == NULL && gmt_M_is_geographic(GMT, GMT_IN))
+ header_work->ProjRefPROJ4 = strdup("+proj=lonlat");
+
  /* Assign the Conf structure members */
  Conf->Grid = Grid_proj;
  Conf->Image = Img_proj;

src/psconvert.c

@@ -2786,7 +2786,6 @@ EXTERN_MSC int GMT_psconvert (void *V_API, int mode, void *args) {
  else if (Ctrl->W.warp && !proj4_cmd)
  GMT_Report (API, GMT_MSG_ERROR, "Could not find the Proj4 command in the PS file. No conversion performed.\n");
  }
-
  else if (Ctrl->W.kml) { /* Write a basic kml file */
  char kml_file[PATH_MAX] = "";
  if (Ctrl->D.active)
@@ -2860,7 +2859,6 @@ EXTERN_MSC int GMT_psconvert (void *V_API, int mode, void *args) {
  GMT_Report (API, GMT_MSG_INFORMATION, "Wrote KML file %s\n", kml_file);
  }
  }
-
  else if (Ctrl->W.active && !found_proj) {
  GMT_Report (API, GMT_MSG_ERROR, "Could not find the 'PROJ' tag in the PS file. No world file created.\n");
  GMT_Report (API, GMT_MSG_ERROR, "This situation occurs when one of the two next cases is true:\n");