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