GRASS GIS 7 Programmer's Manual  7.5.svn(2018)-r72645
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
sample.c
Go to the documentation of this file.
1 /*!
2  \file lib/raster/sample.c
3 
4  \brief Raster library - Sampling methods (extract a cell value from
5  raster map)
6 
7  1/2006: moved to libgis from v.sample/v.drape for clone removal
8 
9  (C) 2001-2009 by the GRASS Development Team
10 
11  This program is free software under the GNU General Public License
12  (>=v2). Read the file COPYING that comes with GRASS for details.
13 
14  \author James Darrell McCauley <darrell mccauley-usa.com>, http://mccauley-usa.com/
15  */
16 
17 #include <string.h>
18 #include <unistd.h>
19 #include <math.h>
20 
21 #include <grass/gis.h>
22 #include <grass/raster.h>
23 #include <grass/glocale.h>
24 
25 /* prototypes */
26 static double scancatlabel(const char *);
27 
28 /*!
29  * \brief Extract a cell value from raster map.
30  *
31  * Extract a cell value from raster map at given northing and easting
32  * with a sampled 3x3 window using a specified interpolation method.
33  *
34  * - NEAREST neighbor interpolation
35  * - BILINEAR bilinear interpolation
36  * - CUBIC cubic interpolation
37  *
38  * \param fd file descriptor
39  * \param window region settings
40  * \param cats categories
41  * \param north northing position
42  * \param east easting position
43  * \param usedesc flag to scan category label
44  * \param itype interpolation method
45  *
46  * \return cell value at given position
47  */
49  const struct Cell_head *window,
50  struct Categories *cats,
51  double north, double east,
52  int usedesc, INTERP_TYPE itype)
53 {
54  double retval;
55 
56  switch (itype) {
57  case INTERP_NEAREST:
58  retval =
59  Rast_get_sample_nearest(fd, window, cats, north, east, usedesc);
60  break;
61  case INTERP_BILINEAR:
62  retval =
63  Rast_get_sample_bilinear(fd, window, cats, north, east, usedesc);
64  break;
65  case INTERP_BICUBIC:
66  retval =
67  Rast_get_sample_cubic(fd, window, cats, north, east, usedesc);
68  break;
69  default:
70  G_fatal_error("Rast_get_sample: %s", _("Unknown interpolation type"));
71  }
72 
73  return retval;
74 }
75 
76 /*!
77  * \brief Extract a cell value from raster map (neighbor interpolation)
78  *
79  * Extract a cell value from raster map at given northing and easting
80  * with a sampled 3x3 window using a neighbor interpolation.
81  *
82  * \param fd file descriptor
83  * \param window region settings
84  * \param cats categories
85  * \param north northing position
86  * \param east easting position
87  * \param usedesc flag to scan category label
88  *
89  * \return cell value at given position
90  */
92  const struct Cell_head * window,
93  struct Categories * cats,
94  double north, double east, int usedesc)
95 {
96  int row, col;
97  DCELL result;
98  DCELL *maprow = Rast_allocate_d_buf();
99 
100  /* convert northing and easting to row and col, resp */
101  row = (int)floor(Rast_northing_to_row(north, window));
102  col = (int)floor(Rast_easting_to_col(east, window));
103 
104  if (row < 0 || row >= Rast_window_rows() ||
105  col < 0 || col >= Rast_window_cols()) {
106  Rast_set_d_null_value(&result, 1);
107  goto done;
108  }
109 
110  Rast_get_d_row(fd, maprow, row);
111 
112  if (Rast_is_d_null_value(&maprow[col])) {
113  Rast_set_d_null_value(&result, 1);
114  goto done;
115  }
116 
117  if (usedesc) {
118  char *buf = Rast_get_c_cat((CELL *) & (maprow[col]), cats);
119 
120  G_squeeze(buf);
121  result = scancatlabel(buf);
122  }
123  else
124  result = maprow[col];
125 
126  done:
127  G_free(maprow);
128 
129  return result;
130 }
131 
132 
133 /*!
134  * \brief Extract a cell value from raster map (bilinear interpolation).
135  *
136  * Extract a cell value from raster map at given northing and easting
137  * with a sampled 3x3 window using a bilinear interpolation.
138  *
139  * \param fd file descriptor
140  * \param window region settings
141  * \param cats categories
142  * \param north northing position
143  * \param east easting position
144  * \param usedesc flag to scan category label
145  *
146  * \return cell value at given position
147  */
149  const struct Cell_head * window,
150  struct Categories * cats,
151  double north, double east, int usedesc)
152 {
153  int row, col;
154  double grid[2][2];
155  DCELL *arow = Rast_allocate_d_buf();
156  DCELL *brow = Rast_allocate_d_buf();
157  double frow, fcol, trow, tcol;
158  DCELL result;
159 
160  frow = Rast_northing_to_row(north, window);
161  fcol = Rast_easting_to_col(east, window);
162 
163  /* convert northing and easting to row and col, resp */
164  row = (int)floor(frow - 0.5);
165  col = (int)floor(fcol - 0.5);
166 
167  trow = frow - row - 0.5;
168  tcol = fcol - col - 0.5;
169 
170  if (row < 0 || row + 1 >= Rast_window_rows() ||
171  col < 0 || col + 1 >= Rast_window_cols()) {
172  Rast_set_d_null_value(&result, 1);
173  goto done;
174  }
175 
176  Rast_get_d_row(fd, arow, row);
177  Rast_get_d_row(fd, brow, row + 1);
178 
179  if (Rast_is_d_null_value(&arow[col]) ||
180  Rast_is_d_null_value(&arow[col + 1]) ||
181  Rast_is_d_null_value(&brow[col]) ||
182  Rast_is_d_null_value(&brow[col + 1])) {
183  Rast_set_d_null_value(&result, 1);
184  goto done;
185  }
186 
187  /*-
188  * now were ready to do bilinear interpolation over
189  * arow[col], arow[col+1],
190  * brow[col], brow[col+1]
191  */
192 
193  if (usedesc) {
194  char *buf;
195 
196  G_squeeze(buf = Rast_get_c_cat((int *)&(arow[col]), cats));
197  grid[0][0] = scancatlabel(buf);
198  G_squeeze(buf = Rast_get_c_cat((CELL *) & (arow[col + 1]), cats));
199  grid[0][1] = scancatlabel(buf);
200  G_squeeze(buf = Rast_get_c_cat((CELL *) & (brow[col]), cats));
201  grid[1][0] = scancatlabel(buf);
202  G_squeeze(buf = Rast_get_c_cat((CELL *) & (brow[col + 1]), cats));
203  grid[1][1] = scancatlabel(buf);
204  }
205  else {
206  grid[0][0] = arow[col];
207  grid[0][1] = arow[col + 1];
208  grid[1][0] = brow[col];
209  grid[1][1] = brow[col + 1];
210  }
211 
212  result = Rast_interp_bilinear(tcol, trow,
213  grid[0][0], grid[0][1], grid[1][0],
214  grid[1][1]);
215 
216  done:
217  G_free(arow);
218  G_free(brow);
219 
220  return result;
221 }
222 
223 /*!
224  * \brief Extract a cell value from raster map (cubic interpolation).
225  *
226  * Extract a cell value from raster map at given northing and easting
227  * with a sampled 3x3 window using a cubic interpolation.
228  *
229  * \param fd file descriptor
230  * \param window region settings
231  * \param cats categories
232  * \param north northing position
233  * \param east easting position
234  * \param usedesc flag to scan category label
235  *
236  * \return cell value at given position
237  */
239  const struct Cell_head * window,
240  struct Categories * cats,
241  double north, double east, int usedesc)
242 {
243  int i, j, row, col;
244  double grid[4][4];
245  DCELL *rows[4];
246  double frow, fcol, trow, tcol;
247  DCELL result;
248 
249  for (i = 0; i < 4; i++)
250  rows[i] = Rast_allocate_d_buf();
251 
252  frow = Rast_northing_to_row(north, window);
253  fcol = Rast_easting_to_col(east, window);
254 
255  /* convert northing and easting to row and col, resp */
256  row = (int)floor(frow - 1.5);
257  col = (int)floor(fcol - 1.5);
258 
259  trow = frow - row - 1.5;
260  tcol = fcol - col - 1.5;
261 
262  if (row < 0 || row + 3 >= Rast_window_rows() ||
263  col < 0 || col + 3 >= Rast_window_cols()) {
264  Rast_set_d_null_value(&result, 1);
265  goto done;
266  }
267 
268  for (i = 0; i < 4; i++)
269  Rast_get_d_row(fd, rows[i], row + i);
270 
271  for (i = 0; i < 4; i++)
272  for (j = 0; j < 4; j++)
273  if (Rast_is_d_null_value(&rows[i][col + j])) {
274  Rast_set_d_null_value(&result, 1);
275  goto done;
276  }
277 
278  /*
279  * now were ready to do cubic interpolation over
280  * arow[col], arow[col+1], arow[col+2], arow[col+3],
281  * brow[col], brow[col+1], brow[col+2], brow[col+3],
282  * crow[col], crow[col+1], crow[col+2], crow[col+3],
283  * drow[col], drow[col+1], drow[col+2], drow[col+3],
284  */
285 
286  if (usedesc) {
287  char *buf;
288 
289  for (i = 0; i < 4; i++) {
290  for (j = 0; j < 4; j++) {
291  G_squeeze(buf =
292  Rast_get_c_cat((CELL *) & (rows[i][col + j]),
293  cats));
294  grid[i][j] = scancatlabel(buf);
295  }
296  }
297  }
298  else {
299  for (i = 0; i < 4; i++)
300  for (j = 0; j < 4; j++)
301  grid[i][j] = rows[i][col + j];
302  }
303 
304  result = Rast_interp_bicubic(tcol, trow,
305  grid[0][0], grid[0][1], grid[0][2],
306  grid[0][3], grid[1][0], grid[1][1],
307  grid[1][2], grid[1][3], grid[2][0],
308  grid[2][1], grid[2][2], grid[2][3],
309  grid[3][0], grid[3][1], grid[3][2],
310  grid[3][3]);
311 
312  done:
313  for (i = 0; i < 4; i++)
314  G_free(rows[i]);
315 
316  return result;
317 }
318 
319 
320 static double scancatlabel(const char *str)
321 {
322  double val;
323 
324  if (strcmp(str, "no data") != 0)
325  sscanf(str, "%lf", &val);
326  else {
327  G_warning(_("\"no data\" label found; setting to zero"));
328  val = 0.0;
329  }
330 
331  return val;
332 }
333 
void G_free(void *buf)
Free allocated memory.
Definition: gis/alloc.c:149
DCELL Rast_get_sample_cubic(int fd, const struct Cell_head *window, struct Categories *cats, double north, double east, int usedesc)
Extract a cell value from raster map (cubic interpolation).
Definition: sample.c:238
2D/3D raster map header (used also for region)
Definition: gis.h:390
double DCELL
Definition: gis.h:581
int Rast_window_rows(void)
Number of rows in active window.
Definition: raster/window.c:85
char * Rast_get_c_cat(CELL *rast, struct Categories *pcats)
Get a raster category label (CELL)
Definition: raster/cats.c:323
DCELL Rast_get_sample_nearest(int fd, const struct Cell_head *window, struct Categories *cats, double north, double east, int usedesc)
Extract a cell value from raster map (neighbor interpolation)
Definition: sample.c:91
fd
Definition: d/range.c:69
void G_fatal_error(const char *msg,...)
Print a fatal error message to stderr.
Definition: gis/error.c:160
DCELL Rast_interp_bicubic(double u, double v, DCELL c00, DCELL c01, DCELL c02, DCELL c03, DCELL c10, DCELL c11, DCELL c12, DCELL c13, DCELL c20, DCELL c21, DCELL c22, DCELL c23, DCELL c30, DCELL c31, DCELL c32, DCELL c33)
Definition: interp.c:42
DCELL Rast_interp_bilinear(double u, double v, DCELL c00, DCELL c01, DCELL c10, DCELL c11)
Definition: interp.c:26
#define INTERP_NEAREST
Definition: raster.h:20
void Rast_get_d_row(int fd, DCELL *buf, int row)
Get raster row (DCELL type)
void G_squeeze(char *line)
Remove superfluous white space.
Definition: strings.c:399
int Rast_window_cols(void)
Number of columns in active window.
void Rast_set_d_null_value(DCELL *dcellVals, int numVals)
To set a number of DCELL raster values to NULL.
Definition: null_val.c:155
DCELL * Rast_allocate_d_buf(void)
Allocates memory for a raster map of type DCELL.
Definition: alloc_cell.c:108
int INTERP_TYPE
Definition: raster.h:28
double Rast_northing_to_row(double north, const struct Cell_head *window)
Northing to row.
int CELL
Definition: gis.h:580
#define INTERP_BICUBIC
Definition: raster.h:22
#define _(str)
Definition: glocale.h:13
DCELL Rast_get_sample(int fd, const struct Cell_head *window, struct Categories *cats, double north, double east, int usedesc, INTERP_TYPE itype)
Extract a cell value from raster map.
Definition: sample.c:48
double Rast_easting_to_col(double east, const struct Cell_head *window)
Easting to column.
int
Reads the categories file for map name in mapset and stores the categories in the pcats structure...
#define INTERP_BILINEAR
Definition: raster.h:21
void G_warning(const char *msg,...)
Print a warning message to stderr.
Definition: gis/error.c:204
int Rast_is_d_null_value(const DCELL *dcellVal)
To check if a DCELL raster value is set to NULL.
Definition: null_val.c:261
DCELL Rast_get_sample_bilinear(int fd, const struct Cell_head *window, struct Categories *cats, double north, double east, int usedesc)
Extract a cell value from raster map (bilinear interpolation).
Definition: sample.c:148