[GRASS-user] query values within different angles

Dear GRASS users,

I need to query raster values around a city for different angles. The result I am aiming at should be the amount of landcover surrounding the city within 10 degrees angles and different distances to the city.

So far I tried r.buffer and v.buffer but that are not providing what I am lookin for. Using r.mapcalc and addressing each pixel separately would be an option but I hope that GRASS combined with perhaps PGSQL and Postgis functions provides an easier way to achieve that.

I would highly appreciate any hints how to derive buffers for specific angles only.

kind regards, Tim

P.S.: using GRASS 6.4.1. under Windows, QGIS 1.6.0 and the newest R version

Tim wrote:

I need to query raster values around a city for different
angles. The result I am aiming at should be the amount of
landcover surrounding the city within 10 degrees angles and
different distances to the city.

So far I tried r.buffer and v.buffer but that are not
providing what I am lookin for. Using r.mapcalc and
addressing each pixel separately would be an option but I
hope that GRASS combined with perhaps PGSQL and Postgis
functions provides an easier way to achieve that.

I would highly appreciate any hints how to derive buffers
for specific angles only.

Hi,

in years past I had a version of r.los for GRASS 5.0 where instead
of reporting the vertical angle to the source point I had it
report horizontal compass heading. I'm sure it's on an old hard
drive somewhere, but in hindsight using r.mapcalc and some simple
trig would be a heck of a lot easier and more efficient than
anything else. No other tools are needed beyond that one amazing
module.

#spearfish
g.region -d
g.region -c
  center easting: 599490.000000
  center northing: 4920855.000000

Ce=599490
Cn=4920855

r.mapcalc "angle_to_pt.cartesian = atan(x() - $Ce, y() - $Cn)"
r.mapcalc "angle_to_pt.compass.1 = 90 - angle_to_pt.cartesian"
r.mapcalc "angle_to_pt.compass = if(angle_to_pt.compass.1 < 0, \
   angle_to_pt.compass.1 + 360, angle_to_pt.compass.1)"
g.remove angle_to_pt.compass.1

MIN_ANGLE=0
MAX_ANGLE=10

r.mapcalc "MASK = if(angle_to_pt.compass >= $MIN_ANGLE \
   && angle_to_pt.compass < $MAX_ANGLE, 1, null() )"

d.redraw

Hamish

Thank you Hamish,

worked pretty good for my purposes. But as I am a beginner in GRASS I need to ask what the $ sign stands for. Because when I applied your solution with the $-sign and just copied it the resulting maps were rather coarse in resolution. So I tried and left the $ sign away. The result was a map in the same resolution as the original map.

Thanks for quick help!

Cheers Tim

Am 31.05.2011 12:33, schrieb Hamish:

Tim wrote:

I need to query raster values around a city for different
angles. The result I am aiming at should be the amount of
landcover surrounding the city within 10 degrees angles and
different distances to the city.

So far I tried r.buffer and v.buffer but that are not
providing what I am lookin for. Using r.mapcalc and
addressing each pixel separately would be an option but I
hope that GRASS combined with perhaps PGSQL and Postgis
functions provides an easier way to achieve that.

I would highly appreciate any hints how to derive buffers
for specific angles only.

Hi,

in years past I had a version of r.los for GRASS 5.0 where instead
of reporting the vertical angle to the source point I had it
report horizontal compass heading. I'm sure it's on an old hard
drive somewhere, but in hindsight using r.mapcalc and some simple
trig would be a heck of a lot easier and more efficient than
anything else. No other tools are needed beyond that one amazing
module.

#spearfish
g.region -d
g.region -c
   center easting: 599490.000000
   center northing: 4920855.000000

Ce=599490
Cn=4920855

r.mapcalc "angle_to_pt.cartesian = atan(x() - $Ce, y() - $Cn)"
r.mapcalc "angle_to_pt.compass.1 = 90 - angle_to_pt.cartesian"
r.mapcalc "angle_to_pt.compass = if(angle_to_pt.compass.1< 0, \
    angle_to_pt.compass.1 + 360, angle_to_pt.compass.1)"
g.remove angle_to_pt.compass.1

MIN_ANGLE=0
MAX_ANGLE=10

r.mapcalc "MASK = if(angle_to_pt.compass>= $MIN_ANGLE \
    && angle_to_pt.compass< $MAX_ANGLE, 1, null() )"

d.redraw

Hamish