<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Script for computing and analyzing viewsheds
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
Step1:  We first generate some random input points within a specified region using v.random (we use a fixed seed here for reproducibility)
Step2:  Workflow for 1 point
Before we start to loop over all vector points, we will first try out the workflow with a single point.
Step3:  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).
Step4:  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).
Step5:  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.
Step6:  We can now get the number of points using vector_info_topo, a wrapper function around v.info.
Step7:  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.
Step8:  Workflow for multiple points
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.
Step9:  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.
Step10:  Finally we will make sure that the new vector layer was created and attributes properly written
<ASSISTANT_TASK:>
Python Code:
import os
os.environ['GRASS_OVERWRITE'] = '1'
Explanation: Script for computing and analyzing viewsheds
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.
End of explanation
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)
Explanation: We first generate some random input points within a specified region using v.random (we use a fixed seed here for reproducibility):
End of explanation
from grass.pygrass.vector.geometry import Point
point = Point(638104, 223048)
Explanation: Workflow for 1 point
Before we start to loop over all vector points, we will first try out the workflow with a single point.
End of explanation
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))
Explanation: 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).
End of explanation
cells = gscript.parse_command('r.univar', map=viewshed_name,
                              flags='g')['n']
print(cells)
print gscript.region()
area = float(cells) * gscript.region()['nsres'] * gscript.region()['nsres']
print(area)
Explanation: 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).
End of explanation
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)
Explanation: 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.
End of explanation
n_points_visible = gscript.vector_info_topo(visible_points)['points']
print n_points_visible
Explanation: We can now get the number of points using vector_info_topo, a wrapper function around v.info.
End of explanation
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)
Explanation: 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.
End of explanation
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)
    
Explanation: Workflow for multiple points
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.
End of explanation
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)
        
    
Explanation: 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.
End of explanation
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()
Explanation: Finally we will make sure that the new vector layer was created and attributes properly written:
End of explanation
<END_TASK>