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