# 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

# 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

r.angle_.sh_.txt