← Back to branch summary

~ubuntu-branches/ubuntu/wily/grass/wily

« back to all changes in this revision

Viewing changes to raster/r.viewshed/grass.cpp

Tags: 7.0.0~rc1+ds1-1~exp1
http://bugs.debian.org/776349
* New upstream release candidate.
* Repack upstream tarball, remove precompiled Python objects.
* Add upstream metadata.
* Update gbp.conf and Vcs-Git URL to use the experimental branch.
* Update watch file for GRASS 7.0.
* Drop build dependencies for Tcl/Tk, add build dependencies:
  python-numpy, libnetcdf-dev, netcdf-bin, libblas-dev, liblapack-dev
* Update Vcs-Browser URL to use cgit instead of gitweb.
* Update paths to use grass70.
* Add configure options: --with-netcdf, --with-blas, --with-lapack,
  remove --with-tcltk-includes.
* Update patches for GRASS 7.
* Update copyright file, changes:
  - Update copyright years
  - Group files by license
  - Remove unused license sections
* Add patches for various typos.
* Fix desktop file with patch instead of d/rules.
* Use minimal dh rules.
* Bump Standards-Version to 3.9.6, no changes.
* Use dpkg-maintscript-helper to replace directories with symlinks.
  (closes: #776349)
* Update my email to use @debian.org address.

Show diffs side-by-side

added added

removed removed

Lines of Context:
raster/r.viewshed/grass.cpp
 
1
 
 
2
/****************************************************************************
 
3
 *
 
4
 * MODULE:       r.viewshed
 
5
 *
 
6
 * AUTHOR(S):    Laura Toma, Bowdoin College - ltoma@bowdoin.edu
 
7
 *               Yi Zhuang - yzhuang@bowdoin.edu
 
8
 *
 
9
 *               Ported to GRASS by William Richard -
 
10
 *               wkrichar@bowdoin.edu or willster3021@gmail.com
 
11
 *               Markus Metz: surface interpolation
 
12
 *
 
13
 * Date:         april 2011 
 
14
 * 
 
15
 * PURPOSE: To calculate the viewshed (the visible cells in the
 
16
 * raster) for the given viewpoint (observer) location.  The
 
17
 * visibility model is the following: Two points in the raster are
 
18
 * considered visible to each other if the cells where they belong are
 
19
 * visible to each other.  Two cells are visible to each other if the
 
20
 * line-of-sight that connects their centers does not intersect the
 
21
 * terrain. The terrain is NOT viewed as a tesselation of flat cells, 
 
22
 * i.e. if the line-of-sight does not pass through the cell center, 
 
23
 * elevation is determined using bilinear interpolation.
 
24
 * The viewshed algorithm is efficient both in
 
25
 * terms of CPU operations and I/O operations. It has worst-case
 
26
 * complexity O(n lg n) in the RAM model and O(sort(n)) in the
 
27
 * I/O-model.  For the algorithm and all the other details see the
 
28
 * paper: "Computing Visibility on * Terrains in External Memory" by
 
29
 * Herman Haverkort, Laura Toma and Yi Zhuang.
 
30
 *
 
31
 * COPYRIGHT: (C) 2008 by the GRASS Development Team
 
32
 *
 
33
 * This program is free software under the GNU General Public License
 
34
 * (>=v2). Read the file COPYING that comes with GRASS for details.
 
35
 *
 
36
 *****************************************************************************/
 
37
 
 
38
 
 
39
#include <stdlib.h>
 
40
#include <stdio.h>
 
41
 
 
42
extern "C"
 
43
{
 
44
#include <grass/config.h>
 
45
#include <grass/gis.h>
 
46
#include <grass/glocale.h>
 
47
}
 
48
 
 
49
#include "grass.h"
 
50
#include "visibility.h"
 
51
 
 
52
 
 
53
 
 
54
 
 
55
 
 
56
/* ------------------------------------------------------------ */
 
57
/* If viewOptions.doCurv is on then adjust the passed height for
 
58
   curvature of the earth; otherwise return the passed height
 
59
   unchanged.
 
60
   If viewOptions.doRefr is on then adjust the curved height for
 
61
   the effect of atmospheric refraction too.
 
62
 */
 
63
surface_type adjust_for_curvature(Viewpoint vp, double row,
 
64
                           double col, surface_type h,
 
65
                           ViewOptions viewOptions, GridHeader *hd)
 
66
{
 
67
 
 
68
    if (!viewOptions.doCurv)
 
69
        return h;
 
70
 
 
71
    assert(viewOptions.ellps_a != 0);
 
72
 
 
73
    /* distance must be in meters because ellps_a is in meters */
 
74
    double dist = G_distance(Rast_col_to_easting(vp.col + 0.5, &(hd->window)),
 
75
                             Rast_row_to_northing(vp.row + 0.5, &(hd->window)),
 
76
                             Rast_col_to_easting(col + 0.5, &(hd->window)),
 
77
                             Rast_row_to_northing(row + 0.5, &(hd->window)));
 
78
 
 
79
    double adjustment = (dist * dist) / (2.0 * viewOptions.ellps_a);
 
80
 
 
81
    if (!viewOptions.doRefr)
 
82
        return h - adjustment;
 
83
 
 
84
    return h - (adjustment * (1.0 - viewOptions.refr_coef));
 
85
}
 
86
 
 
87
 
 
88
 
 
89
/* ************************************************************ */
 
90
/*return a GridHeader that has all the relevant data filled in */
 
91
GridHeader *read_header(char *rastName, Cell_head * region)
 
92
{
 
93
 
 
94
    assert(rastName);
 
95
 
 
96
    /*allocate the grid header to fill */
 
97
    GridHeader *hd = (GridHeader *) G_malloc(sizeof(GridHeader));
 
98
 
 
99
    assert(hd);
 
100
 
 
101
    /*get the num rows and cols from GRASS */
 
102
    int nrows, ncols;
 
103
 
 
104
    nrows = Rast_window_rows();
 
105
    ncols = Rast_window_cols();
 
106
    /*check for loss of prescion */
 
107
    if (nrows <= maxDimension && ncols <= maxDimension) {
 
108
        hd->nrows = (dimensionType) nrows;
 
109
        hd->ncols = (dimensionType) ncols;
 
110
    }
 
111
    else {
 
112
        G_warning("ERROR: nrows (%d) > maxDimension (%d) AND/OR ncols (%d) > maxDimension (%d)", 
 
113
                   nrows, maxDimension, ncols, maxDimension);
 
114
        G_fatal_error(_("Computational region too large. Use smaller area or lower raster resolution"));
 
115
    }
 
116
 
 
117
 
 
118
    /*fill in rest of header */
 
119
    hd->xllcorner = Rast_col_to_easting(0, region);
 
120
    hd->yllcorner = Rast_row_to_northing(0, region);
 
121
    /* Cell_head stores 2 resolutions, while GridHeader only stores 1 */
 
122
       // make sure the two Cell_head resolutions are equal
 
123
    if (fabs(region->ew_res - region->ns_res) > .001) {
 
124
        G_warning(_("East-west resolution does not equal north-south resolution. "
 
125
                    "The viewshed computation assumes the cells are square, so in "
 
126
                    "this case this may result in innacuracies."));
 
127
        //    exit(EXIT_FAILURE);
 
128
    }
 
129
    hd->ew_res = region->ew_res;
 
130
    hd->ns_res = region->ns_res;
 
131
    //store the null value of the map
 
132
    Rast_set_null_value(&(hd->nodata_value), 1, G_SURFACE_TYPE);
 
133
    G_verbose_message("Nodata value set to %f", hd->nodata_value);
 
134
    
 
135
    
 
136
    
 
137
    return hd;
 
138
}
 
139
 
 
140
/* calculate ENTER and EXIT event elevation (bilinear interpolation) */
 
141
surface_type calculate_event_elevation(AEvent e, int nrows, int ncols,
 
142
                                       dimensionType vprow, dimensionType vpcol,
 
143
                                       G_SURFACE_T **inrast, RASTER_MAP_TYPE data_type)
 
144
{
 
145
    int row1, col1;
 
146
    surface_type event_elev;
 
147
    G_SURFACE_T elev1, elev2, elev3, elev4;
 
148
    
 
149
    calculate_event_row_col(e, vprow, vpcol, &row1, &col1);
 
150
    if (row1 >= 0 && row1 < nrows && col1 >= 0 && col1 < ncols) {
 
151
        elev1 = inrast[row1 - e.row + 1][col1];
 
152
        elev2 = inrast[row1 - e.row + 1][e.col];
 
153
        elev3 = inrast[1][col1];
 
154
        elev4 = inrast[1][e.col];
 
155
        if (Rast_is_null_value(&elev1, data_type) ||
 
156
            Rast_is_null_value(&elev2, data_type) ||
 
157
            Rast_is_null_value(&elev3, data_type) ||
 
158
            Rast_is_null_value(&elev4, data_type))
 
159
            event_elev = inrast[1][e.col];
 
160
        else {
 
161
            event_elev = (elev1 + elev2 + elev3 + elev4) / 4.;
 
162
        }
 
163
    }
 
164
    else
 
165
        event_elev = inrast[1][e.col];
 
166
 
 
167
    return event_elev;
 
168
}
 
169
 
 
170
 
 
171
/*  ************************************************************ */
 
172
/* input: an array capable to hold the max number of events, a raster
 
173
   name, a viewpoint and the viewOptions; action: figure out all events
 
174
   in the input file, and write them to the event list. data is
 
175
   allocated and initialized with all the cells on the same row as the
 
176
   viewpoint. it returns the number of events. initialize and fill
 
177
   AEvent* with all the events for the map.  Used when solving in
 
178
   memory, so the AEvent* should fit in memory.  */
 
179
size_t
 
180
init_event_list_in_memory(AEvent * eventList, char *rastName,
 
181
                                Viewpoint * vp, GridHeader * hd,
 
182
                                ViewOptions viewOptions, surface_type ***data,
 
183
                                MemoryVisibilityGrid * visgrid)
 
184
{
 
185
 
 
186
    G_message(_("Computing events..."));
 
187
    assert(eventList && vp && visgrid);
 
188
    //GRASS should be defined 
 
189
 
 
190
    /*alloc data ; data is used to store all the cells on the same row
 
191
       as the viewpoint. */
 
192
    *data = (surface_type **)G_malloc(3 * sizeof(surface_type *));
 
193
    assert(*data);
 
194
    (*data)[0] = (surface_type *)G_malloc(3 * Rast_window_cols() * sizeof(surface_type));
 
195
    assert((*data)[0]);
 
196
    (*data)[1] = (*data)[0] + Rast_window_cols();
 
197
    (*data)[2] = (*data)[1] + Rast_window_cols();
 
198
 
 
199
    /*get the mapset name */
 
200
    const char *mapset;
 
201
 
 
202
    mapset = G_find_raster(rastName, "");
 
203
    if (mapset == NULL)
 
204
        G_fatal_error(_("Raster map [%s] not found"), rastName);
 
205
 
 
206
    /*open map */
 
207
    int infd;
 
208
 
 
209
    if ((infd = Rast_open_old(rastName, mapset)) < 0)
 
210
        G_fatal_error(_("Cannot open raster file [%s]"), rastName);
 
211
 
 
212
    /*get the data_type */
 
213
    RASTER_MAP_TYPE data_type;
 
214
 
 
215
    /* data_type = G_raster_map_type(rastName, mapset); */
 
216
    data_type = G_SURFACE_TYPE;
 
217
 
 
218
    /*buffer to hold 3 rows */
 
219
    G_SURFACE_T **inrast;
 
220
    int nrows = Rast_window_rows();
 
221
    int ncols = Rast_window_cols();
 
222
 
 
223
    inrast = (G_SURFACE_T **)G_malloc(3 * sizeof(G_SURFACE_T *));
 
224
    assert(inrast);
 
225
    inrast[0] = (G_SURFACE_T *)Rast_allocate_buf(data_type);
 
226
    assert(inrast[0]);
 
227
    inrast[1] = (G_SURFACE_T *)Rast_allocate_buf(data_type);
 
228
    assert(inrast[1]);
 
229
    inrast[2] = (G_SURFACE_T *)Rast_allocate_buf(data_type);
 
230
    assert(inrast[2]);
 
231
    
 
232
    Rast_set_null_value(inrast[0], ncols, data_type);
 
233
    Rast_set_null_value(inrast[1], ncols, data_type);
 
234
    Rast_set_null_value(inrast[2], ncols, data_type);
 
235
 
 
236
    int isnull = 0;
 
237
 
 
238
    /*keep track of the number of events added, to be returned later */
 
239
    size_t nevents = 0;
 
240
 
 
241
    /*scan through the raster data */
 
242
    dimensionType i, j;
 
243
    double ax, ay;
 
244
    AEvent e;
 
245
    
 
246
    /* read first row */
 
247
    Rast_get_row(infd, inrast[2], 0, data_type);
 
248
 
 
249
    e.angle = -1;
 
250
    for (i = 0; i < nrows; i++) {
 
251
        /*read in the raster row */
 
252
        
 
253
        G_SURFACE_T *tmprast = inrast[0];
 
254
        inrast[0] = inrast[1];
 
255
        inrast[1] = inrast[2];
 
256
        inrast[2] = tmprast;
 
257
 
 
258
        if (i < nrows - 1)
 
259
            Rast_get_row(infd, inrast[2], i + 1, data_type);
 
260
        else
 
261
            Rast_set_null_value(inrast[2], ncols, data_type);
 
262
 
 
263
        G_percent(i, nrows, 2);
 
264
 
 
265
        /*fill event list with events from this row */
 
266
        for (j = 0; j < Rast_window_cols(); j++) {
 
267
            e.row = i;
 
268
            e.col = j;
 
269
 
 
270
            /*read the elevation value into the event */
 
271
            isnull = Rast_is_null_value(&(inrast[1][j]), data_type);
 
272
            e.elev[1] = inrast[1][j];
 
273
 
 
274
            /* adjust for curvature */
 
275
            e.elev[1] = adjust_for_curvature(*vp, i, j, e.elev[1], viewOptions, hd);
 
276
 
 
277
            /*write it into the row of data going through the viewpoint */
 
278
            if (i == vp->row) {
 
279
                (*data)[0][j] = e.elev[1];
 
280
                (*data)[1][j] = e.elev[1];
 
281
                (*data)[2][j] = e.elev[1];
 
282
            }
 
283
 
 
284
            /* set the viewpoint, and don't insert it into eventlist */
 
285
            if (i == vp->row && j == vp->col) {
 
286
                set_viewpoint_elev(vp, e.elev[1] + viewOptions.obsElev);
 
287
                if (viewOptions.tgtElev > 0)
 
288
                    vp->target_offset = viewOptions.tgtElev;
 
289
                else
 
290
                    vp->target_offset = 0.;
 
291
                if (isnull) {
 
292
                    /*what to do when viewpoint is NODATA ? */
 
293
                    G_warning(_("Viewpoint is NODATA."));
 
294
                    G_message(_("Will assume its elevation is = %f"),
 
295
                              vp->elev);
 
296
                }
 
297
 
 
298
                add_result_to_inmem_visibilitygrid(visgrid, i, j,
 
299
                                                   180);
 
300
                continue;
 
301
            }
 
302
 
 
303
            /*don't insert in eventlist nodata cell events */
 
304
            if (isnull) {
 
305
                /* record this cell as being NODATA; this is necessary so
 
306
                   that we can distingush invisible events, from nodata
 
307
                   events in the output */
 
308
                add_result_to_inmem_visibilitygrid(visgrid, i, j,
 
309
                                                   hd->nodata_value);
 
310
                continue;
 
311
            }
 
312
 
 
313
            /* if point is outside maxDist, do NOT include it as an
 
314
               event */
 
315
            if (is_point_outside_max_dist
 
316
                (*vp, *hd, i, j, viewOptions.maxDist))
 
317
                continue;
 
318
 
 
319
            /* if it got here it is not the viewpoint, not NODATA, and
 
320
               within max distance from viewpoint; generate its 3 events
 
321
               and insert them */
 
322
 
 
323
            /* get ENTER elevation */
 
324
            e.eventType = ENTERING_EVENT;
 
325
            e.elev[0] = calculate_event_elevation(e, nrows, ncols,
 
326
                                       vp->row, vp->col, inrast, data_type);
 
327
            /* adjust for curvature */
 
328
            if (viewOptions.doCurv) {
 
329
                calculate_event_position(e, vp->row, vp->col, &ay, &ax);
 
330
                e.elev[0] = adjust_for_curvature(*vp, ay, ax, e.elev[0], viewOptions, hd);
 
331
            }
 
332
 
 
333
            /* get EXIT elevation */
 
334
            e.eventType = EXITING_EVENT;
 
335
            e.elev[2] = calculate_event_elevation(e, nrows, ncols,
 
336
                                       vp->row, vp->col, inrast, data_type);
 
337
            /* adjust for curvature */
 
338
            if (viewOptions.doCurv) {
 
339
                calculate_event_position(e, vp->row, vp->col, &ay, &ax);
 
340
                e.elev[2] = adjust_for_curvature(*vp, ay, ax, e.elev[2], viewOptions, hd);
 
341
            }
 
342
 
 
343
            /*write adjusted elevation into the row of data going through the viewpoint */
 
344
            if (i == vp->row) {
 
345
                (*data)[0][j] = e.elev[0];
 
346
                (*data)[1][j] = e.elev[1];
 
347
                (*data)[2][j] = e.elev[2];
 
348
            }
 
349
 
 
350
            /*put event into event list */
 
351
            e.eventType = ENTERING_EVENT;
 
352
            calculate_event_position(e, vp->row, vp->col, &ay, &ax);
 
353
            e.angle = calculate_angle(ax, ay, vp->col, vp->row);
 
354
            eventList[nevents] = e;
 
355
            nevents++;
 
356
 
 
357
            e.eventType = CENTER_EVENT;
 
358
            calculate_event_position(e, vp->row, vp->col, &ay, &ax);
 
359
            e.angle = calculate_angle(ax, ay, vp->col, vp->row);
 
360
            eventList[nevents] = e;
 
361
            nevents++;
 
362
 
 
363
            e.eventType = EXITING_EVENT;
 
364
            calculate_event_position(e, vp->row, vp->col, &ay, &ax);
 
365
            e.angle = calculate_angle(ax, ay, vp->col, vp->row);
 
366
            eventList[nevents] = e;
 
367
            nevents++;
 
368
 
 
369
        }
 
370
    }
 
371
    G_percent(nrows, nrows, 2);
 
372
 
 
373
    Rast_close(infd);
 
374
 
 
375
    G_free(inrast[0]);
 
376
    G_free(inrast[1]);
 
377
    G_free(inrast[2]);
 
378
    G_free(inrast);
 
379
 
 
380
    return nevents;
 
381
}
 
382
 
 
383
 
 
384
 
 
385
 
 
386
 
 
387
/* ************************************************************ */
 
388
/* input: an arcascii file, a grid header and a viewpoint; action:
 
389
   figure out all events in the input file, and write them to the
 
390
   stream.  It assumes the file pointer is positioned rigth after the
 
391
   grid header so that this function can read all grid elements.
 
392
 
 
393
   if data is not NULL, it creates an array that stores all events on
 
394
   the same row as the viewpoint. 
 
395
 */
 
396
AMI_STREAM < AEvent > *init_event_list(char *rastName, Viewpoint * vp,
 
397
                                             GridHeader * hd,
 
398
                                             ViewOptions viewOptions,
 
399
                                             surface_type ***data,
 
400
                                             IOVisibilityGrid * visgrid)
 
401
{
 
402
 
 
403
    G_message(_("Computing events..."));
 
404
    assert(rastName && vp && hd && visgrid);
 
405
 
 
406
    if (data != NULL) {
 
407
        /*data is used to store all the cells on the same row as the
 
408
           //viewpoint. */
 
409
        *data = (surface_type **)G_malloc(3 * sizeof(surface_type *));
 
410
        assert(*data);
 
411
        (*data)[0] = (surface_type *)G_malloc(3 * Rast_window_cols() * sizeof(surface_type));
 
412
        assert((*data)[0]);
 
413
        (*data)[1] = (*data)[0] + Rast_window_cols();
 
414
        (*data)[2] = (*data)[1] + Rast_window_cols();
 
415
    }
 
416
 
 
417
    /*create the event stream that will hold the events */
 
418
    AMI_STREAM < AEvent > *eventList = new AMI_STREAM < AEvent > ();
 
419
    assert(eventList);
 
420
 
 
421
    /*determine which mapset we are in */
 
422
    const char *mapset;
 
423
 
 
424
    mapset = G_find_raster(rastName, "");
 
425
    if (mapset == NULL)
 
426
        G_fatal_error(_("Raster map [%s] not found"), rastName);
 
427
 
 
428
    /*open map */
 
429
    int infd;
 
430
 
 
431
    if ((infd = Rast_open_old(rastName, mapset)) < 0)
 
432
        G_fatal_error(_("Cannot open raster file [%s]"), rastName);
 
433
 
 
434
    RASTER_MAP_TYPE data_type;
 
435
 
 
436
    /* data_type = G_raster_map_type(rastName, mapset); */
 
437
    data_type = G_SURFACE_TYPE;
 
438
    G_SURFACE_T **inrast;
 
439
    int nrows = Rast_window_rows();
 
440
    int ncols = Rast_window_cols();
 
441
 
 
442
    inrast = (G_SURFACE_T **)G_malloc(3 * sizeof(G_SURFACE_T *));
 
443
    assert(inrast);
 
444
    inrast[0] = (G_SURFACE_T *)Rast_allocate_buf(data_type);
 
445
    assert(inrast[0]);
 
446
    inrast[1] = (G_SURFACE_T *)Rast_allocate_buf(data_type);
 
447
    assert(inrast[1]);
 
448
    inrast[2] = (G_SURFACE_T *)Rast_allocate_buf(data_type);
 
449
    assert(inrast[2]);
 
450
    
 
451
    Rast_set_null_value(inrast[0], ncols, data_type);
 
452
    Rast_set_null_value(inrast[1], ncols, data_type);
 
453
    Rast_set_null_value(inrast[2], ncols, data_type);
 
454
 
 
455
    /*scan through the raster data */
 
456
    int isnull = 0;
 
457
    dimensionType i, j;
 
458
    double ax, ay;
 
459
    AEvent e;
 
460
 
 
461
    /* read first row */
 
462
    Rast_get_row(infd, inrast[2], 0, data_type);
 
463
 
 
464
    e.angle = -1;
 
465
 
 
466
    /*start scanning through the grid */
 
467
    for (i = 0; i < nrows; i++) {
 
468
        
 
469
        G_percent(i, nrows, 2);
 
470
 
 
471
        /*read in the raster row */
 
472
        
 
473
        G_SURFACE_T *tmprast = inrast[0];
 
474
        inrast[0] = inrast[1];
 
475
        inrast[1] = inrast[2];
 
476
        inrast[2] = tmprast;
 
477
 
 
478
        if (i < nrows - 1)
 
479
            Rast_get_row(infd, inrast[2], i + 1, data_type);
 
480
        else
 
481
            Rast_set_null_value(inrast[2], ncols, data_type);
 
482
 
 
483
        /*fill event list with events from this row */
 
484
        for (j = 0; j < ncols; j++) {
 
485
 
 
486
            e.row = i;
 
487
            e.col = j;
 
488
 
 
489
            /*read the elevation value into the event */
 
490
            isnull = Rast_is_null_value(&(inrast[1][j]), data_type);
 
491
            e.elev[1] = inrast[1][j];
 
492
 
 
493
            /* adjust for curvature */
 
494
            e.elev[1] = adjust_for_curvature(*vp, i, j, e.elev[1], viewOptions, hd);
 
495
 
 
496
            if (data != NULL) {
 
497
 
 
498
                /*write the row of data going through the viewpoint */
 
499
                if (i == vp->row) {
 
500
                    (*data)[0][j] = e.elev[1];
 
501
                    (*data)[1][j] = e.elev[1];
 
502
                    (*data)[2][j] = e.elev[1];
 
503
                }
 
504
            }
 
505
 
 
506
            /* set the viewpoint */
 
507
            if (i == vp->row && j == vp->col) {
 
508
                set_viewpoint_elev(vp, e.elev[1] + viewOptions.obsElev);
 
509
                /*what to do when viewpoint is NODATA */
 
510
                if (is_nodata(hd, e.elev[1])) {
 
511
                    G_warning("Viewpoint is NODATA.");
 
512
                    G_message("Will assume its elevation is %.f", e.elev[1]);
 
513
                };
 
514
                if (viewOptions.tgtElev > 0)
 
515
                    vp->target_offset = viewOptions.tgtElev;
 
516
                else
 
517
                    vp->target_offset = 0.;
 
518
 
 
519
                /* add viewpoint to visibility grid */
 
520
                VisCell visCell = { i, j, 180 };
 
521
                add_result_to_io_visibilitygrid(visgrid, &visCell);
 
522
 
 
523
                /*don't insert viewpoint into eventlist */
 
524
                continue;
 
525
            }
 
526
 
 
527
            /*don't insert the nodata cell events */
 
528
            if (is_nodata(hd, e.elev[1])) {
 
529
                /* record this cell as being NODATA. ; this is necessary so
 
530
                   that we can distingush invisible events, from nodata
 
531
                   events in the output */
 
532
                VisCell visCell = { i, j, hd->nodata_value };
 
533
                add_result_to_io_visibilitygrid(visgrid, &visCell);
 
534
                continue;
 
535
            }
 
536
 
 
537
            /* if point is outside maxDist, do NOT include it as an
 
538
               event */
 
539
            if (is_point_outside_max_dist
 
540
                (*vp, *hd, i, j, viewOptions.maxDist))
 
541
                continue;
 
542
 
 
543
            /* get ENTER elevation */
 
544
            e.eventType = ENTERING_EVENT;
 
545
            e.elev[0] = calculate_event_elevation(e, nrows, ncols,
 
546
                                       vp->row, vp->col, inrast, data_type);
 
547
            /* adjust for curvature */
 
548
            if (viewOptions.doCurv) {
 
549
                calculate_event_position(e, vp->row, vp->col, &ay, &ax);
 
550
                e.elev[0] = adjust_for_curvature(*vp, ay, ax, e.elev[0], viewOptions, hd);
 
551
            }
 
552
 
 
553
            /* get EXIT elevation */
 
554
            e.eventType = EXITING_EVENT;
 
555
            e.elev[2] = calculate_event_elevation(e, nrows, ncols,
 
556
                                       vp->row, vp->col, inrast, data_type);
 
557
            /* adjust for curvature */
 
558
            if (viewOptions.doCurv) {
 
559
                calculate_event_position(e, vp->row, vp->col, &ay, &ax);
 
560
                e.elev[2] = adjust_for_curvature(*vp, ay, ax, e.elev[2], viewOptions, hd);
 
561
            }
 
562
 
 
563
            if (data != NULL) {
 
564
 
 
565
                /*write the row of data going through the viewpoint */
 
566
                if (i == vp->row) {
 
567
                    (*data)[0][j] = e.elev[0];
 
568
                    (*data)[1][j] = e.elev[1];
 
569
                    (*data)[2][j] = e.elev[2];
 
570
                }
 
571
            }
 
572
 
 
573
            /*put event into event list */
 
574
            e.eventType = ENTERING_EVENT;
 
575
            calculate_event_position(e, vp->row, vp->col, &ay, &ax);
 
576
            e.angle = calculate_angle(ax, ay, vp->col, vp->row);
 
577
            eventList->write_item(e);
 
578
 
 
579
            e.eventType = CENTER_EVENT;
 
580
            calculate_event_position(e, vp->row, vp->col, &ay, &ax);
 
581
            e.angle = calculate_angle(ax, ay, vp->col, vp->row);
 
582
            eventList->write_item(e);
 
583
 
 
584
            e.eventType = EXITING_EVENT;
 
585
            calculate_event_position(e, vp->row, vp->col, &ay, &ax);
 
586
            e.angle = calculate_angle(ax, ay, vp->col, vp->row);
 
587
            eventList->write_item(e);
 
588
        }                       /* for j */
 
589
 
 
590
    }                           /* for i */
 
591
    G_percent(nrows, nrows, 2);
 
592
 
 
593
    Rast_close(infd);
 
594
 
 
595
    G_free(inrast[0]);
 
596
    G_free(inrast[1]);
 
597
    G_free(inrast[2]);
 
598
    G_free(inrast);
 
599
 
 
600
    G_debug(1, "nbEvents = %lu", (unsigned long)eventList->stream_len());
 
601
    G_debug(1, "Event stream length: %lu x %dB (%lu MB)",
 
602
              (unsigned long)eventList->stream_len(), (int)sizeof(AEvent),
 
603
              (unsigned long)(((long long)(eventList->stream_len() *
 
604
                                           sizeof(AEvent))) >> 20));
 
605
 
 
606
    return eventList;
 
607
}
 
608
 
 
609
 
 
610
 
 
611
 
 
612
 
 
613
 
 
614
/* ************************************************************ */
 
615
/*  saves the grid into a GRASS raster.  Loops through all elements x
 
616
   in row-column order and writes fun(x) to file. */
 
617
void
 
618
save_grid_to_GRASS(Grid * grid, char *filename, RASTER_MAP_TYPE type,
 
619
                   OutputMode mode)
 
620
{
 
621
 
 
622
    G_important_message(_("Writing output raster map..."));
 
623
    assert(grid && filename);
 
624
 
 
625
    /*open the new raster  */
 
626
    int outfd;
 
627
 
 
628
    outfd = Rast_open_new(filename, type);
 
629
 
 
630
    /*get the buffer to store values to read and write to each row */
 
631
    void *outrast;
 
632
 
 
633
    outrast = Rast_allocate_buf(type);
 
634
    assert(outrast);
 
635
 
 
636
    dimensionType i, j;
 
637
 
 
638
    for (i = 0; i < Rast_window_rows(); i++) {
 
639
        G_percent(i, Rast_window_rows(), 5);
 
640
        for (j = 0; j < Rast_window_cols(); j++) {
 
641
            if (is_invisible_nodata(grid->grid_data[i][j])) {
 
642
                writeNodataValue(outrast, j, type);
 
643
            }
 
644
            else if (mode == OUTPUT_BOOL) {
 
645
                ((CELL *) outrast)[j] = (CELL) booleanVisibilityOutput(grid->grid_data[i][j]);
 
646
            }
 
647
            else if (mode == OUTPUT_ANGLE) {
 
648
                if (is_visible(grid->grid_data[i][j])) {
 
649
                    ((FCELL *) outrast)[j] = (FCELL) angleVisibilityOutput(grid->grid_data[i][j]);
 
650
                }
 
651
                else {
 
652
                    writeNodataValue(outrast, j, FCELL_TYPE);
 
653
                }
 
654
            }
 
655
        }                       /* for j */
 
656
        Rast_put_row(outfd, outrast, type);
 
657
    }                           /* for i */
 
658
    G_percent(1, 1, 1);
 
659
 
 
660
    G_free(outrast);
 
661
    Rast_close(outfd);
 
662
    return;
 
663
}
 
664
 
 
665
 
 
666
 
 
667
 
 
668
 
 
669
/* ************************************************************ */
 
670
/*  using the visibility information recorded in visgrid, it creates an
 
671
   output viewshed raster with name outfname; for every point p that
 
672
   is visible in the grid, the corresponding value in the output
 
673
   raster is elevation(p) - viewpoint_elevation(p); the elevation
 
674
   values are read from elevfname raster */
 
675
 
 
676
void
 
677
save_vis_elev_to_GRASS(Grid * visgrid, char *elevfname, char *visfname,
 
678
                       float vp_elev)
 
679
{
 
680
 
 
681
    G_message(_("Saving grid to <%s>"), visfname);
 
682
    assert(visgrid && elevfname && visfname);
 
683
 
 
684
    /*get the mapset name */
 
685
    const char *mapset;
 
686
 
 
687
    mapset = G_find_raster(elevfname, "");
 
688
    if (mapset == NULL)
 
689
        G_fatal_error(_("Raster map [%s] not found"), elevfname);
 
690
 
 
691
    /*open elevation map */
 
692
    int elevfd;
 
693
 
 
694
    if ((elevfd = Rast_open_old(elevfname, mapset)) < 0)
 
695
        G_fatal_error(_("Cannot open raster file [%s]"), elevfname);
 
696
 
 
697
    /*get elevation data_type */
 
698
    RASTER_MAP_TYPE elev_data_type;
 
699
 
 
700
    elev_data_type = Rast_map_type(elevfname, mapset);
 
701
 
 
702
    /* create the visibility raster of same type */
 
703
    int visfd;
 
704
 
 
705
    visfd = Rast_open_new(visfname, elev_data_type);
 
706
 
 
707
    /*get the buffers to store each row */
 
708
    void *elevrast;
 
709
 
 
710
    elevrast = Rast_allocate_buf(elev_data_type);
 
711
    assert(elevrast);
 
712
    void *visrast;
 
713
 
 
714
    visrast = Rast_allocate_buf(elev_data_type);
 
715
    assert(visrast);
 
716
 
 
717
    dimensionType i, j;
 
718
    double elev = 0, viewshed_value;
 
719
 
 
720
    for (i = 0; i < Rast_window_rows(); i++) {
 
721
        /* get the row from elevation */
 
722
        Rast_get_row(elevfd, elevrast, i, elev_data_type);
 
723
 
 
724
        for (j = 0; j < Rast_window_cols(); j++) {
 
725
 
 
726
            /* read the current elevation value */
 
727
            int isNull = 0;
 
728
 
 
729
            switch (elev_data_type) {
 
730
            case CELL_TYPE:
 
731
                isNull = Rast_is_c_null_value(&((CELL *) elevrast)[j]);
 
732
                elev = (double)(((CELL *) elevrast)[j]);
 
733
                break;
 
734
            case FCELL_TYPE:
 
735
                isNull = Rast_is_f_null_value(&((FCELL *) elevrast)[j]);
 
736
                elev = (double)(((FCELL *) elevrast)[j]);
 
737
                break;
 
738
            case DCELL_TYPE:
 
739
                isNull = Rast_is_d_null_value(&((DCELL *) elevrast)[j]);
 
740
                elev = (double)(((DCELL *) elevrast)[j]);
 
741
                break;
 
742
            }
 
743
 
 
744
            if (is_visible(visgrid->grid_data[i][j])) {
 
745
                /* elevation cannot be null */
 
746
                assert(!isNull);
 
747
                /* write elev - viewpoint_elevation */
 
748
                viewshed_value = elev - vp_elev;
 
749
                writeValue(visrast, j, viewshed_value, elev_data_type);
 
750
            }
 
751
            else if (is_invisible_not_nodata(visgrid->grid_data[i][j])) {
 
752
                /* elevation cannot be null */
 
753
                assert(!isNull);
 
754
                /* write INVISIBLE */
 
755
                /*
 
756
                viewshed_value = INVISIBLE;
 
757
                writeValue(visrast, j, viewshed_value, elev_data_type);
 
758
                */
 
759
                /* write  NODATA */
 
760
                writeNodataValue(visrast, j, elev_data_type);
 
761
            }
 
762
            else {
 
763
                /* nodata */
 
764
                assert(isNull);
 
765
                /* write  NODATA */
 
766
                writeNodataValue(visrast, j, elev_data_type);
 
767
            }
 
768
 
 
769
 
 
770
        }                       /* for j */
 
771
        Rast_put_row(visfd, visrast, elev_data_type);
 
772
    }                           /* for i */
 
773
 
 
774
    Rast_close(elevfd);
 
775
    Rast_close(visfd);
 
776
    return;
 
777
}
 
778
 
 
779
 
 
780
 
 
781
 
 
782
/* helper function to deal with GRASS writing to a row buffer */
 
783
void writeValue(void *bufrast, int j, double x, RASTER_MAP_TYPE data_type)
 
784
{
 
785
 
 
786
    switch (data_type) {
 
787
    case CELL_TYPE:
 
788
        ((CELL *) bufrast)[j] = (CELL) x;
 
789
        break;
 
790
    case FCELL_TYPE:
 
791
        ((FCELL *) bufrast)[j] = (FCELL) x;
 
792
        break;
 
793
    case DCELL_TYPE:
 
794
        ((DCELL *) bufrast)[j] = (DCELL) x;
 
795
        break;
 
796
    default:
 
797
        G_fatal_error(_("Unknown data type"));
 
798
    }
 
799
}
 
800
 
 
801
void writeNodataValue(void *bufrast, int j, RASTER_MAP_TYPE data_type)
 
802
{
 
803
 
 
804
    switch (data_type) {
 
805
    case CELL_TYPE:
 
806
        Rast_set_c_null_value(&((CELL *) bufrast)[j], 1);
 
807
        break;
 
808
    case FCELL_TYPE:
 
809
        Rast_set_f_null_value(&((FCELL *) bufrast)[j], 1);
 
810
        break;
 
811
    case DCELL_TYPE:
 
812
        Rast_set_d_null_value(&((DCELL *) bufrast)[j], 1);
 
813
        break;
 
814
    default:
 
815
        G_fatal_error(_("Unknown data type"));
 
816
    }
 
817
}
 
818
 
 
819
 
 
820
/* ************************************************************ */
 
821
/* write the visibility grid to GRASS. assume all cells that are not
 
822
   in stream are NOT visible. assume stream is sorted in (i,j) order.
 
823
   for each value x it writes to grass fun(x) */
 
824
void
 
825
save_io_visibilitygrid_to_GRASS(IOVisibilityGrid * visgrid,
 
826
                                char *fname, RASTER_MAP_TYPE type,
 
827
                                float (*fun) (float))
 
828
{
 
829
 
 
830
    G_message(_("Saving grid to <%s>"), fname);
 
831
    assert(fname && visgrid);
 
832
 
 
833
    /* open the output raster  and set up its row buffer */
 
834
    int visfd;
 
835
 
 
836
    visfd = Rast_open_new(fname, type);
 
837
    void *visrast;
 
838
 
 
839
    visrast = Rast_allocate_buf(type);
 
840
    assert(visrast);
 
841
 
 
842
    /*set up reading data from visibility stream */
 
843
    AMI_STREAM < VisCell > *vstr = visgrid->visStr;
 
844
    off_t streamLen, counter = 0;
 
845
 
 
846
    streamLen = vstr->stream_len();
 
847
    vstr->seek(0);
 
848
 
 
849
    /*read the first element */
 
850
    AMI_err ae;
 
851
    VisCell *curResult = NULL;
 
852
 
 
853
    if (streamLen > 0) {
 
854
        ae = vstr->read_item(&curResult);
 
855
        assert(ae == AMI_ERROR_NO_ERROR);
 
856
        counter++;
 
857
    }
 
858
 
 
859
    dimensionType i, j;
 
860
 
 
861
    for (i = 0; i < Rast_window_rows(); i++) {
 
862
        for (j = 0; j < Rast_window_cols(); j++) {
 
863
 
 
864
            if (curResult->row == i && curResult->col == j) {
 
865
                /*cell is recodred in the visibility stream: it must be
 
866
                   either visible, or NODATA  */
 
867
                if (is_visible(curResult->angle))
 
868
                    writeValue(visrast, j, fun(curResult->angle), type);
 
869
                else
 
870
                    writeNodataValue(visrast, j, type);
 
871
 
 
872
                /*read next element of stream */
 
873
                if (counter < streamLen) {
 
874
                    ae = vstr->read_item(&curResult);
 
875
                    assert(ae == AMI_ERROR_NO_ERROR);
 
876
                    counter++;
 
877
                }
 
878
            }
 
879
            else {
 
880
                /*  this cell is not in stream, so it is invisible */
 
881
                writeNodataValue(visrast, j, type);
 
882
            }
 
883
        }                       /* for j */
 
884
 
 
885
        Rast_put_row(visfd, visrast, type);
 
886
    }                           /* for i */
 
887
 
 
888
    Rast_close(visfd);
 
889
}
 
890
 
 
891
 
 
892
 
 
893
 
 
894
 
 
895
 
 
896
 
 
897
/* ************************************************************ */
 
898
/*  using the visibility information recorded in visgrid, it creates
 
899
   an output viewshed raster with name outfname; for every point p
 
900
   that is visible in the grid, the corresponding value in the output
 
901
   raster is elevation(p) - viewpoint_elevation(p); the elevation
 
902
   values are read from elevfname raster. assume stream is sorted in
 
903
   (i,j) order. */
 
904
void
 
905
save_io_vis_and_elev_to_GRASS(IOVisibilityGrid * visgrid, char *elevfname,
 
906
                              char *visfname, float vp_elev)
 
907
{
 
908
 
 
909
    G_message(_("Saving grid to <%s>"), visfname);
 
910
    assert(visfname && visgrid);
 
911
 
 
912
    /*get mapset name and data type */
 
913
    const char *mapset;
 
914
 
 
915
    mapset = G_find_raster(elevfname, "");
 
916
    if (mapset == NULL)
 
917
        G_fatal_error(_("Opening <%s>: cannot find raster"), elevfname);
 
918
 
 
919
    /*open elevation map */
 
920
    int elevfd;
 
921
 
 
922
    if ((elevfd = Rast_open_old(elevfname, mapset)) < 0)
 
923
        G_fatal_error(_("Cannot open raster file [%s]"), elevfname);
 
924
 
 
925
    /*get elevation data_type */
 
926
    RASTER_MAP_TYPE elev_data_type;
 
927
 
 
928
    elev_data_type = Rast_map_type(elevfname, mapset);
 
929
 
 
930
    /* open visibility raster */
 
931
    int visfd;
 
932
 
 
933
    visfd = Rast_open_new(visfname, elev_data_type);
 
934
 
 
935
    /*get the buffers to store each row */
 
936
    void *elevrast;
 
937
 
 
938
    elevrast = Rast_allocate_buf(elev_data_type);
 
939
    assert(elevrast);
 
940
    void *visrast;
 
941
 
 
942
    visrast = Rast_allocate_buf(elev_data_type);
 
943
    assert(visrast);
 
944
 
 
945
    /*set up stream reading stuff */
 
946
    off_t streamLen, counter = 0;
 
947
    AMI_err ae;
 
948
    VisCell *curResult = NULL;
 
949
 
 
950
    AMI_STREAM < VisCell > *vstr = visgrid->visStr;
 
951
    streamLen = vstr->stream_len();
 
952
    vstr->seek(0);
 
953
 
 
954
    /*read the first element */
 
955
    if (streamLen > 0) {
 
956
        ae = vstr->read_item(&curResult);
 
957
        assert(ae == AMI_ERROR_NO_ERROR);
 
958
        counter++;
 
959
    }
 
960
 
 
961
    dimensionType i, j;
 
962
    double elev = 0, viewshed_value;
 
963
 
 
964
    for (i = 0; i < Rast_window_rows(); i++) {
 
965
 
 
966
        Rast_get_row(elevfd, elevrast, i, elev_data_type);
 
967
 
 
968
        for (j = 0; j < Rast_window_cols(); j++) {
 
969
 
 
970
            /* read the current elevation value */
 
971
            int isNull = 0;
 
972
 
 
973
            switch (elev_data_type) {
 
974
            case CELL_TYPE:
 
975
                isNull = Rast_is_c_null_value(&((CELL *) elevrast)[j]);
 
976
                elev = (double)(((CELL *) elevrast)[j]);
 
977
                break;
 
978
            case FCELL_TYPE:
 
979
                isNull = Rast_is_f_null_value(&((FCELL *) elevrast)[j]);
 
980
                elev = (double)(((FCELL *) elevrast)[j]);
 
981
                break;
 
982
            case DCELL_TYPE:
 
983
                isNull = Rast_is_d_null_value(&((DCELL *) elevrast)[j]);
 
984
                elev = (double)(((DCELL *) elevrast)[j]);
 
985
                break;
 
986
            }
 
987
 
 
988
 
 
989
            if (curResult->row == i && curResult->col == j) {
 
990
                /*cell is recodred in the visibility stream: it must be
 
991
                   either visible, or NODATA  */
 
992
                if (is_visible(curResult->angle))
 
993
                    writeValue(visrast, j, elev - vp_elev, elev_data_type);
 
994
                else
 
995
                    writeNodataValue(visrast, j, elev_data_type);
 
996
                /*read next element of stream */
 
997
                if (counter < streamLen) {
 
998
                    ae = vstr->read_item(&curResult);
 
999
                    assert(ae == AMI_ERROR_NO_ERROR);
 
1000
                    counter++;
 
1001
                }
 
1002
            }
 
1003
            else {
 
1004
                /*  this cell is not in stream, so it is  invisible */
 
1005
                    writeNodataValue(visrast, j, elev_data_type);
 
1006
            }
 
1007
        }                       /* for j */
 
1008
 
 
1009
        Rast_put_row(visfd, visrast, elev_data_type);
 
1010
    }                           /* for i */
 
1011
 
 
1012
    Rast_close(elevfd);
 
1013
    Rast_close(visfd);
 
1014
    return;
 
1015
}