<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Compute naturalness areas
This notebook computes the naturalness areas for each point in a 'sample_points' vector map based on a 'naturalness' vector map. The computed areas are exported to a csv file.
This notebook uses GRASS GIS (7.0.4), and must be run inside of a GRASS environment (start the jupyter notebook server from the GRASS command line).
The code in this notebook was modified from examples in "How to write a Python GRASS GIS 7 addon" developed for a FOSS4G Europe 2015 workshop. The code is contained in the python-grass-addon repository on GitHub.com.
Required packages
numpy <br />
pandas <br />
pyprind
Variable definitions
naturalness – vector map representing naturalness <br />
naturalness_value_field – name of the field containing the naturalness values <br />
sample_points – vector map containing observer points <br />
radius – radius in map units around each sample point to calculate land cover areas within <br />
use_viewshed – consider only 'visible' land cover, based on a viewshed from each sample point <br />
Step1:  get viewshed suffix (for filenames)
Step2:  Import statements
Step3:  GRASS import statements
Step4:  Function declarations
connect to attribute table
Step5:  extract point from vector map
Step6:  create buffer around point
Step7:  convert viewshed from raster to vector map
Step8:  overlay a vector map on an underlying vector map using 'and' selection operator
Step9:  add area column to vector map
Step10:  create table showing total area by landcover type
Step11:  Utility functions
export vector map to a shapefile
Step12:  get info about a vector map
Step13:  Compute landcover areas
connect to 'sample_points' attribute table
Step14:  loop through sample points
Step15:  Export naturalness area table
create table with naturalness areas by type for each point
<ASSISTANT_TASK:>
Python Code:
naturalness = 'naturalness'
naturalness_value_field = 'value'
n_types = 7
sample_points = 'sample_points_field'
radius = 500
use_viewshed = False
Explanation: Compute naturalness areas
This notebook computes the naturalness areas for each point in a 'sample_points' vector map based on a 'naturalness' vector map. The computed areas are exported to a csv file.
This notebook uses GRASS GIS (7.0.4), and must be run inside of a GRASS environment (start the jupyter notebook server from the GRASS command line).
The code in this notebook was modified from examples in "How to write a Python GRASS GIS 7 addon" developed for a FOSS4G Europe 2015 workshop. The code is contained in the python-grass-addon repository on GitHub.com.
Required packages
numpy <br />
pandas <br />
pyprind
Variable definitions
naturalness – vector map representing naturalness <br />
naturalness_value_field – name of the field containing the naturalness values <br />
sample_points – vector map containing observer points <br />
radius – radius in map units around each sample point to calculate land cover areas within <br />
use_viewshed – consider only 'visible' land cover, based on a viewshed from each sample point <br />
End of explanation
def getViewshedSuffix():
    if use_viewshed:
        viewshed_suffix = '_viewshed'
    else:
        viewshed_suffix = ''
    return viewshed_suffix
naturalness_area_table_filename = "/home/ubuntu/naturalness_areas_{0}m{1}.csv".format(radius, getViewshedSuffix())
Explanation: get viewshed suffix (for filenames)
End of explanation
import numpy as np
import pandas
import pyprind
Explanation: Import statements
End of explanation
import grass.script as gscript
from grass.pygrass.vector.geometry import Point
from grass.pygrass.vector import Vector
from grass.pygrass.vector import VectorTopo
from grass.pygrass.vector.table import DBlinks
Explanation: GRASS import statements
End of explanation
def connectToAttributeTable(map):
    vector = VectorTopo(map)
    vector.open(mode='r')
    dblinks = DBlinks(vector.c_mapinfo)
    link = dblinks[0]
    return link.table()
Explanation: Function declarations
connect to attribute table
End of explanation
def extractPoint(input, ID, output):
    where = 'ID = {0}'.format(ID)
    type = 'point'
    gscript.read_command('v.extract',
                         input=input,
                         where=where,
                         output=output,
                         type=type, 
                         overwrite=True)
Explanation: extract point from vector map
End of explanation
def bufferPoint(input, output, radius):
    gscript.read_command('v.buffer',
                         input=input,
                         output=output,
                         type='point',
                         distance=radius,
                        overwrite=True)
Explanation: create buffer around point
End of explanation
def vectorizeViewshed(input, ID, output):
    type = 'area'
    column = 'visible'
    gscript.read_command('r.to.vect', 
                         input=input, 
                         output=output,
                         type=type,
                         column=column,
                         overwrite=True)
Explanation: convert viewshed from raster to vector map
End of explanation
def overlay(overlay, underlay, output):
    operator='and'
    gscript.read_command('v.overlay',
                         ainput=overlay,
                         binput=underlay,
                         operator=operator,
                         output=output,
                         overwrite=True)
Explanation: overlay a vector map on an underlying vector map using 'and' selection operator
End of explanation
def calculateAreas(map):
    
    #add new area column
    gscript.read_command('v.db.addcolumn',
                         map=map,
                         columns="area_square_meters DOUBLE PRECISION")
    
    #compute area and insert into area column
    gscript.read_command('v.to.db',
                         map=map,
                         type='centroid',
                         option='area',
                         columns='area_square_meters',
                         unit='meters')
Explanation: add area column to vector map
End of explanation
def getNaturalnessAreaByValue(map):
    #get area data
    table = connectToAttributeTable(map=map)
    table.filters.select()
    columns = table.columns.names()
    cursor = table.execute()
    result = np.array(cursor.fetchall())
    cursor.close()
    data = pandas.DataFrame(result, columns=columns).set_index('cat')
    
    #make sure naturalness_value_field is a numeric data type
    data['b_' + naturalness_value_field] = pandas.to_numeric(data['b_' + naturalness_value_field])
    
    #create naturalness column
    data['naturalness'] = np.nan
    
    #define naturalness value categories (round values)
    for index, row in data.iterrows():
        naturalness = np.round(row['b_' + naturalness_value_field])
        data.set_value(index, 'naturalness', int(naturalness))
    
    #sum areas by naturalness value
    data['area_square_meters'] = pandas.to_numeric(data['area_square_meters'])
    areas = data[['naturalness', 'area_square_meters']].groupby(by='naturalness').sum()
    
    #calculate mean
    total_area = data['area_square_meters'].sum()
    percent_area = data['area_square_meters'] / total_area
    weighted_area = percent_area * data['b_' + naturalness_value_field]
    mean = weighted_area.sum()
    
    #add to areas dataframe
    areas = areas.set_value('mean', 'area_square_meters', mean)
    
    return areas
Explanation: create table showing total area by landcover type
End of explanation
def exportVectorToShapefile(map, output):
    gscript.read_command('v.out.ogr',
                         input=map,
                         format='ESRI_Shapefile',
                         output=output,
                         flags='e',
                         overwrite=True)
Explanation: Utility functions
export vector map to a shapefile
End of explanation
def getVectorMapInfo(map):
    return gscript.read_command('v.info', map=map)
Explanation: get info about a vector map
End of explanation
point_table = connectToAttributeTable(map=sample_points)
point_table.filters.select()
columns = point_table.columns.names()
cursor = point_table.execute()
result = np.array(cursor.fetchall())
cursor.close()
point_data = pandas.DataFrame(result, columns=columns).set_index('cat')
Explanation: Compute landcover areas
connect to 'sample_points' attribute table
End of explanation
with Vector(sample_points, mode='r') as points:
    
    #setup progress bar
    progress_bar = pyprind.ProgBar(points.n_lines, bar_char='█', title='Naturalness analysis progress', monitor=True, stream=1, width=50)
    
    #iterate through points
    for point in points:
        
        #get point ID (SiteID)
        ID = point_data['ID'][point.cat-1]
        
        #update progress bar
        progress_bar.update(item_id=ID)
        
        #buffer current point
        extractPoint(input='sample_points_field', ID=ID, output='tmp_buffer_point')
        bufferPoint(input='tmp_buffer_point', output='tmp_point_buffer', radius=radius)
        
        #set buffer as overlay
        overlay_input = 'tmp_point_buffer'
        #consider only visible naturalness if 'use_viewshed' = True
        if use_viewshed:
            viewshed = 'vect_{0}_viewshed'.format(ID)
            visible_viewshed = 'vect_{0}_viewshed_{1}m'.format(ID, radius)
            #vectorize viewshed
            vectorizeViewshed(input='{0}_viewshed'.format(ID), ID=ID, output=viewshed)
            #overlay buffer on viewshed
            overlay(overlay=overlay_input,
                    underlay=viewshed,
                    output=visible_viewshed)
            #set overlay to the visible viewshed
            overlay_input = visible_viewshed
        overlay_output = 'vect_{0}_naturalness_{1}m{2}'.format(ID, radius, getViewshedSuffix())
        #overlay naturalness
        overlay(overlay=overlay_input,
                         underlay=naturalness,
                         output=overlay_output)
        
        #calculate naturalness area
        calculateAreas(map=overlay_output)
Explanation: loop through sample points
End of explanation
#create table
index_start = 0
''' set the first index number,
    allowing easier insertion into a database table that already contains 
    area calculations with other parameters
    (i.e radius and use_viewshed)'''
columns = ['ID', 'SiteID', 'IncludedArea']
columns = columns + [ str(n) for n in range(1,n_types+1) ] + ['mean']
area_table = pandas.DataFrame(columns=columns)
#set naming variables
included_area = '{0}m{1}'.format(radius, getViewshedSuffix())  #0=radius, 1=viewshed_suffix
map_pattern = 'vect_{0}_naturalness_{1}m{2}'  #0=ID, 1=radius, 2=viewshed_suffix
#iterate through points
for index, point in point_data.iterrows():
    ID = point['ID']
    map = map_pattern.format(ID, radius, getViewshedSuffix())
    #initiate row
    row = {'ID':"{0:.3g}".format(int(index) + index_start),
           'SiteID': str(ID),
           'IncludedArea': included_area}
    #get naturalness areas
    areas = getNaturalnessAreaByValue(map)
    
    #iterate through area types
    for index, area in areas.iterrows():
        #add area to row
        try:
            row["{0:.3g}".format(int(index))] = area['area_square_meters']
        except ValueError:
            row[index] = area['area_square_meters']
    #append row to table
    area_table = area_table.append(row, ignore_index=True)
area_table.set_index('ID', inplace=True)
#export table to file
area_table.to_csv(naturalness_area_table_filename, header=False)
Explanation: Export naturalness area table
create table with naturalness areas by type for each point
End of explanation
<END_TASK>