diff --git a/src/gmt_init.c b/src/gmt_init.c index 1895dffc57c..a5226ec184d 100644 --- a/src/gmt_init.c +++ b/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]; diff --git a/src/gmt_plot.c b/src/gmt_plot.c index d84a4180b31..b94f7db1e79 100644 --- a/src/gmt_plot.c +++ b/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); diff --git a/src/gmt_prototypes.h b/src/gmt_prototypes.h index a41cf69b10a..b512395bcbf 100644 --- a/src/gmt_prototypes.h +++ b/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); diff --git a/src/grdblend.c b/src/grdblend.c index 4f9b3d6c8f4..5f80b963c4a 100644 --- a/src/grdblend.c +++ b/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); } diff --git a/src/grdimage.c b/src/grdimage.c index 0d531779e8f..8054f03ad00 100644 --- a/src/grdimage.c +++ b/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; diff --git a/src/psconvert.c b/src/psconvert.c index d0b168675f5..da56e649c71 100644 --- a/src/psconvert.c +++ b/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");