We will now put together what we have learned so far and write a script which computes viewsheds from points. Subsequently, we'll analyze the viewshed properties: viewshed size, number of visible points and the distance to the closest visible point. The dataset is the North Carolina sample database.
Both PyGRASS and GRASS GIS Python Scripting Library will be used in this session.
We will compute viewsheds on raster elevation (in PERMANENT mapset) from vector points.
To simplify the re-running of examples, we set this environmental variable, which allows direct overwriting of results from previous runs, bypassing the overwrite checks.
In [ ]:
import os
os.environ['GRASS_OVERWRITE'] = '1'
We first generate some random input points within a specified region using v.random (we use a fixed seed here for reproducibility):
In [ ]:
import grass.script as gscript
gscript.run_command('g.region', n=225200, s=222500, w=637500, e=640000, raster='elevation')
gscript.run_command('v.random', output='input_points', npoints=20, seed=2, quiet=True)
Before we start to loop over all vector points, we will first try out the workflow with a single point.
In [ ]:
from grass.pygrass.vector.geometry import Point
point = Point(638104, 223048)
We now compute the viewshed from this point using r.viewshed and then change all invisible cells to null values (no data) and all visible cells to value 1 using raster algebra r.mapcalc. This is needed for converting it to the vector data model further on (vectorization).
In [ ]:
elevation = 'elevation'
input_points = 'input_points'
viewshed_name = 'viewshed'
tmp_viewshed_name = 'tmp_viewshed'
viewshed_id = 1
gscript.run_command('r.viewshed', input=elevation, observer_elevation=3,
output=tmp_viewshed_name, coordinates=point.coords())
gscript.mapcalc(exp="{viewshed} = if({tmp}, {vid}, null())".format(viewshed=viewshed_name,
tmp=tmp_viewshed_name,
vid=viewshed_id))
As a first property of the viewshed, compute its area. The area is computed with r.univar as the number of non-null cells times the raster cell size. The result is in reported map units squared (i.e., square meters in this case).
In [ ]:
cells = gscript.parse_command('r.univar', map=viewshed_name,
flags='g')['n']
print(cells)
print gscript.region()
In [ ]:
area = float(cells) * gscript.region()['nsres'] * gscript.region()['nsres']
print(area)
The next task is to find and count the number of points contained in the input vector layer which are visible from the current point. One way to do this is to derive the vector layer of visible points by spatially overlapping the input points with the viewshed, for this see v.select. The viewshed must be first converted to the vector data model with r.to.vect.
In [ ]:
visible_points = 'tmp_points'
gscript.run_command('r.to.vect', input=viewshed_name, output=viewshed_name,
type='area')
gscript.run_command('v.select', ainput=input_points, atype='point',
binput=viewshed_name, btype='area',
operator='overlap', flags='t', output=visible_points)
We can now get the number of points using vector_info_topo, a wrapper function around v.info.
In [ ]:
n_points_visible = gscript.vector_info_topo(visible_points)['points']
print n_points_visible
The last viewshed property we want to compute is the distance from the current point to the closest visible point. Since v.distance requires vector input, we first save the current point as a separate vector layer using v.in.ascii.
In [ ]:
tmp_point = 'tmp_current_point'
if float(n_points_visible) >= 1:
gscript.write_command('v.in.ascii', input='-', stdin='%s|%s' % (point.x, point.y),
output=tmp_point)
distance = gscript.read_command('v.distance', from_=tmp_point, from_type='point', flags='p',
to=visible_points, to_type='point', upload='dist', dmin=1).strip()
distance = float(distance.splitlines()[1].split('|')[1])
else:
distance = 0
print(distance)
Now we will put the previous snippets together and compute the viewshed properties in a loop for all given input points.
Notice that we subtract 1 from the number of visible points since the current point is one of the visible ones.
In [ ]:
from grass.pygrass.vector import Vector
import grass.script as gscript
elevation = 'elevation'
input_points = 'input_points'
tmp_viewshed_name = 'tmp_viewshed'
tmp_visible_points = 'tmp_points'
tmp_point = 'tmp_current_point'
with Vector(input_points, mode='r') as points:
for point in points:
viewshed_id = str(point.cat)
viewshed_name = 'viewshed_' + viewshed_id
gscript.run_command('r.viewshed', input=elevation, observer_elevation=3,
output=tmp_viewshed_name, coordinates=point.coords())
gscript.mapcalc(exp="{viewshed} = if({tmp}, {vid}, null())".format(viewshed=viewshed_name,
tmp=tmp_viewshed_name,
vid=viewshed_id))
# viewshed size
cells = gscript.parse_command('r.univar', map=viewshed_name,
flags='g')['n']
area = float(cells) * gscript.region()['nsres'] * gscript.region()['nsres']
# visible points
gscript.run_command('r.to.vect', input=viewshed_name, output=viewshed_name,
type='area')
gscript.run_command('v.select', ainput=input_points, atype='point',
binput=viewshed_name, btype='area',
operator='overlap', flags='t', output=tmp_visible_points)
n_points_visible = gscript.vector_info_topo(tmp_visible_points)['points'] - 1
# distance to closest visible point
if float(n_points_visible) >= 1:
gscript.write_command('v.in.ascii', input='-', stdin='%s|%s' % (point.x, point.y),
output=tmp_point)
distance = gscript.read_command('v.distance', from_=tmp_point, from_type='point', flags='p',
to=tmp_visible_points, to_type='point', upload='dist', dmin=1).strip()
distance = float(distance.splitlines()[1].split('|')[1])
else:
distance = 0
print "%s, %d, %s, %.2f" % (viewshed_id, area, n_points_visible, distance)
Instead of printing the resulting properties on standard output, we save them into the attribute table of a new output vector layer. This we create by opening it in write mode and by passing as the parameters the columns of the attribute table.
In [ ]:
from grass.pygrass.vector import Vector
import grass.script as gscript
elevation = 'elevation'
input_points = 'input_points'
#
# output vector
#
output_points = 'output_points'
tmp_viewshed_name = 'tmp_viewshed'
tmp_visible_points = 'tmp_points'
tmp_point = 'tmp_current_point'
#
# define columns of the attribute table of the output vector
#
columns = [('cat', 'INTEGER'),
('area', 'DOUBLE PRECISION'),
('n_points_visible', 'INTEGER'),
('distance_to_closest', 'DOUBLE PRECISION')]
#
# we can open the input vector and create and open the output vector at once
#
with Vector(input_points, mode='r') as points, \
Vector(output_points, mode='w', tab_cols=columns) as output:
for point in points:
viewshed_id = str(point.cat)
viewshed_name = 'viewshed_' + viewshed_id
gscript.run_command('r.viewshed', input=elevation, observer_elevation=3,
output=tmp_viewshed_name, coordinates=point.coords())
gscript.mapcalc(exp="{viewshed} = if({tmp}, {vid}, null())".format(viewshed=viewshed_name,
tmp=tmp_viewshed_name,
vid=viewshed_id))
# viewshed size
cells = gscript.parse_command('r.univar', map=viewshed_name,
flags='g')['n']
area = float(cells) * gscript.region()['nsres'] * gscript.region()['nsres']
# visible points
gscript.run_command('r.to.vect', input=viewshed_name, output=viewshed_name,
type='area')
gscript.run_command('v.select', ainput=input_points, atype='point',
binput=viewshed_name, btype='area',
operator='overlap', flags='t', output=tmp_visible_points)
n_points_visible = gscript.vector_info_topo(tmp_visible_points)['points'] - 1
# distance to closest visible point
if float(n_points_visible) >= 1:
gscript.write_command('v.in.ascii', input='-', stdin='%s|%s' % (point.x, point.y),
output=tmp_point)
distance = gscript.read_command('v.distance', from_=tmp_point, from_type='point', flags='p',
to=tmp_visible_points, to_type='point', upload='dist', dmin=1).strip()
distance = float(distance.splitlines()[1].split('|')[1])
else:
distance = 0
#
# write each point with its attributes
#
output.write(point, (area, n_points_visible, distance))
output.table.conn.commit()
print "%s, %d, %s, %.2f" % (viewshed_id, area, n_points_visible, distance)
Finally we will make sure that the new vector layer was created and attributes properly written:
In [ ]:
with Vector(output_points, mode='r') as points:
# we can filter/sort the results
points.table.filters.select().order_by(u'area').get_sql()
print points.table.execute().fetchall()
In [ ]: