File: grdcut.c

package info (click to toggle)
gmt 3.4-3
  • links: PTS
  • area: main
  • in suites: woody
  • size: 3,528 kB
  • ctags: 3,140
  • sloc: ansic: 54,081; sh: 2,552; makefile: 404; asm: 38
file content (164 lines) | stat: -rw-r--r-- 5,335 bytes parent folder | download
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
/*--------------------------------------------------------------------
 *	$Id: grdcut.c,v 1.3 2001/03/19 18:15:09 pwessel Exp $
 *
 *	Copyright (c) 1991-2001 by P. Wessel and W. H. F. Smith
 *	See COPYING file for copying and redistribution conditions.
 *
 *	This program is free software; you can redistribute it and/or modify
 *	it under the terms of the GNU General Public License as published by
 *	the Free Software Foundation; version 2 of the License.
 *
 *	This program is distributed in the hope that it will be useful,
 *	but WITHOUT ANY WARRANTY; without even the implied warranty of
 *	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *	GNU General Public License for more details.
 *
 *	Contact info: gmt.soest.hawaii.edu
 *--------------------------------------------------------------------*/
/*
 * grdcut.c reads a grdfile and writes a portion within it
 * to a new file.
 *
 * Author:	Walter Smith
 * Date:	5 august, 1988
 * Version:	3.4
 *
 */
 
#include "gmt.h"

float	*grd;

main (int argc, char **argv)
{
	BOOLEAN	error = FALSE, global = FALSE;

	int	i, nx_old, ny_old, nx_new, ny_new, one_or_zero;

	double	w_new = 0.0, e_new = 0.0, s_new = 0.0, n_new = 0.0;
	double	w_old, e_old, s_old, n_old;

	char *grd_in, *grd_out, format[BUFSIZ];

	struct GRD_HEADER header, test_header;

	argc = GMT_begin (argc, argv);
	
	grd_in = grd_out = CNULL;
	
	/* Check and interpret the command line arguments */
	
	for (i =1; i < argc; i++) {
		if (argv[i][0] == '-') {
			switch(argv[i][1]) {
				/* Common parameters */
			
				case 'R':
				case 'V':
				case '\0':
					error += GMT_get_common_args (argv[i], &w_new, &e_new, &s_new, &n_new);
					break;
					
	 			case 'G':
	 				grd_out = &argv[i][2];
					break;
				default:		/* Options not recognized */
					error = TRUE;
					GMT_default_error (argv[i][1]);
					break;
			}
		}
		else
	 		grd_in = argv[i];
	}
	
	if (argc == 1 || GMT_quick) {
		fprintf (stderr,"grdcut %s - Extract subsets from grdfiles\n\n", GMT_VERSION);
		fprintf (stderr, "usage: grdcut <input_grd> -G<output_grd> -R<west/east/south/north> [-V]\n");
		
		if (GMT_quick) exit (EXIT_FAILURE);
		
		fprintf (stderr, "\t<input_grd> is file to extract a subset from.\n");
		fprintf (stderr, "\t-G specifies output grdfile\n");
		GMT_explain_option ('R');
		fprintf (stderr, "\t   Obviously, the WESN you specify must be within the WESN of the input file.\n");
		fprintf (stderr, "\t   If in doubt, run grdinfo first and check range of old file.\n");
		fprintf (stderr, "\n\tOPTIONS:\n");
		GMT_explain_option ('V');
		exit (EXIT_FAILURE);
	}

	/* Check that the options selected make sense */
	
	if (!project_info.region_supplied) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR:  Must specify -R option\n", GMT_program);
		error++;
	}
	if (!grd_out) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR -G option:  Must specify output file\n", GMT_program);
		error++;
	}
	if (!grd_in) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR:  Must specify input file\n", GMT_program);
		error++;
	}
	if (error) exit (EXIT_FAILURE);

	GMT_put_history (argc, argv);	/* Update .gmtcommands */


	if (GMT_read_grd_info (grd_in, &header)) {
		fprintf (stderr, "%s: Error opening file %s\n", GMT_program, grd_in);
		exit (EXIT_FAILURE);
	}
		
	if (s_new < header.y_min || s_new > header.y_max) error = TRUE;
	if (n_new < header.y_min || n_new > header.y_max) error = TRUE;
	global = (fabs (header.x_max - header.x_min) == 360.0);
		
	if ( !global && ((w_new < header.x_min) || (e_new > header.x_max)) ) error = TRUE;
	if (error) {
		fprintf (stderr, "%s: Subset exceeds data domain!\n", GMT_program);
		exit (EXIT_FAILURE);
	}

	/* Make sure output grid is kosher */

	test_header.x_min = w_new;	test_header.x_max = e_new;	test_header.x_inc = header.x_inc;
	test_header.y_min = s_new;	test_header.y_max = n_new;	test_header.y_inc = header.y_inc;
	GMT_grd_RI_verify (&test_header, 1);

	GMT_grd_init (&header, argc, argv, TRUE);

	w_old = header.x_min;	e_old = header.x_max;
	s_old = header.y_min;	n_old = header.y_max;
	nx_old = header.nx;	ny_old = header.ny;
	one_or_zero = (header.node_offset) ? 0 : 1;
	nx_new = irint ((e_new - w_new) / header.x_inc) + one_or_zero;
	ny_new = irint ((n_new - s_new) / header.y_inc) + one_or_zero;
			
	grd = (float *) GMT_memory (VNULL, (size_t)(nx_new * ny_new), sizeof (float), GMT_program);
	if (GMT_read_grd (grd_in, &header, grd, w_new, e_new, s_new, n_new, GMT_pad, FALSE)) {
		fprintf (stderr, "%s: Error reading file %s\n", GMT_program, grd_in);
		exit (EXIT_FAILURE);
	}

	if (gmtdefs.verbose) {
		sprintf (format, "\t%s\t%s\t%s\t%s\t%s\t%s\t%%d\t%%d\n\0", gmtdefs.d_format, gmtdefs.d_format,
			gmtdefs.d_format, gmtdefs.d_format, gmtdefs.d_format, gmtdefs.d_format);
		fprintf (stderr, "%s: File spec:\tW E S N dx dy nx ny:\n", GMT_program);
		fprintf (stderr, "%s: Old:", GMT_program);
		fprintf (stderr, format, w_old, e_old, s_old, n_old, header.x_inc, header.y_inc, nx_old, ny_old);
		fprintf (stderr, "%s: New:", GMT_program);
		fprintf (stderr, format, w_new, e_new, s_new, n_new, header.x_inc, header.y_inc, nx_new, ny_new);
	}

	if (GMT_write_grd (grd_out, &header, grd, 0.0, 0.0, 0.0, 0.0, GMT_pad, FALSE)) {
		fprintf (stderr, "%s: Error writing file %s\n", GMT_program, grd_out);
		exit (EXIT_FAILURE);
	}
	
	GMT_free ((void *) grd);
	
	GMT_end (argc, argv);
}