GRASS Programmer's Manual  6.4.2(2012)
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
sample.c
Go to the documentation of this file.
1 
21 #include <string.h>
22 #include <unistd.h>
23 #include <math.h>
24 #include <grass/gis.h>
25 #include <grass/glocale.h>
26 
27 /* prototypes */
28 static double scancatlabel(const char *);
29 static void raster_row_error(const struct Cell_head *window, double north,
30  double east);
31 
32 
51  int fd,
52  const struct Cell_head *window,
53  struct Categories *cats,
54  double north, double east,
55  int usedesc, INTERP_TYPE itype)
56 {
57  double retval;
58 
59  switch (itype) {
60  case NEAREST:
61  retval = G_get_raster_sample_nearest(fd, window, cats, north, east, usedesc);
62  break;
63  case BILINEAR:
64  retval = G_get_raster_sample_bilinear(fd, window, cats, north, east, usedesc);
65  break;
66  case CUBIC:
67  retval = G_get_raster_sample_cubic(fd, window, cats, north, east, usedesc);
68  break;
69  default:
70  G_fatal_error("G_get_raster_sample: %s",
71  _("Unknown interpolation type"));
72  }
73 
74  return retval;
75 }
76 
77 
79  int fd,
80  const struct Cell_head *window,
81  struct Categories *cats,
82  double north, double east, int usedesc)
83 {
84  int row, col;
85  DCELL result;
86  DCELL *maprow = G_allocate_d_raster_buf();
87 
88  /* convert northing and easting to row and col, resp */
89  row = (int)floor(G_northing_to_row(north, window));
90  col = (int)floor(G_easting_to_col(east, window));
91 
92  if (row < 0 || row >= G_window_rows() ||
93  col < 0 || col >= G_window_cols()) {
94  G_set_d_null_value(&result, 1);
95  goto done;
96  }
97 
98  if (G_get_d_raster_row(fd, maprow, row) < 0)
99  raster_row_error(window, north, east);
100 
101  if (G_is_d_null_value(&maprow[col])) {
102  G_set_d_null_value(&result, 1);
103  goto done;
104  }
105 
106  if (usedesc) {
107  char *buf = G_get_cat(maprow[col], cats);
108 
109  G_squeeze(buf);
110  result = scancatlabel(buf);
111  }
112  else
113  result = maprow[col];
114 
115 done:
116  G_free(maprow);
117 
118  return result;
119 }
120 
121 
123  int fd,
124  const struct Cell_head *window,
125  struct Categories *cats,
126  double north, double east, int usedesc)
127 {
128  int row, col;
129  double grid[2][2];
130  DCELL *arow = G_allocate_d_raster_buf();
131  DCELL *brow = G_allocate_d_raster_buf();
132  double frow, fcol, trow, tcol;
133  DCELL result;
134 
135  frow = G_northing_to_row(north, window);
136  fcol = G_easting_to_col(east, window);
137 
138  /* convert northing and easting to row and col, resp */
139  row = (int)floor(frow - 0.5);
140  col = (int)floor(fcol - 0.5);
141 
142  trow = frow - row - 0.5;
143  tcol = fcol - col - 0.5;
144 
145  if (row < 0 || row + 1 >= G_window_rows() ||
146  col < 0 || col + 1 >= G_window_cols()) {
147  G_set_d_null_value(&result, 1);
148  goto done;
149  }
150 
151  if (G_get_d_raster_row(fd, arow, row) < 0)
152  raster_row_error(window, north, east);
153  if (G_get_d_raster_row(fd, brow, row + 1) < 0)
154  raster_row_error(window, north, east);
155 
156  if (G_is_d_null_value(&arow[col]) || G_is_d_null_value(&arow[col + 1]) ||
157  G_is_d_null_value(&brow[col]) || G_is_d_null_value(&brow[col + 1])) {
158  G_set_d_null_value(&result, 1);
159  goto done;
160  }
161 
162  /*-
163  * now were ready to do bilinear interpolation over
164  * arow[col], arow[col+1],
165  * brow[col], brow[col+1]
166  */
167 
168  if (usedesc) {
169  char *buf;
170 
171  G_squeeze(buf = G_get_cat((int)arow[col], cats));
172  grid[0][0] = scancatlabel(buf);
173  G_squeeze(buf = G_get_cat((int)arow[col + 1], cats));
174  grid[0][1] = scancatlabel(buf);
175  G_squeeze(buf = G_get_cat((int)brow[col], cats));
176  grid[1][0] = scancatlabel(buf);
177  G_squeeze(buf = G_get_cat((int)brow[col + 1], cats));
178  grid[1][1] = scancatlabel(buf);
179  }
180  else {
181  grid[0][0] = arow[col];
182  grid[0][1] = arow[col + 1];
183  grid[1][0] = brow[col];
184  grid[1][1] = brow[col + 1];
185  }
186 
187  result = G_interp_bilinear(tcol, trow,
188  grid[0][0], grid[0][1], grid[1][0], grid[1][1]);
189 
190 done:
191  G_free(arow);
192  G_free(brow);
193 
194  return result;
195 }
196 
198  int fd,
199  const struct Cell_head *window,
200  struct Categories *cats,
201  double north, double east, int usedesc)
202 {
203  int i, j, row, col;
204  double grid[4][4];
205  DCELL *rows[4];
206  double frow, fcol, trow, tcol;
207  DCELL result;
208 
209  for (i = 0; i < 4; i++)
210  rows[i] = G_allocate_d_raster_buf();
211 
212  frow = G_northing_to_row(north, window);
213  fcol = G_easting_to_col(east, window);
214 
215  /* convert northing and easting to row and col, resp */
216  row = (int)floor(frow - 1.5);
217  col = (int)floor(fcol - 1.5);
218 
219  trow = frow - row - 1.5;
220  tcol = fcol - col - 1.5;
221 
222  if (row < 0 || row + 3 >= G_window_rows() ||
223  col < 0 || col + 3 >= G_window_cols()) {
224  G_set_d_null_value(&result, 1);
225  goto done;
226  }
227 
228  for (i = 0; i < 4; i++)
229  if (G_get_d_raster_row(fd, rows[i], row + i) < 0)
230  raster_row_error(window, north, east);
231 
232  for (i = 0; i < 4; i++)
233  for (j = 0; j < 4; j++)
234  if (G_is_d_null_value(&rows[i][col + j])) {
235  G_set_d_null_value(&result, 1);
236  goto done;
237  }
238 
239  /*
240  * now were ready to do cubic interpolation over
241  * arow[col], arow[col+1], arow[col+2], arow[col+3],
242  * brow[col], brow[col+1], brow[col+2], brow[col+3],
243  * crow[col], crow[col+1], crow[col+2], crow[col+3],
244  * drow[col], drow[col+1], drow[col+2], drow[col+3],
245  */
246 
247  if (usedesc) {
248  char *buf;
249 
250  for (i = 0; i < 4; i++) {
251  for (j = 0; j < 4; j++) {
252  G_squeeze(buf = G_get_cat(rows[i][col + j], cats));
253  grid[i][j] = scancatlabel(buf);
254  }
255  }
256  }
257  else {
258  for (i = 0; i < 4; i++)
259  for (j = 0; j < 4; j++)
260  grid[i][j] = rows[i][col + j];
261  }
262 
263  result = G_interp_bicubic(tcol, trow,
264  grid[0][0], grid[0][1], grid[0][2], grid[0][3],
265  grid[1][0], grid[1][1], grid[1][2], grid[1][3],
266  grid[2][0], grid[2][1], grid[2][2], grid[2][3],
267  grid[3][0], grid[3][1], grid[3][2], grid[3][3]);
268 
269 done:
270  for (i = 0; i < 4; i++)
271  G_free(rows[i]);
272 
273  return result;
274 }
275 
276 
277 static double scancatlabel(const char *str)
278 {
279  double val;
280 
281  if (strcmp(str, "no data") != 0)
282  sscanf(str, "%lf", &val);
283  else {
284  G_warning(_("\"no data\" label found; setting to zero"));
285  val = 0.0;
286  }
287 
288  return val;
289 }
290 
291 
292 static void raster_row_error(const struct Cell_head *window, double north,
293  double east)
294 {
295  G_debug(3, "DIAG: \tRegion is: n=%g s=%g e=%g w=%g",
296  window->north, window->south, window->east, window->west);
297  G_debug(3, " \tData point is north=%g east=%g", north, east);
298 
299  G_fatal_error(_("Problem reading raster map"));
300 }