GRASS GIS 7 Programmer's Manual  7.9.dev(2021)-e5379bbd7
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 void void void G_fatal_error(const char *,...) __attribute__((format(printf
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:412
#define Rast_is_d_null_value(dcellVal)
Definition: defs/raster.h:420
double DCELL
Definition: gis.h:603
void Rast_get_d_row(int, DCELL *, int)
Get raster row (DCELL type)
double Rast_northing_to_row(double, const struct Cell_head *)
Northing to row.
void G_free(void *)
Free allocated memory.
Definition: gis/alloc.c:149
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
double Rast_easting_to_col(double, const struct Cell_head *)
Easting to column.
char * Rast_get_c_cat(CELL *, struct Categories *)
Get a raster category label (CELL)
Definition: raster/cats.c:323
DCELL * Rast_allocate_d_buf(void)
Allocates memory for a raster map of type DCELL.
Definition: alloc_cell.c:108
#define INTERP_NEAREST
Definition: raster.h:20
DCELL Rast_interp_bicubic(double, double, DCELL, DCELL, DCELL, DCELL, DCELL, DCELL, DCELL, DCELL, DCELL, DCELL, DCELL, DCELL, DCELL, DCELL, DCELL, DCELL)
Definition: interp.c:42
int INTERP_TYPE
Definition: raster.h:28
int Rast_window_rows(void)
Number of rows in active window.
Definition: raster/window.c:85
int Rast_window_cols(void)
Number of columns in active window.
void G_warning(const char *,...) __attribute__((format(printf
int CELL
Definition: gis.h:602
#define INTERP_BICUBIC
Definition: raster.h:22
#define _(str)
Definition: glocale.h:10
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
#define INTERP_BILINEAR
Definition: raster.h:21
DCELL Rast_interp_bilinear(double, double, DCELL, DCELL, DCELL, DCELL)
Definition: interp.c:26
void G_squeeze(char *)
Remove superfluous white space.
Definition: strings.c:441
void Rast_set_d_null_value(DCELL *, int)
To set a number of DCELL raster values to NULL.
Definition: null_val.c:155
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