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