Function GMT_grdsample

src/grdsample.c:186–415

#define Return(code) {Free_Ctrl (GMT, Ctrl); gmt_end_module (GMT, GMT_cpy); bailout (code);}

EXTERN_MSC int GMT_grdsample (void *V_API, int mode, void *args) {

	int error = 0;
	openmp_int row, col;
	unsigned int registration;

	uint64_t ij;

	char format[GMT_BUFSIZ];

	double *lon = NULL, lat, wesn_i[4], wesn_o[4], inc[2], ince[2];

	struct GRDSAMPLE_CTRL *Ctrl = NULL;
	struct GMT_GRID *Gin = NULL, *Gout = NULL;
	struct GMT_GRID_HEADER_HIDDEN *HH = NULL;
	struct GMT_CTRL *GMT = NULL, *GMT_cpy = NULL;
	struct GMT_OPTION *options = NULL;
	struct GMTAPI_CTRL *API = gmt_get_api_ptr (V_API);	/* Cast from void to GMTAPI_CTRL pointer */

	/*----------------------- Standard module initialization and parsing ----------------------*/

	if (API == NULL) return (GMT_NOT_A_SESSION);
	if (mode == GMT_MODULE_PURPOSE) return (usage (API, GMT_MODULE_PURPOSE));	/* Return the purpose of program */
	options = GMT_Create_Options (API, mode, args);	if (API->error) return (API->error);	/* Set or get option list */

	if ((error = gmt_report_usage (API, options, 0, usage)) != GMT_NOERROR) bailout (error);	/* Give usage if requested */

	/* Parse the command-line arguments */

	if ((GMT = gmt_init_module (API, THIS_MODULE_LIB, THIS_MODULE_CLASSIC_NAME, THIS_MODULE_KEYS, THIS_MODULE_NEEDS, module_kw, &options, &GMT_cpy)) == NULL) bailout (API->error); /* Save current state */
	if (GMT_Parse_Common (API, THIS_MODULE_OPTIONS, options)) Return (API->error);
	Ctrl = New_Ctrl (GMT);	/* Allocate and initialize a new control structure */
	if ((error = parse (GMT, Ctrl, options)) != 0) Return (error);

	/*---------------------------- This is the grdsample main code ----------------------------*/

	gmt_grd_set_datapadding (GMT, true);	/* Turn on gridpadding when reading a subset */

	gmt_enable_threads (GMT);	/* Set number of active threads, if supported */
	GMT_Report (API, GMT_MSG_INFORMATION, "Processing input grid\n");
	gmt_set_pad (GMT, 2U);	/* Ensure space for BCs in case an API passed pad == 0 */
	if ((Gin = GMT_Read_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_CONTAINER_ONLY, NULL, Ctrl->In.file, NULL)) == NULL) {	/* Get header only */
		Return (API->error);
	}

	/* There are two separate regions that we need to worry about and ensure they are consistent when a user
	 * selects a subset. This is because the grid region and the desired region may not only be different, but
	 * may be phase-shifted. That means a single w/e/s/n cannot be used both to specify a subset for READING and
	 * also as an OUTPUT region.  Thus, below we create two regions:
	 * wesn_i: This is initially the input grid's domain.  If no -R is given then that is what we will read.
	 * wesn_o: This is the region given by -R.  If no -R then it is the same as wesn_i.
	 * If no -R is given then wesn_i == wesn_o and (presumably) there is only differences in inc and registration.
	 * If there is a -R, then we adjust wesn_i and wesn_o thus:
	 *   wesn_i: We move the boundaries inwards by the grid's own increments until the bounds are equal to or inside the -R.
	 *			 Then, if we are inside and we are able to move one step outwards we do so to ensure we cover the output range.
	 *   wesn_o: We move the boundaries inwards by the -Iinc spacing as long as we are outside the input grid.
	 */