"G" == Glynn Clements <glynn.clements@virgin.net> writes:
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?
G> Probably the simplest method would be to create a vector map
G> consisting of a single polygon, then convert it to a raster map with
G> v.to.rast.
I wouldn't know where to begin there. Alternatively?
G> Alternatively, you could compute the heading relative to the centre
G> of the map (or any other point) using the (undocumented)
G> two-argument form of r.mapcalc's atan() function, e.g.
G> r.mapcalc "eval(heading = atan(x() - $xc, y() - $yc), ...
G> Note: the argument order is the reverse of ANSI C's atan2() function.
$ xc=235577 yc=2675353 ;\
r.mapcalc "eval(heading = atan(x() - $xc, y() - $yc), \
nantou.triang = heading>140 && heading<170)"
you have confused me
I even tried if(heading...)
--
http://jidanni.org/ Taiwan(04)25854780
On 11 Nov 2002, Dan Jacobson wrote:
>>>>> "G" == Glynn Clements <glynn.clements@virgin.net> writes:
G> Alternatively, you could compute the heading relative to the centre
G> of the map (or any other point) using the (undocumented)
G> two-argument form of r.mapcalc's atan() function, e.g.
G> r.mapcalc "eval(heading = atan(x() - $xc, y() - $yc), ...
G> Note: the argument order is the reverse of ANSI C's atan2() function.
$ xc=235577 yc=2675353 ;\
r.mapcalc "eval(heading = atan(x() - $xc, y() - $yc), \
nantou.triang = heading>140 && heading<170)"
you have confused me
I even tried if(heading...)
What about
xc=235577 yc=2675353 ;\
r.mapcalc "nantou.triang = eval(heading = atan(x() - $xc, y() - $yc), \
heading>140 && heading<170)"
So the last expression in the eval is the one that is assigned to every
cell in the output map. This worked for me.
--
http://jidanni.org/ Taiwan(04)25854780
Dan Jacobson wrote:
G> Alternatively, you could compute the heading relative to the centre
G> of the map (or any other point) using the (undocumented)
G> two-argument form of r.mapcalc's atan() function, e.g.
G> r.mapcalc "eval(heading = atan(x() - $xc, y() - $yc), ...
G> Note: the argument order is the reverse of ANSI C's atan2() function.
$ xc=235577 yc=2675353 ;\
r.mapcalc "eval(heading = atan(x() - $xc, y() - $yc), \
nantou.triang = heading>140 && heading<170)"
you have confused me
Sorry, it should have been:
r.mapcalc "$out = eval(heading = atan(x() - $xc, y() - $yc), ...
Paul's suggestion should work.
The input to r.mapcalc has to be of the form "name = expression", or a
sequence of such e.g. "name1 = expression1 ; name2 = expression2".
Assignments at the top level are treated differently to assignments
within an expression. An assignment at the top level creates a new
map, while assignments within an expression don't.
--
Glynn Clements <glynn.clements@virgin.net>
Actually by trial and error I found I had to use
r.mapcalc "pie=eval(heading=360+90-atan(x()-$xc,y()-$yc),heading>140&&heading<170)"
to make the 30 degree pie shaped mask in the SE quadrant like I wanted.
But even using a 2 degree wide sliver still didn't stop r.los from
taking hours and hours. So I had to put r.los out of its misery.
Apparently if you use small square masks, r.los might even finish sooner or
later, But the man page gives no strategies, and the program gives no
progress indications.
It seems that nice parameters for r.los might be starting_azimuth and
ending_azimuth, but I was told r.los is set in stone and will not be
enhanced.
And r.los doesn't have a switch to turn off the computing the angles
of sight for each cell if one does not need them.
And r.cva is not in GRASS officially etc. and it and grass 5.0.0
itself are not ready in a .deb ...
Some mail:
I'm happy for you to include r.cva if you think it's appropriate.
But how urgent is this? A couple of months back I converted r.los to
make proper use of INT, FCELL or DCELL and so should be able to do
something similar for r.cva over Christmas (incidentally, you might
also like the revised r.los code - let me know?).
Markus> Yes, the revised r.los code is interesting for us/me 
Markus> Whenever the r.cva is updated, please let me know.
--
http://jidanni.org/ Taiwan(04)25854780
In my recent use of r.mapcalc for 1. bending of DTMs to simulate
curvature of the earth, 2. making layers based on distance for
coloring by distance instead of elevation, I used trig functions and sqrt
producing floating number results. Well, I will wrap the formulae
next time in round(), to get integer results, as 1. the DTMs were
integers to begin with, 2. I recall there are several advantages with
integer files.
So for folks reading this thread: remember to use round().
--
http://jidanni.org/ Taiwan(04)25854780
We have considered the earth's curvature, but um, "The Effect Of
Atmospheric Refraction On The Observed Elevation Angles Of Peaks"
http://tchester.org/sgm/analysis/peaks/refraction.html indeed affected
my results. Without considering it some peaks visible in binoculars
would sink behind other peaks in d.3d's drawing. So now I'm OK with
$ xc=235480 yc=2675355
$ r.mapcalc "nantou.bend=round(nantou-((x()-$xc)^2+(y()-$yc)^2)*.0000000675)"
which I found in a local textbook.
I suppose I won't fret about height and temperature, as they seem to
cancel out on the above webpage, though I could easily utilize FROM
and TO heights if there was a simple formula available.
--
http://jidanni.org/ Taiwan(04)25854780