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
compute_window_row(int fd,int row,int * cellRow)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
do_reclass_int(int fd,void * cell,int null_is_zero)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
read_data_fp_compressed(int fd,int row,unsigned char * data_buf,int * nbytes)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
rle_decompress(unsigned char * dst,const unsigned char * src,int nbytes,int size)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
read_data_compressed(int fd,int row,unsigned char * data_buf,int * nbytes)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
read_data_uncompressed(int fd,int row,unsigned char * data_buf,int * nbytes)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
read_data_gdal(int fd,int row,unsigned char * data_buf,int * nbytes)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
read_data(int fd,int row,unsigned char * data_buf,int * nbytes)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 */
cell_values_int(int fd,const unsigned char * data,const COLUMN_MAPPING * cmap,int nbytes,void * cell,int n)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
cell_values_float(int fd,const unsigned char * data,const COLUMN_MAPPING * cmap,int nbytes,void * cell,int n)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
cell_values_double(int fd,const unsigned char * data,const COLUMN_MAPPING * cmap,int nbytes,void * cell,int n)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
gdal_values_int(int fd,const unsigned char * data,const COLUMN_MAPPING * cmap,int nbytes,CELL * cell,int n)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
gdal_values_float(int fd,const float * data,const COLUMN_MAPPING * cmap,int nbytes,FCELL * cell,int n)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
gdal_values_double(int fd,const double * data,const COLUMN_MAPPING * cmap,int nbytes,DCELL * cell,int n)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 */
transfer_to_cell_XX(int fd,void * cell)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
transfer_to_cell_fi(int fd,void * cell)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
transfer_to_cell_di(int fd,void * cell)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
transfer_to_cell_if(int fd,void * cell)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
transfer_to_cell_df(int fd,void * cell)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
transfer_to_cell_id(int fd,void * cell)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
transfer_to_cell_fd(int fd,void * cell)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 */
get_map_row_nomask(int fd,void * rast,int row,RASTER_MAP_TYPE data_type)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
get_map_row_no_reclass(int fd,void * rast,int row,RASTER_MAP_TYPE data_type,int null_is_zero,int with_mask)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
get_map_row(int fd,void * rast,int row,RASTER_MAP_TYPE data_type,int null_is_zero,int with_mask)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 */
Rast_get_row_nomask(int fd,void * buf,int row,RASTER_MAP_TYPE data_type)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 */
Rast_get_c_row_nomask(int fd,CELL * buf,int row)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 */
Rast_get_f_row_nomask(int fd,FCELL * buf,int row)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 */
Rast_get_d_row_nomask(int fd,DCELL * buf,int row)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 */
Rast_get_row(int fd,void * buf,int row,RASTER_MAP_TYPE data_type)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 */
Rast_get_c_row(int fd,CELL * buf,int row)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 */
Rast_get_f_row(int fd,FCELL * buf,int row)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 */
Rast_get_d_row(int fd,DCELL * buf,int row)820 void Rast_get_d_row(int fd, DCELL * buf, int row)
821 {
822 Rast_get_row(fd, buf, row, DCELL_TYPE);
823 }
824
read_null_bits_compressed(int null_fd,unsigned char * flags,int row,size_t size,int fd)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
Rast__read_null_bits(int fd,int row,unsigned char * flags)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
get_null_value_row_nomask(int fd,char * flags,int row)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
get_null_value_row_gdal(int fd,char * flags,int row)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
embed_mask(char * flags,int row)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
999 if (R__.fileinfo[R__.mask_fd].reclass_flag) {
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
get_null_value_row(int fd,char * flags,int row,int with_mask)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
embed_nulls(int fd,void * buf,int row,RASTER_MAP_TYPE map_type,int null_is_zero,int with_mask)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 */
Rast_get_null_value_row(int fd,char * flags,int row)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 }
1092