[GRASSLIST:4845] Re: 100KM raduis viewshed

[Glynn Clements]

Dan Jacobson wrote:

> Dear Sirs, what about curvature of the earth vs. line of sight?  The
> r.cva man page mentions that it and r.los don't take curvature of the
> earth into account.
> 
> You see, I'm trying to find out which mountains are visible from here
> in a 100KM radius.  Mainly the ones on the first photo on
> http://jidanni.org/location/view_en.html
> 
> I have a 40 meter DTM.  Shall I run it thru some mapcalc function of
> radius to get curvature of the earth accounted for?

That should work; 100km represents an angle of ~0.9 degrees, so
assuming that sin(x) == x results in a horizontal error of ~4.1
metres.

The following script seems about right.

#!/bin/sh
in=elevation.dem
out=corrected.dem
e=`g.region -p | fgrep east:  | awk '{print $2}'`
w=`g.region -p | fgrep west:  | awk '{print $2}'`
xc=`echo "($e+$w)/2.0" | bc`
n=`g.region -p | fgrep north: | awk '{print $2}'`
s=`g.region -p | fgrep south: | awk '{print $2}'`
yc=`echo "($n+$s)/2.0" | bc`
r.mapcalc "$out = eval(r = 6378137.0 , d = sqrt((x() - $xc)^2 + (y() - $yc)^2) , a = d/r , dy = r * (cos(a * 57.29577951308232) - 1.0) , $in + dy)"

> Tried r.los but I killed it after a half an hour of waiting.  I don't
> care about the angle of sight (90 or 91 degrees usually), by the way.
> I will try r.cva after I install the whole grass source just to get it
> compiled, someday.
> 
> I want to concentrate on just 30 degrees of my horizon.  That means I
> must use mapcalc to make a patt_map mask that any cell with the
> azimuths I want to be =1. A pie slice shaped mask. Hmmm. How to
> proceed?

Probably the simplest method would be to create a vector map
consisting of a single polygon, then convert it to a raster map with
v.to.rast.

Alternatively, you could compute the heading relative to the centre of
the map (or any other point) using the (undocumented) two-argument
form of r.mapcalc's atan() function, e.g.

	r.mapcalc "eval(heading = atan(x() - $xc, y() - $yc), ...

Note: the argument order is the reverse of ANSI C's atan2() function.

-- 
Glynn Clements <glynn.clements at virgin.net>