GRASS GIS 7 Programmer's Manual  7.5.svn(2017)-r71933
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
raster/get_row.c
Go to the documentation of this file.
1 /*!
2  \file lib/raster/get_row.c
3 
4  \brief Raster library - Get raster row
5 
6  (C) 2003-2009 by the GRASS Development Team
7 
8  This program is free software under the GNU General Public License
9  (>=v2). Read the file COPYING that comes with GRASS for details.
10 
11  \author Original author CERL
12  */
13 
14 #include <string.h>
15 #include <unistd.h>
16 #include <sys/types.h>
17 
18 #include <grass/config.h>
19 #include <grass/raster.h>
20 #include <grass/glocale.h>
21 
22 #include "R.h"
23 
24 static void embed_nulls(int, void *, int, RASTER_MAP_TYPE, int, int);
25 
26 static int compute_window_row(int fd, int row, int *cellRow)
27 {
28  struct fileinfo *fcb = &R__.fileinfo[fd];
29  double f;
30  int r;
31 
32  /* check for row in window */
33  if (row < 0 || row >= R__.rd_window.rows) {
34  G_fatal_error(_("Reading raster map <%s@%s> request for row %d is outside region"),
35  fcb->name, fcb->mapset, row);
36  }
37 
38  /* convert window row to cell file row */
39  f = row * fcb->C1 + fcb->C2;
40  r = (int)f;
41  if (f < r) /* adjust for rounding up of negatives */
42  r--;
43 
44  if (r < 0 || r >= fcb->cellhd.rows)
45  return 0;
46 
47  *cellRow = r;
48 
49  return 1;
50 }
51 
52 static void do_reclass_int(int fd, void *cell, int null_is_zero)
53 {
54  struct fileinfo *fcb = &R__.fileinfo[fd];
55  CELL *c = cell;
56  CELL *reclass_table = fcb->reclass.table;
57  CELL min = fcb->reclass.min;
58  CELL max = fcb->reclass.max;
59  int i;
60 
61  for (i = 0; i < R__.rd_window.cols; i++) {
62  if (Rast_is_c_null_value(&c[i])) {
63  if (null_is_zero)
64  c[i] = 0;
65  continue;
66  }
67 
68  if (c[i] < min || c[i] > max) {
69  if (null_is_zero)
70  c[i] = 0;
71  else
72  Rast_set_c_null_value(&c[i], 1);
73  continue;
74  }
75 
76  c[i] = reclass_table[c[i] - min];
77 
78  if (null_is_zero && Rast_is_c_null_value(&c[i]))
79  c[i] = 0;
80  }
81 }
82 
83 static void read_data_fp_compressed(int fd, int row, unsigned char *data_buf,
84  int *nbytes)
85 {
86  struct fileinfo *fcb = &R__.fileinfo[fd];
87  off_t t1 = fcb->row_ptr[row];
88  off_t t2 = fcb->row_ptr[row + 1];
89  size_t readamount = t2 - t1;
90  size_t bufsize = fcb->cellhd.cols * fcb->nbytes;
91 
92  if (lseek(fcb->data_fd, t1, SEEK_SET) < 0)
93  G_fatal_error(_("Error reading raster data for row %d of <%s>"),
94  row, fcb->name);
95 
96  *nbytes = fcb->nbytes;
97 
98  if ((size_t) G_read_compressed(fcb->data_fd, readamount, data_buf,
99  bufsize, fcb->cellhd.compressed) != bufsize)
100  G_fatal_error(_("Error uncompressing raster data for row %d of <%s>"),
101  row, fcb->name);
102 }
103 
104 static void rle_decompress(unsigned char *dst, const unsigned char *src,
105  int nbytes, int size)
106 {
107  int pairs = size / (nbytes + 1);
108  int i;
109 
110  for (i = 0; i < pairs; i++) {
111  int repeat = *src++;
112  int j;
113 
114  for (j = 0; j < repeat; j++) {
115  memcpy(dst, src, nbytes);
116  dst += nbytes;
117  }
118 
119  src += nbytes;
120  }
121 }
122 
123 static void read_data_compressed(int fd, int row, unsigned char *data_buf,
124  int *nbytes)
125 {
126  struct fileinfo *fcb = &R__.fileinfo[fd];
127  off_t t1 = fcb->row_ptr[row];
128  off_t t2 = fcb->row_ptr[row + 1];
129  ssize_t readamount = t2 - t1;
130  size_t bufsize;
131  unsigned char *cmp, *cmp2;
132  int n;
133 
134  if (lseek(fcb->data_fd, t1, SEEK_SET) < 0)
135  G_fatal_error(_("Error reading raster data for row %d of <%s>"),
136  row, fcb->name);
137 
138  cmp = G_malloc(readamount);
139 
140  if (read(fcb->data_fd, cmp, readamount) != readamount) {
141  G_free(cmp);
142  G_fatal_error(_("Error reading raster data for row %d of <%s>"),
143  row, fcb->name);
144  }
145 
146  /* save cmp for free below */
147  cmp2 = cmp;
148 
149  /* Now decompress the row */
150  if (fcb->cellhd.compressed > 0) {
151  /* one byte is nbyte count */
152  n = *nbytes = *cmp++;
153  readamount--;
154  }
155  else
156  /* pre 3.0 compression */
157  n = *nbytes = fcb->nbytes;
158 
159  bufsize = n * fcb->cellhd.cols;
160  if (fcb->cellhd.compressed < 0 || readamount < bufsize) {
161  if (fcb->cellhd.compressed == 1)
162  rle_decompress(data_buf, cmp, n, readamount);
163  else {
164  if (G_expand(cmp, readamount, data_buf, bufsize,
165  fcb->cellhd.compressed) != bufsize)
166  G_fatal_error(_("Error uncompressing raster data for row %d of <%s>"),
167  row, fcb->name);
168  }
169  }
170  else
171  memcpy(data_buf, cmp, readamount);
172 
173  G_free(cmp2);
174 }
175 
176 static void read_data_uncompressed(int fd, int row, unsigned char *data_buf,
177  int *nbytes)
178 {
179  struct fileinfo *fcb = &R__.fileinfo[fd];
180  ssize_t bufsize = fcb->cellhd.cols * fcb->nbytes;
181 
182  *nbytes = fcb->nbytes;
183 
184  if (lseek(fcb->data_fd, (off_t) row * bufsize, SEEK_SET) == -1)
185  G_fatal_error(_("Error reading raster data for row %d of <%s>"),
186  row, fcb->name);
187 
188  if (read(fcb->data_fd, data_buf, bufsize) != bufsize)
189  G_fatal_error(_("Error reading raster data for row %d of <%s>"),
190  row, fcb->name);
191 }
192 
193 #ifdef HAVE_GDAL
194 static void read_data_gdal(int fd, int row, unsigned char *data_buf,
195  int *nbytes)
196 {
197  struct fileinfo *fcb = &R__.fileinfo[fd];
198  unsigned char *buf;
199  CPLErr err;
200 
201  *nbytes = fcb->nbytes;
202 
203  if (fcb->gdal->vflip)
204  row = fcb->cellhd.rows - 1 - row;
205 
206  buf = fcb->gdal->hflip ? G_malloc(fcb->cellhd.cols * fcb->cur_nbytes)
207  : data_buf;
208 
209  err =
210  Rast_gdal_raster_IO(fcb->gdal->band, GF_Read, 0, row,
211  fcb->cellhd.cols, 1, buf, fcb->cellhd.cols, 1,
212  fcb->gdal->type, 0, 0);
213 
214  if (fcb->gdal->hflip) {
215  int i;
216 
217  for (i = 0; i < fcb->cellhd.cols; i++)
218  memcpy(data_buf + i * fcb->cur_nbytes,
219  buf + (fcb->cellhd.cols - 1 - i) * fcb->cur_nbytes,
220  fcb->cur_nbytes);
221  G_free(buf);
222  }
223 
224  if (err != CE_None)
225  G_fatal_error(_("Error reading raster data via GDAL for row %d of <%s>"),
226  row, fcb->name);
227 }
228 #endif
229 
230 static void read_data(int fd, int row, unsigned char *data_buf, int *nbytes)
231 {
232  struct fileinfo *fcb = &R__.fileinfo[fd];
233 
234 #ifdef HAVE_GDAL
235  if (fcb->gdal) {
236  read_data_gdal(fd, row, data_buf, nbytes);
237  return;
238  }
239 #endif
240 
241  if (!fcb->cellhd.compressed)
242  read_data_uncompressed(fd, row, data_buf, nbytes);
243  else if (fcb->map_type == CELL_TYPE)
244  read_data_compressed(fd, row, data_buf, nbytes);
245  else
246  read_data_fp_compressed(fd, row, data_buf, nbytes);
247 }
248 
249 /* copy cell file data to user buffer translated by window column mapping */
250 static void cell_values_int(int fd, const unsigned char *data,
251  const COLUMN_MAPPING * cmap, int nbytes,
252  void *cell, int n)
253 {
254  CELL *c = cell;
255  COLUMN_MAPPING cmapold = 0;
256  int big = (size_t) nbytes >= sizeof(CELL);
257  int i;
258 
259  for (i = 0; i < n; i++) {
260  const unsigned char *d;
261  int neg;
262  CELL v;
263  int j;
264 
265  if (!cmap[i]) {
266  c[i] = 0;
267  continue;
268  }
269 
270  if (cmap[i] == cmapold) {
271  c[i] = c[i - 1];
272  continue;
273  }
274 
275  d = data + (cmap[i] - 1) * nbytes;
276 
277  if (big && (*d & 0x80)) {
278  neg = 1;
279  v = *d++ & 0x7f;
280  }
281  else {
282  neg = 0;
283  v = *d++;
284  }
285 
286  for (j = 1; j < nbytes; j++)
287  v = (v << 8) + *d++;
288 
289  c[i] = neg ? -v : v;
290 
291  cmapold = cmap[i];
292  }
293 }
294 
295 static void cell_values_float(int fd, const unsigned char *data,
296  const COLUMN_MAPPING * cmap, int nbytes,
297  void *cell, int n)
298 {
299  struct fileinfo *fcb = &R__.fileinfo[fd];
300  const float *work_buf = (const float *) fcb->data;
301  FCELL *c = cell;
302  int i;
303 
304  for (i = 0; i < n; i++) {
305  if (!cmap[i]) {
306  c[i] = 0;
307  continue;
308  }
309 
310  G_xdr_get_float(&c[i], &work_buf[cmap[i]-1]);
311  }
312 }
313 
314 static void cell_values_double(int fd, const unsigned char *data,
315  const COLUMN_MAPPING * cmap, int nbytes,
316  void *cell, int n)
317 {
318  struct fileinfo *fcb = &R__.fileinfo[fd];
319  const double *work_buf = (const double *) fcb->data;
320  DCELL *c = cell;
321  int i;
322 
323  for (i = 0; i < n; i++) {
324  if (!cmap[i]) {
325  c[i] = 0;
326  continue;
327  }
328 
329  G_xdr_get_double(&c[i], &work_buf[cmap[i]-1]);
330  }
331 }
332 
333 #ifdef HAVE_GDAL
334 static void gdal_values_int(int fd, const unsigned char *data,
335  const COLUMN_MAPPING * cmap, int nbytes,
336  CELL * cell, int n)
337 {
338  struct fileinfo *fcb = &R__.fileinfo[fd];
339  const unsigned char *d;
340  COLUMN_MAPPING cmapold = 0;
341  int i;
342 
343  for (i = 0; i < n; i++) {
344  if (!cmap[i]) {
345  cell[i] = 0;
346  continue;
347  }
348 
349  if (cmap[i] == cmapold) {
350  cell[i] = cell[i - 1];
351  continue;
352  }
353 
354  d = data + (cmap[i] - 1) * nbytes;
355 
356  switch (fcb->gdal->type) {
357  case GDT_Byte:
358  cell[i] = *(GByte *) d;
359  break;
360  case GDT_Int16:
361  cell[i] = *(GInt16 *) d;
362  break;
363  case GDT_UInt16:
364  cell[i] = *(GUInt16 *) d;
365  break;
366  case GDT_Int32:
367  cell[i] = *(GInt32 *) d;
368  break;
369  case GDT_UInt32:
370  cell[i] = *(GUInt32 *) d;
371  break;
372  default:
373  /* shouldn't happen */
374  Rast_set_c_null_value(&cell[i], 1);
375  break;
376  }
377 
378  cmapold = cmap[i];
379  }
380 }
381 
382 static void gdal_values_float(int fd, const float *data,
383  const COLUMN_MAPPING * cmap, int nbytes,
384  FCELL * cell, int n)
385 {
386  COLUMN_MAPPING cmapold = 0;
387  int i;
388 
389  for (i = 0; i < n; i++) {
390  if (!cmap[i]) {
391  cell[i] = 0;
392  continue;
393  }
394 
395  if (cmap[i] == cmapold) {
396  cell[i] = cell[i - 1];
397  continue;
398  }
399 
400  cell[i] = data[cmap[i] - 1];
401 
402  cmapold = cmap[i];
403  }
404 }
405 
406 static void gdal_values_double(int fd, const double *data,
407  const COLUMN_MAPPING * cmap, int nbytes,
408  DCELL * cell, int n)
409 {
410  COLUMN_MAPPING cmapold = 0;
411  int i;
412 
413  for (i = 0; i < n; i++) {
414  if (!cmap[i]) {
415  cell[i] = 0;
416  continue;
417  }
418 
419  if (cmap[i] == cmapold) {
420  cell[i] = cell[i - 1];
421  continue;
422  }
423 
424  cell[i] = data[cmap[i] - 1];
425 
426  cmapold = cmap[i];
427  }
428 }
429 #endif
430 
431 /* transfer_to_cell_XY takes bytes from fcb->data, converts these bytes with
432  the appropriate procedure (e.g. XDR or byte reordering) into type X
433  values which are put into array work_buf.
434  finally the values in work_buf are converted into
435  type Y and put into 'cell'.
436  if type X == type Y the intermediate step of storing the values in
437  work_buf might be omitted. check the appropriate function for XY to
438  determine the procedure of conversion.
439  */
440 static void transfer_to_cell_XX(int fd, void *cell)
441 {
442  static void (*cell_values_type[3]) () = {
443  cell_values_int, cell_values_float, cell_values_double};
444 #ifdef HAVE_GDAL
445  static void (*gdal_values_type[3]) () = {
446  gdal_values_int, gdal_values_float, gdal_values_double};
447 #endif
448  struct fileinfo *fcb = &R__.fileinfo[fd];
449 
450 #ifdef HAVE_GDAL
451  if (fcb->gdal)
452  (gdal_values_type[fcb->map_type]) (fd, fcb->data, fcb->col_map,
453  fcb->cur_nbytes, cell,
454  R__.rd_window.cols);
455  else
456 #endif
457  (cell_values_type[fcb->map_type]) (fd, fcb->data, fcb->col_map,
458  fcb->cur_nbytes, cell,
459  R__.rd_window.cols);
460 }
461 
462 static void transfer_to_cell_fi(int fd, void *cell)
463 {
464  struct fileinfo *fcb = &R__.fileinfo[fd];
465  FCELL *work_buf = G_malloc(R__.rd_window.cols * sizeof(FCELL));
466  int i;
467 
468  transfer_to_cell_XX(fd, work_buf);
469 
470  for (i = 0; i < R__.rd_window.cols; i++)
471  ((CELL *) cell)[i] = (fcb->col_map[i] == 0)
472  ? 0 : Rast_quant_get_cell_value(&fcb->quant, work_buf[i]);
473 
474  G_free(work_buf);
475 }
476 
477 static void transfer_to_cell_di(int fd, void *cell)
478 {
479  struct fileinfo *fcb = &R__.fileinfo[fd];
480  DCELL *work_buf = G_malloc(R__.rd_window.cols * sizeof(DCELL));
481  int i;
482 
483  transfer_to_cell_XX(fd, work_buf);
484 
485  for (i = 0; i < R__.rd_window.cols; i++)
486  ((CELL *) cell)[i] = (fcb->col_map[i] == 0)
487  ? 0 : Rast_quant_get_cell_value(&fcb->quant, work_buf[i]);
488 
489  G_free(work_buf);
490 }
491 
492 static void transfer_to_cell_if(int fd, void *cell)
493 {
494  CELL *work_buf = G_malloc(R__.rd_window.cols * sizeof(CELL));
495  int i;
496 
497  transfer_to_cell_XX(fd, work_buf);
498 
499  for (i = 0; i < R__.rd_window.cols; i++)
500  ((FCELL *) cell)[i] = work_buf[i];
501 
502  G_free(work_buf);
503 }
504 
505 static void transfer_to_cell_df(int fd, void *cell)
506 {
507  DCELL *work_buf = G_malloc(R__.rd_window.cols * sizeof(DCELL));
508  int i;
509 
510  transfer_to_cell_XX(fd, work_buf);
511 
512  for (i = 0; i < R__.rd_window.cols; i++)
513  ((FCELL *) cell)[i] = work_buf[i];
514 
515  G_free(work_buf);
516 }
517 
518 static void transfer_to_cell_id(int fd, void *cell)
519 {
520  CELL *work_buf = G_malloc(R__.rd_window.cols * sizeof(CELL));
521  int i;
522 
523  transfer_to_cell_XX(fd, work_buf);
524 
525  for (i = 0; i < R__.rd_window.cols; i++)
526  ((DCELL *) cell)[i] = work_buf[i];
527 
528  G_free(work_buf);
529 }
530 
531 static void transfer_to_cell_fd(int fd, void *cell)
532 {
533  FCELL *work_buf = G_malloc(R__.rd_window.cols * sizeof(FCELL));
534  int i;
535 
536  transfer_to_cell_XX(fd, work_buf);
537 
538  for (i = 0; i < R__.rd_window.cols; i++)
539  ((DCELL *) cell)[i] = work_buf[i];
540 
541  G_free(work_buf);
542 }
543 
544 /*
545  * works for all map types and doesn't consider
546  * null row corresponding to the requested row
547  */
548 static int get_map_row_nomask(int fd, void *rast, int row,
549  RASTER_MAP_TYPE data_type)
550 {
551  static void (*transfer_to_cell_FtypeOtype[3][3]) () = {
552  {
553  transfer_to_cell_XX, transfer_to_cell_if, transfer_to_cell_id}, {
554  transfer_to_cell_fi, transfer_to_cell_XX, transfer_to_cell_fd}, {
555  transfer_to_cell_di, transfer_to_cell_df, transfer_to_cell_XX}
556  };
557  struct fileinfo *fcb = &R__.fileinfo[fd];
558  int r;
559  int row_status;
560 
561  row_status = compute_window_row(fd, row, &r);
562 
563  if (!row_status) {
564  fcb->cur_row = -1;
565  Rast_zero_input_buf(rast, data_type);
566  return 0;
567  }
568 
569  /* read cell file row if not in memory */
570  if (r != fcb->cur_row) {
571  fcb->cur_row = r;
572  read_data(fd, fcb->cur_row, fcb->data, &fcb->cur_nbytes);
573  }
574 
575  (transfer_to_cell_FtypeOtype[fcb->map_type][data_type]) (fd, rast);
576 
577  return 1;
578 }
579 
580 static void get_map_row_no_reclass(int fd, void *rast, int row,
581  RASTER_MAP_TYPE data_type, int null_is_zero,
582  int with_mask)
583 {
584  get_map_row_nomask(fd, rast, row, data_type);
585  embed_nulls(fd, rast, row, data_type, null_is_zero, with_mask);
586 }
587 
588 static void get_map_row(int fd, void *rast, int row, RASTER_MAP_TYPE data_type,
589  int null_is_zero, int with_mask)
590 {
591  struct fileinfo *fcb = &R__.fileinfo[fd];
592  int size = Rast_cell_size(data_type);
593  CELL *temp_buf = NULL;
594  void *buf;
595  int type;
596  int i;
597 
598  if (fcb->reclass_flag && data_type != CELL_TYPE) {
599  temp_buf = G_malloc(R__.rd_window.cols * sizeof(CELL));
600  buf = temp_buf;
601  type = CELL_TYPE;
602  }
603  else {
604  buf = rast;
605  type = data_type;
606  }
607 
608  get_map_row_no_reclass(fd, buf, row, type, null_is_zero, with_mask);
609 
610  if (!fcb->reclass_flag)
611  return;
612 
613  /* if the map is reclass table, get and
614  reclass CELL row and copy results to needed type */
615 
616  do_reclass_int(fd, buf, null_is_zero);
617 
618  if (data_type == CELL_TYPE)
619  return;
620 
621  for (i = 0; i < R__.rd_window.cols; i++) {
622  Rast_set_c_value(rast, temp_buf[i], data_type);
623  rast = G_incr_void_ptr(rast, size);
624  }
625 
626  if (fcb->reclass_flag && data_type != CELL_TYPE) {
627  G_free(temp_buf);
628  }
629 }
630 
631 /*!
632  * \brief Read raster row without masking
633  *
634  * This routine reads the specified <em>row</em> from the raster map
635  * open on file descriptor <em>fd</em> into the <em>buf</em> buffer
636  * like Rast_get_c_row() does. The difference is that masking is
637  * suppressed. If the user has a mask set, Rast_get_c_row() will apply
638  * the mask but Rast_get_c_row_nomask() will ignore it. This routine
639  * prints a diagnostic message and returns -1 if there is an error
640  * reading the raster map. Otherwise a nonnegative value is returned.
641  *
642  * <b>Note.</b> Ignoring the mask is not generally acceptable. Users
643  * expect the mask to be applied. However, in some cases ignoring the
644  * mask is justified. For example, the GRASS modules
645  * <i>r.describe</i>, which reads the raster map directly to report
646  * all data values in a raster map, and <i>r.slope.aspect</i>, which
647  * produces slope and aspect from elevation, ignore both the mask and
648  * the region. However, the number of GRASS modules which do this
649  * should be minimal. See Mask for more information about the mask.
650  *
651  * \param fd file descriptor for the opened raster map
652  * \param buf buffer for the row to be placed into
653  * \param row data row desired
654  * \param data_type data type
655  *
656  * \return void
657  */
658 void Rast_get_row_nomask(int fd, void *buf, int row, RASTER_MAP_TYPE data_type)
659 {
660  get_map_row(fd, buf, row, data_type, 0, 0);
661 }
662 
663 /*!
664  * \brief Read raster row without masking (CELL type)
665  *
666  * Same as Rast_get_c_row() except no masking occurs.
667  *
668  * \param fd file descriptor for the opened raster map
669  * \param buf buffer for the row to be placed into
670  * \param row data row desired
671  *
672  * \return void
673  */
674 void Rast_get_c_row_nomask(int fd, CELL * buf, int row)
675 {
676  Rast_get_row_nomask(fd, buf, row, CELL_TYPE);
677 }
678 
679 /*!
680  * \brief Read raster row without masking (FCELL type)
681  *
682  * Same as Rast_get_f_row() except no masking occurs.
683  *
684  * \param fd file descriptor for the opened raster map
685  * \param buf buffer for the row to be placed into
686  * \param row data row desired
687  *
688  * \return void
689  */
690 void Rast_get_f_row_nomask(int fd, FCELL * buf, int row)
691 {
692  Rast_get_row_nomask(fd, buf, row, FCELL_TYPE);
693 }
694 
695 /*!
696  * \brief Read raster row without masking (DCELL type)
697  *
698  * Same as Rast_get_d_row() except no masking occurs.
699  *
700  * \param fd file descriptor for the opened raster map
701  * \param buf buffer for the row to be placed into
702  * \param row data row desired
703  *
704  * \return void
705  */
706 void Rast_get_d_row_nomask(int fd, DCELL * buf, int row)
707 {
708  Rast_get_row_nomask(fd, buf, row, DCELL_TYPE);
709 }
710 
711 /*!
712  * \brief Get raster row
713  *
714  * If <em>data_type</em> is
715  * - CELL_TYPE, calls Rast_get_c_row()
716  * - FCELL_TYPE, calls Rast_get_f_row()
717  * - DCELL_TYPE, calls Rast_get_d_row()
718  *
719  * Reads appropriate information into the buffer <em>buf</em> associated
720  * with the requested row <em>row</em>. <em>buf</em> is associated with the
721  * current window.
722  *
723  * Note, that the type of the data in <em>buf</em> (say X) is independent of
724  * the type of the data in the file described by <em>fd</em> (say Y).
725  *
726  * - Step 1: Read appropriate raw map data into a intermediate buffer.
727  * - Step 2: Convert the data into a CPU readable format, and subsequently
728  * resample the data. the data is stored in a second intermediate
729  * buffer (the type of the data in this buffer is Y).
730  * - Step 3: Convert this type Y data into type X data and store it in
731  * buffer "buf". Conversion is performed in functions
732  * "transfer_to_cell_XY". (For details of the conversion between
733  * two particular types check the functions).
734  * - Step 4: read or simmulate null value row and zero out cells corresponding
735  * to null value cells. The masked out cells are set to null when the
736  * mask exists. (the MASK is taken care of by null values
737  * (if the null file doesn't exist for this map, then the null row
738  * is simulated by assuming that all zero are nulls *** in case
739  * of Rast_get_row() and assuming that all data is valid
740  * in case of G_get_f/d_raster_row(). In case of deprecated function
741  * Rast_get_c_row() all nulls are converted to zeros (so there are
742  * no embedded nulls at all). Also all masked out cells become zeros.
743  *
744  * \param fd file descriptor for the opened raster map
745  * \param buf buffer for the row to be placed into
746  * \param row data row desired
747  * \param data_type data type
748  *
749  * \return void
750  */
751 void Rast_get_row(int fd, void *buf, int row, RASTER_MAP_TYPE data_type)
752 {
753  get_map_row(fd, buf, row, data_type, 0, 1);
754 }
755 
756 /*!
757  * \brief Get raster row (CELL type)
758  *
759  * Reads a row of raster data and leaves the NULL values intact. (As
760  * opposed to the deprecated function Rast_get_c_row() which
761  * converts NULL values to zero.)
762  *
763  * <b>NOTE.</b> When the raster map is old and null file doesn't
764  * exist, it is assumed that all 0-cells are no-data. When map is
765  * floating point, uses quant rules set explicitly by
766  * Rast_set_quant_rules() or stored in map's quant file to convert floats
767  * to integers.
768  *
769  * \param fd file descriptor for the opened raster map
770  * \param buf buffer for the row to be placed into
771  * \param row data row desired
772  *
773  * \return void
774  */
775 void Rast_get_c_row(int fd, CELL * buf, int row)
776 {
777  Rast_get_row(fd, buf, row, CELL_TYPE);
778 }
779 
780 /*!
781  * \brief Get raster row (FCELL type)
782  *
783  * Read a row from the raster map open on <em>fd</em> into the
784  * <tt>float</tt> array <em>fcell</em> performing type conversions as
785  * necessary based on the actual storage type of the map. Masking,
786  * resampling into the current region. NULL-values are always
787  * embedded in <tt>fcell</tt> (<em>never converted to a value</em>).
788  *
789  * \param fd file descriptor for the opened raster map
790  * \param buf buffer for the row to be placed into
791  * \param row data row desired
792  *
793  * \return void
794  */
795 void Rast_get_f_row(int fd, FCELL * buf, int row)
796 {
797  Rast_get_row(fd, buf, row, FCELL_TYPE);
798 }
799 
800 /*!
801  * \brief Get raster row (DCELL type)
802  *
803  * Same as Rast_get_f_row() except that the array <em>dcell</em>
804  * is <tt>double</tt>.
805  *
806  * \param fd file descriptor for the opened raster map
807  * \param buf buffer for the row to be placed into
808  * \param row data row desired
809  *
810  * \return void
811  */
812 void Rast_get_d_row(int fd, DCELL * buf, int row)
813 {
814  Rast_get_row(fd, buf, row, DCELL_TYPE);
815 }
816 
817 static int read_null_bits_compressed(int null_fd, unsigned char *flags,
818  int row, size_t size, int fd)
819 {
820  struct fileinfo *fcb = &R__.fileinfo[fd];
821  off_t t1 = fcb->null_row_ptr[row];
822  off_t t2 = fcb->null_row_ptr[row + 1];
823  size_t readamount = t2 - t1;
824  unsigned char *compressed_buf;
825 
826  if (lseek(null_fd, t1, SEEK_SET) < 0)
827  G_fatal_error(_("Error reading null data for row %d of <%s>"),
828  row, fcb->name);
829 
830  if (readamount == size) {
831  if (read(null_fd, flags, size) != size) {
832  G_fatal_error(_("Error reading null data for row %d of <%s>"),
833  row, fcb->name);
834  }
835  return 1;
836  }
837 
838  compressed_buf = G_malloc(readamount);
839 
840  if (read(null_fd, compressed_buf, readamount) != readamount) {
841  G_free(compressed_buf);
842  G_fatal_error(_("Error reading null data for row %d of <%s>"),
843  row, fcb->name);
844  }
845 
846  /* null bits file compressed with LZ4, see lib/gis/compress.h */
847  if (G_lz4_expand(compressed_buf, readamount, flags, size) < 1) {
848  G_fatal_error(_("Error uncompressing null data for row %d of <%s>"),
849  row, fcb->name);
850  }
851 
852  G_free(compressed_buf);
853 
854  return 1;
855 }
856 
857 int Rast__read_null_bits(int fd, int row, unsigned char *flags)
858 {
859  struct fileinfo *fcb = &R__.fileinfo[fd];
860  int null_fd = fcb->null_fd;
861  int cols = fcb->cellhd.cols;
862  off_t offset;
863  ssize_t size;
864  int R;
865 
866  if (compute_window_row(fd, row, &R) <= 0) {
867  Rast__init_null_bits(flags, cols);
868  return 1;
869  }
870 
871  if (null_fd < 0)
872  return 0;
873 
874  size = Rast__null_bitstream_size(cols);
875 
876  if (fcb->null_row_ptr)
877  return read_null_bits_compressed(null_fd, flags, R, size, fd);
878 
879  offset = (off_t) size * R;
880 
881  if (lseek(null_fd, offset, SEEK_SET) < 0)
882  G_fatal_error(_("Error reading null row %d for <%s>"), R, fcb->name);
883 
884  if (read(null_fd, flags, size) != size)
885  G_fatal_error(_("Error reading null row %d for <%s>"), R, fcb->name);
886 
887  return 1;
888 }
889 
890 #define check_null_bit(flags, bit_num) ((flags)[(bit_num)>>3] & ((unsigned char)0x80>>((bit_num)&7)) ? 1 : 0)
891 
892 static void get_null_value_row_nomask(int fd, char *flags, int row)
893 {
894  struct fileinfo *fcb = &R__.fileinfo[fd];
895  int j;
896 
897  if (row > R__.rd_window.rows || row < 0) {
898  G_warning(_("Reading raster map <%s@%s> request for row %d is outside region"),
899  fcb->name, fcb->mapset, row);
900  for (j = 0; j < R__.rd_window.cols; j++)
901  flags[j] = 1;
902  return;
903  }
904 
905  if (row != fcb->null_cur_row) {
906  if (!Rast__read_null_bits(fd, row, fcb->null_bits)) {
907  fcb->null_cur_row = -1;
908  if (fcb->map_type == CELL_TYPE) {
909  /* If can't read null row, assume that all map 0's are nulls */
910  CELL *mask_buf = G_malloc(R__.rd_window.cols * sizeof(CELL));
911 
912  get_map_row_nomask(fd, mask_buf, row, CELL_TYPE);
913  for (j = 0; j < R__.rd_window.cols; j++)
914  flags[j] = (mask_buf[j] == 0);
915 
916  G_free(mask_buf);
917  }
918  else { /* fp map */
919  /* if can't read null row, assume that all data is valid */
920  G_zero(flags, sizeof(char) * R__.rd_window.cols);
921  /* the flags row is ready now */
922  }
923 
924  return;
925  } /*if no null file */
926  else
927  fcb->null_cur_row = row;
928  }
929 
930  /* copy null row to flags row translated by window column mapping */
931  for (j = 0; j < R__.rd_window.cols; j++) {
932  if (!fcb->col_map[j])
933  flags[j] = 1;
934  else
935  flags[j] = check_null_bit(fcb->null_bits, fcb->col_map[j] - 1);
936  }
937 }
938 
939 /*--------------------------------------------------------------------------*/
940 
941 #ifdef HAVE_GDAL
942 
943 static void get_null_value_row_gdal(int fd, char *flags, int row)
944 {
945  struct fileinfo *fcb = &R__.fileinfo[fd];
946  DCELL *tmp_buf = Rast_allocate_d_input_buf();
947  int i;
948 
949  if (get_map_row_nomask(fd, tmp_buf, row, DCELL_TYPE) <= 0) {
950  memset(flags, 1, R__.rd_window.cols);
951  G_free(tmp_buf);
952  return;
953  }
954 
955  for (i = 0; i < R__.rd_window.cols; i++)
956  /* note: using == won't work if the null value is NaN */
957  flags[i] = !fcb->col_map[i] ||
958  tmp_buf[i] == fcb->gdal->null_val ||
959  tmp_buf[i] != tmp_buf[i];
960 
961  G_free(tmp_buf);
962 }
963 
964 #endif
965 
966 /*--------------------------------------------------------------------------*/
967 
968 /*--------------------------------------------------------------------------*/
969 
970 static void embed_mask(char *flags, int row)
971 {
972  CELL *mask_buf = G_malloc(R__.rd_window.cols * sizeof(CELL));
973  int i;
974 
975  if (R__.auto_mask <= 0) {
976  G_free(mask_buf);
977  return;
978  }
979 
980  if (get_map_row_nomask(R__.mask_fd, mask_buf, row, CELL_TYPE) < 0) {
981  G_free(mask_buf);
982  return;
983  }
984 
986  embed_nulls(R__.mask_fd, mask_buf, row, CELL_TYPE, 0, 0);
987  do_reclass_int(R__.mask_fd, mask_buf, 1);
988  }
989 
990  for (i = 0; i < R__.rd_window.cols; i++)
991  if (mask_buf[i] == 0 || Rast_is_c_null_value(&mask_buf[i]))
992  flags[i] = 1;
993 
994  G_free(mask_buf);
995 }
996 
997 static void get_null_value_row(int fd, char *flags, int row, int with_mask)
998 {
999 #ifdef HAVE_GDAL
1000  struct fileinfo *fcb = &R__.fileinfo[fd];
1001 
1002  if (fcb->gdal)
1003  get_null_value_row_gdal(fd, flags, row);
1004  else
1005 #endif
1006  get_null_value_row_nomask(fd, flags, row);
1007 
1008  if (with_mask)
1009  embed_mask(flags, row);
1010 }
1011 
1012 static void embed_nulls(int fd, void *buf, int row, RASTER_MAP_TYPE map_type,
1013  int null_is_zero, int with_mask)
1014 {
1015  struct fileinfo *fcb = &R__.fileinfo[fd];
1016  size_t size = Rast_cell_size(map_type);
1017  char *null_buf;
1018  int i;
1019 
1020  /* this is because without null file the nulls can be only due to 0's
1021  in data row or mask */
1022  if (null_is_zero && !fcb->null_file_exists
1023  && (R__.auto_mask <= 0 || !with_mask))
1024  return;
1025 
1026  null_buf = G_malloc(R__.rd_window.cols);
1027 
1028  get_null_value_row(fd, null_buf, row, with_mask);
1029 
1030  for (i = 0; i < R__.rd_window.cols; i++) {
1031  /* also check for nulls which might be already embedded by quant
1032  rules in case of fp map. */
1033  if (null_buf[i] || Rast_is_null_value(buf, map_type)) {
1034  /* G__set_[f/d]_null_value() sets it to 0 is the embedded mode
1035  is not set and calls G_set_[f/d]_null_value() otherwise */
1036  Rast__set_null_value(buf, 1, null_is_zero, map_type);
1037  }
1038  buf = G_incr_void_ptr(buf, size);
1039  }
1040 
1041  G_free(null_buf);
1042 }
1043 
1044 /*!
1045  \brief Read or simulate null value row
1046 
1047  Read or simulate null value row and set the cells corresponding
1048  to null value to 1. The masked out cells are set to null when the
1049  mask exists. (the MASK is taken care of by null values
1050  (if the null file doesn't exist for this map, then the null row
1051  is simulated by assuming that all zeros in raster map are nulls.
1052  Also all masked out cells become nulls.
1053 
1054  \param fd file descriptor for the opened map
1055  \param buf buffer for the row to be placed into
1056  \param flags
1057  \param row data row desired
1058 
1059  \return void
1060  */
1061 void Rast_get_null_value_row(int fd, char *flags, int row)
1062 {
1063  struct fileinfo *fcb = &R__.fileinfo[fd];
1064 
1065  if (!fcb->reclass_flag)
1066  get_null_value_row(fd, flags, row, 1);
1067  else {
1068  CELL *buf = G_malloc(R__.rd_window.cols * sizeof(CELL));
1069  int i;
1070 
1071  Rast_get_c_row(fd, buf, row);
1072  for (i = 0; i < R__.rd_window.cols; i++)
1073  flags[i] = Rast_is_c_null_value(&buf[i]) ? 1 : 0;
1074 
1075  G_free(buf);
1076  }
1077 }
#define CELL_TYPE
Definition: raster.h:11
struct Cell_head rd_window
Definition: R.h:83
void Rast_get_row_nomask(int fd, void *buf, int row, RASTER_MAP_TYPE data_type)
Read raster row without masking.
int G_read_compressed(int fd, int rbytes, unsigned char *dst, int nbytes, int number)
Definition: compress.c:248
int nbytes
Definition: R.h:58
int G_lz4_expand(unsigned char *src, int src_sz, unsigned char *dst, int dst_sz)
Definition: cmprlz4.c:136
void G_free(void *buf)
Free allocated memory.
Definition: gis/alloc.c:149
void Rast_get_f_row(int fd, FCELL *buf, int row)
Get raster row (FCELL type)
RASTER_MAP_TYPE map_type
Definition: R.h:59
void Rast_get_null_value_row(int fd, char *flags, int row)
Read or simulate null value row.
Definition: R.h:72
char * name
Definition: R.h:63
CELL min
Definition: raster.h:37
int data_fd
Definition: R.h:68
off_t * row_ptr
Definition: R.h:49
#define min(x, y)
Definition: draw2.c:31
double DCELL
Definition: gis.h:581
CELL max
Definition: raster.h:38
void Rast_get_c_row(int fd, CELL *buf, int row)
Get raster row (CELL type)
int G_expand(unsigned char *src, int src_sz, unsigned char *dst, int dst_sz, int number)
Definition: compress.c:237
unsigned char * data
Definition: R.h:55
char * dst
Definition: lz4.h:354
void Rast_get_row(int fd, void *buf, int row, RASTER_MAP_TYPE data_type)
Get raster row.
COLUMN_MAPPING * col_map
Definition: R.h:50
void Rast__set_null_value(void *rast, int numVals, int null_is_zero, RASTER_MAP_TYPE data_type)
To set one or more raster values to null.
Definition: null_val.c:80
int Rast_is_c_null_value(const CELL *cellVal)
To check if a CELL raster value is set to NULL.
Definition: null_val.c:209
#define NULL
Definition: ccmath.h:32
unsigned char * null_bits
Definition: R.h:57
struct GDAL_link * gdal
Definition: R.h:67
#define max(x, y)
Definition: draw2.c:32
struct Reclass reclass
Definition: R.h:43
fd
Definition: d/range.c:69
CELL Rast_quant_get_cell_value(struct Quant *q, DCELL dcellVal)
Returns a CELL category for the floating-point value based on the quantization rules in q...
Definition: quant.c:601
void G_fatal_error(const char *msg,...)
Print a fatal error message to stderr.
Definition: gis/error.c:159
void Rast_get_f_row_nomask(int fd, FCELL *buf, int row)
Read raster row without masking (FCELL type)
SYMBOL * err(FILE *fp, SYMBOL *s, char *msg)
Definition: symbol/read.c:220
double C1
Definition: R.h:51
int compressed
Compression mode (raster header only)
Definition: gis.h:403
#define DCELL_TYPE
Definition: raster.h:13
size_t Rast_cell_size(RASTER_MAP_TYPE data_type)
Returns size of a raster cell in bytes.
Definition: alloc_cell.c:39
void Rast_set_c_null_value(CELL *cellVals, int numVals)
To set a number of CELL raster values to NULL.
Definition: null_val.c:124
for(cat=0;;cat++)
void Rast_zero_input_buf(void *rast, RASTER_MAP_TYPE data_type)
Definition: zero_cell.c:33
int null_fd
Definition: R.h:56
void Rast_get_d_row(int fd, DCELL *buf, int row)
Get raster row (DCELL type)
int COLUMN_MAPPING
Definition: R.h:17
void Rast_get_c_row_nomask(int fd, CELL *buf, int row)
Read raster row without masking (CELL type)
char * mapset
Definition: R.h:64
int cur_row
Definition: R.h:52
DCELL * Rast_allocate_d_input_buf(void)
Definition: alloc_cell.c:172
int Rast_is_null_value(const void *rast, RASTER_MAP_TYPE data_type)
To check if a raster value is set to NULL.
Definition: null_val.c:179
int Rast__read_null_bits(int fd, int row, unsigned char *flags)
#define check_null_bit(flags, bit_num)
int mask_fd
Definition: R.h:75
struct fileinfo * fileinfo
Definition: R.h:87
float FCELL
Definition: gis.h:582
int auto_mask
Definition: R.h:76
void Rast_get_d_row_nomask(int fd, DCELL *buf, int row)
Read raster row without masking (DCELL type)
int cur_nbytes
Definition: R.h:54
double C2
Definition: R.h:51
int cols
Number of columns for 2D data.
Definition: gis.h:409
void G_xdr_get_float(float *dst, const void *src)
Definition: gis/xdr.c:73
int null_file_exists
Definition: R.h:62
off_t * null_row_ptr
Definition: R.h:69
struct Cell_head cellhd
Definition: R.h:42
void Rast_set_c_value(void *rast, CELL cval, RASTER_MAP_TYPE data_type)
Places a CELL raster value.
Definition: raster/raster.c:96
int CELL
Definition: gis.h:580
void * G_incr_void_ptr(const void *ptr, size_t size)
Advance void pointer.
Definition: gis/alloc.c:186
Definition: R.h:39
#define _(str)
Definition: glocale.h:13
int RASTER_MAP_TYPE
Definition: raster.h:25
int null_cur_row
Definition: R.h:53
#define FCELL_TYPE
Definition: raster.h:12
struct Quant quant
Definition: R.h:66
int
Reads the categories file for map name in mapset and stores the categories in the pcats structure...
void G_zero(void *buf, int i)
Zero out a buffer, buf, of length i.
Definition: gis/zero.c:23
void G_xdr_get_double(double *dst, const void *src)
Definition: gis/xdr.c:83
int rows
Number of rows for 2D data.
Definition: gis.h:405
void G_warning(const char *msg,...)
Print a warning message to stderr.
Definition: gis/error.c:203
int Rast__null_bitstream_size(int cols)
Determines null bitstream size.
Definition: alloc_cell.c:148
double r
Definition: r_raster.c:39
CELL * table
Definition: raster.h:39
void Rast__init_null_bits(unsigned char *flags, int cols)
?
Definition: null_val.c:495
int reclass_flag
Definition: R.h:48