r.viewshed precision bug

[neteler]

I have added a test script for r.viewshed which illustrates the precision bug:

grassaddons/raster/r.viewshed/testscript.sh

Attached related screenshots. This is a showstopper - util fixed, r.los should not be replaced.

Markus

[neteler]

r.los correct output

[neteler]

r.viewshed incorrect output on same synthetic DEM

[hamish]

Hi,

slightly hijacking the bug report, but fyi for "areas outside of LOS should be NULL, not zero", I fixed that for r.los some months ago; maybe the code is related & the patch near the same? (no idea, but it would be worth a quick look to check)

Hamish

[mmetz]

The precision bug is hopefully fixed in r46199.

Markus M

[mmetz]

Replying to hamish:

Hi, slightly hijacking the bug report, but fyi for "areas outside of LOS should be NULL, not zero", I fixed that for r.los some months ago; maybe the code is related & the patch near the same? (no idea, but it would be worth a quick look to check)

To what kind of output does this "NULL, not zero output" refer to? Binary, i.e. NULL = not visible, 1 = visible instead of zero = not visible, 1 = visible? Or viewing angle or elevation difference?

[hamish]

The precision bug is hopefully fixed in r46199.

excellent...! with that fixed, how close do you think it is to migrating r.viewshed out of addons and into the main source tree? or is there still work to do?

if it's ready, does it mean r.los is now redundant, or does r.los have some purpose/feature worth saving?

To what kind of output does this "NULL, not zero output" refer to?

I think I meant visible/hidden (see r25102), but in hindsight I suspect that you can completely disregard comment:1.

thanks, Hamish

[mmetz]

Replying to hamish:

The precision bug is hopefully fixed in r46199.

excellent...! with that fixed, how close do you think it is to migrating r.viewshed out of addons and into the main source tree? or is there still work to do? if it's ready, does it mean r.los is now redundant, or does r.los have some purpose/feature worth saving?

Maybe the r.los input option patt_map which is used as a mask, but from reading the r.los manual, the difference to setting a MASK is not clear to me.

For GRASS 7, we are free to change module options, for GRASS 6, this patt_map option could be added for compatibility reasons, but the only effect would then be to print out a warning that a MASK should be used instead of patt_map=mask_map?

[hamish]

note compiler warning: (x3)

include/grass/iostream/replacementHeapBlock.h:146: warning: 'str' may be used uninitialized in this function

spearfish test:

g.region -d
r.surf.volcano output=gauss method=gaussian sigma=1
d.where
r.viewshed in=gauss out=gauss.shed coord=597900,4920840 max_dist=7500 --o

###
d.histogram gauss.shed

seems to keep -1 as invisibile not NULL ????? fix:

r.null gauss.shed setnull=-1
r.colors gauss.shed color=bcyr
d.redraw

# 180deg is just @ viewing cell, so for good colors we also need:  (PITA, so maybe we could build that into the module)
r.null gauss.shed setnull=180
r.colors gauss.shed color=bcyr
d.redraw
###

echo "597900,4920840" | v.in.ascii out=view_pt fs=,
nviz gauss color=gauss.shed points=view_pt

coord, -r desc: lat/lon -> easting/northing

(seems to use current e,n LOCATION coords already..)

drop the -r flag now that debugging stage is ending?

tgt_elev=: how is that used if a target coord is never defined?

any interesting coverage/patch map functionality from r.cva which would be interesting to clean-room reimpliement?

timings on my trusty old P4:

G:65> time r.viewshed in=elevation.dem out=elev.shed coord=597900,4920840 max_dist=7500 --o
real    0m3.898s
user    0m3.768s
sys     0m0.084s

G65> time r.los in=elevation.dem out=elev.shed coord=597900,4920840 max_dist=7500
real    0m4.156s
user    0m4.004s
sys     0m0.108s


G65> g.region rast=elevation.10m
G65> time r.viewshed in=elevation.10m out=elev.shed coord=597900,4920840 max_dist=7500 --o
real    0m42.330s
user    0m40.791s
sys     0m0.532s

G65> time r.los in=elevation.10m out=elev.shed coord=597900,4920840 max_dist=7500 --o
real    (finsihed in about 20 minutes)

indeed it seems to scale a bit better than r.los ... :-)

For GRASS 7, we are free to change module options, for GRASS 6, this patt_map option could be added for compatibility reasons, but the only effect would then be to print out a warning that a MASK should be used instead of patt_map=mask_map?

I wasn't thinking to make this as r.los2; as it is a totally diff't code base I'd add it named as r.viewshed, and drop r.los in grass7. but maybe it is prefered to keep the r.los name? if so then at least in GRASS 6.x it would be fine to have the option as "This does nothing and is kept for backwards compatibility reasons." and the warning as you describe. ? since r.los is so useless at modern region array sizes, maybe we should.. alternately, we could add r.viewshed in 6.x with a wrapper script for the old r.los name, like r.cats -> r.category.

Hamish

[hamish]

Hamish:

tgt_elev=: how is that used if a target coord is never defined?

never mind, I get it.

alternately, we could add r.viewshed in 6.x with a wrapper script for the old r.los name, like r.cats -> r.category.

I'm liking this idea. [+ compatibility option for gr6]

any interesting coverage/patch map functionality from r.cva which would be interesting to clean-room reimplemented?

fyi, r.cva:

  r.cva draws heavily on the code for r.los, which was written by
  Kewan Q. Khawaja, Intelligent Engineering Systems Laboratory, M.I.T.

  *****

  PURPOSE

  1) Perform cumulative viewshed (CVA) analysis.

  OPTIONS

  1)  Record the number of cells visible from:
  
  all          every map cell;
  systematic   a systematic sample of map cells;
  random       a random sample of map cells;
  sites        a list of sites.

  2)  Record how often each cell in the map is visible from:
  
  all          every map cell;
  systematic   a systematic sample of map cells;
  random       a random sample of map cells;
  sites        a list of sites.

  3) Random sampling may occur with or without replacement.

  4) All types of sampling may be subject to binary masks such that:

     i)  Not all possible destination cells (cells to be `looked at')
         are of interest;

     ii) Not all possible viewpoints are of interest.

  5) It is possible to specify the height of the observer.

  6) It is possible to specify the height of the target object
     above ground.                                                     

  NOTE

  r.cva includes fixes for two bugs in the current version of r.los.
  These are that the observer elevation is cast to an integer, and that
  the patt_map does not effect the viewpoint and its immediately 
  neighbouring cells                                                   

module options:

input: elev raster
target_mask: binary raster map for target cells
viewpt_mask: binary raster map for viewpoints
sites: vector points layer
output: raster
obs_elev: ...
target_elev: ...
max_dist: ...
seed: for rand num gen
sample: sampling frequency as a percentage of cells  (default 10.0)
type: {all,systematic,random,sites} 
       all cells; systematic sample; random sample; sites map

-f  calculate the `visibility from' rather than `viewsheds of'
-m  differentiate masked cells from data value zero
-n  treat sample size as no. of cells
-r  allow replacement during random sampling

Hamish

[mmetz]

Replying to hamish:

note compiler warning: (x3)

> include/grass/iostream/replacementHeapBlock.h:146: warning: 'str' may be used uninitialized in this function

Hmm, not here, although I usually get warnings about uninitialized variables.

[...]

seems to keep -1 as invisibile not NULL ?????

fixed in r46201, came from original code

> [...]
> 
> # 180deg is just @ viewing cell, so for good colors we also need:  (PITA, so maybe we could build that into the module)
> [...]

Same as r.los, the viewing cell is defined to be straight below the viewpoint = 180 deg.

coord, -r desc: lat/lon -> easting/northing (seems to use current e,n LOCATION coords already..) drop the -r flag now that debugging stage is ending?

Sure, can be dropped.

tgt_elev=: how is that used if a target coord is never defined?

Assume you want to identify the cells where the top of a hypothetical structure (person, building, signal receiver, ...) is visible from the observation point, and it suffices if the top is visible, the base need not be visible.

timings on my trusty old P4:

[...]

indeed it seems to scale a bit better than r.los ... :-)

For small regions, you would mostly test the speed of reading the input map and writing the output map, but not the speed needed for LOS computation...

For GRASS 7, we are free to change module options, for GRASS 6, this patt_map option could be added for compatibility reasons, but the only effect would then be to print out a warning that a MASK should be used instead of patt_map=mask_map?

I wasn't thinking to make this as r.los2; as it is a totally diff't code base I'd add it named as r.viewshed, and drop r.los in grass7.

Makes sense...

but maybe it is prefered to keep the r.los name? if so then at least in GRASS 6.x it would be fine to have the option as "This does nothing and is kept for backwards compatibility reasons." and the warning as you describe. ? since r.los is so useless at modern region array sizes, maybe we should.. alternately, we could add r.viewshed in 6.x with a wrapper script for the old r.los name, like r.cats -> r.category.

In GRASS 7, other names changed too, e.g. v.db.addcol -> v.db.addcolumn

A quick glimpse into the manual with the list of raster modules should reveal the name of module that can do LOS analysis, a new name should be ok considering the vastly different code base.

Markus M

[hamish]

Replying to hamish:

180deg is just @ viewing cell, so for good colors we also need: (PITA, so maybe we could build that into the module)

Replying to mmetz:

Same as r.los, the viewing cell is defined to be straight below the viewpoint = 180 deg.

what my poor wording was trying to suggest: we make the module set color rules based on what's left of the overall range after 180 and -1 are taken out, much the same as r.slope.aspect or v.surf.rst will set color rules for their output maps. Not that we change the data value of 180 at the viewpoint, just to add a special color rule for that.

Hamish

[neteler]

Replying to hamish:

Replying to hamish:

180deg is just @ viewing cell, so for good colors we also need: (PITA, so maybe we could build that into the module)

Replying to mmetz:

Same as r.los, the viewing cell is defined to be straight below the viewpoint = 180 deg.

what my poor wording was trying to suggest: we make the module set color rules based on what's left of the overall range after 180 and -1 are taken out, much the same as r.slope.aspect or v.surf.rst will set color rules for their output maps. Not that we change the data value of 180 at the viewpoint, just to add a special color rule for that.

Any suggestions to get this one closed? "Just" a special color table is missing?

[mmetz]

Replying to neteler:

Any suggestions to get this one closed? "Just" a special color table is missing?

The precision bug is fixed, a nice color table is IMHO an enhancement, unrelated to this ticket. BTW, the module is still an addon for 6.x.

Closing ticket.

Markus M