GRASS Programmer's Manual  6.5.svn(2014)-r66266
 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 }
DCELL G_get_raster_sample_nearest(int fd, const struct Cell_head *window, struct Categories *cats, double north, double east, int usedesc)
Definition: sample.c:78
DCELL G_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:24
int G_get_d_raster_row(int fd, DCELL *buf, int row)
Get raster row (DCELL type)
Definition: gis/get_row.c:965
void G_free(void *buf)
Free allocated memory.
Definition: gis/alloc.c:142
double G_easting_to_col(double east, const struct Cell_head *window)
Easting to column.
Definition: window_map.c:198
FILE * fd
Definition: g3dcolor.c:368
void G_set_d_null_value(DCELL *dcellVals, int numVals)
Definition: null_val.c:176
DCELL G_get_raster_sample_bilinear(int fd, const struct Cell_head *window, struct Categories *cats, double north, double east, int usedesc)
Definition: sample.c:122
DCELL G_interp_bilinear(double u, double v, DCELL c00, DCELL c01, DCELL c10, DCELL c11)
Definition: interp.c:9
unsigned char * grid
DCELL G_get_raster_sample_cubic(int fd, const struct Cell_head *window, struct Categories *cats, double north, double east, int usedesc)
Definition: sample.c:197
char * G_get_cat(CELL num, struct Categories *pcats)
get a category label
Definition: gis/cats.c:583
int G_is_d_null_value(const DCELL *dcellVal)
Returns 1 if dcell is NULL, 0 otherwise. This will test if the value dcell is a NaN. Same test as in G_is_f_null_value().
Definition: null_val.c:306
tuple window
Definition: tools.py:543
int
Definition: g3dcolor.c:48
int G_window_cols(void)
Number of columns in active window.
Definition: window_map.c:306
DCELL G_get_raster_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:50
char buf[GNAME_MAX+sizeof(G3D_DIRECTORY)+2]
Definition: g3drange.c:62
G_warning("category support for [%s] in mapset [%s] %s", name, mapset, type)
int G_debug(int level, const char *msg,...)
Print debugging message.
Definition: gis/debug.c:51
double G_northing_to_row(double north, const struct Cell_head *window)
Northing to row.
Definition: window_map.c:129
char * G_squeeze(char *line)
Remove superfluous white space.
Definition: squeeze.c:45
int G_fatal_error(const char *msg,...)
Print a fatal error message to stderr.
int G_window_rows(void)
Number of rows in active window.
Definition: window_map.c:273
DCELL * G_allocate_d_raster_buf(void)
Allocates memory for a raster map of type DCELL.
Definition: alloc_cell.c:126