Estimates of Local Sky Obstruction (Regions of Problematic GPS Reception)
Some possible approaches
- r.angle.sh (see image below and attached script)
r.mapcalc "$GIS_OPT_output = atan( ($GIS_OPT_elev - $z) / sqrt( (x() - $x)^2 + (y() - $y)^2) )"
- Compute percent of visible terrain that is above 10 degrees vertical elevation at each cell, using r.los.
# setup:
#
# 30 m resolution takes about 4.5 minutes on my laptop
# g.region res=30 -ap
# 10 m resolution takes about 19 minutes on my laptop
# g.region res=10 -ap
# dump x,y,z values from DEM
r.out.xyz elev10m > ~/temp/xyz
# setup coordinate file
echo -n "" > stats.xyp
g.remove --q rast=MASK
# loop through each cell
for cell in `cat ~/temp/xyz`
do
los_cmd=`echo $cell | awk -F"|" '{print "r.los --q --o in=elev10m out=los coordinate="$1","$2}'`
angle_cmd=`echo $cell | awk -F"|" '{print "r.mapcalc \"ang = if( los >= 100, 1, null()) \""}'`
eval $los_cmd
eval $angle_cmd
p=`r.stats -p -N ang | awk '{print $2}' | sed 's/\%//g'`
echo $cell | awk -v p=$p -F"|" '{ if(p != "") {print $1"|"$2"|"p}}' >> stats.xyp
done
# load results when done:
# r.in.xyz in=stats.xyp out=angp --o
Percent of visible terrain obscuring the sky above 10 degrees
