Talk:FUTURES land-change modeling for evaluating innovative conservation scenarios

Python: Loss of farmland and forest

import numpy as np
import grass.script as gscript 
 
# list maps computed by r.futures.parallelpga for one scenario
current_mapset = gscript.gisenv()['MAPSET']
# put the name you are using for the scenarios in the pattern option 
maps = gscript.list_grouped('raster', pattern="final*")[current_mapset]
developed = []
forest = []
agriculture = []
other = []
for each in maps:
    # compute landcover in 2035 by updating landcover 2011 with the simulation results:
    gscript.mapcalc("landuse_2035_{name} = if({name} >=1, 21, landuse_2011)".format(name=each), quiet=True, overwrite=True)
    # get information about how many cells there are of each land use category, use r.stats
    res = gscript.read_command('r.stats', flags='cn', input="landuse_2035_{name}".format(name=each), quiet=True).strip()
    for line in res.strip().splitlines():
        cat, cells = line.split()
        cat = int(cat)
        cells = int(cells)
        # developed
        if 21 <= cat <= 24:
            developed.append(cells)
        # agriculture
        elif cat == 81 or cat == 82:
            agriculture.append(cells)
        # forest
        elif 41 <= cat <= 43:
            forest.append(cells)
        # other land use
        else:
            other.append(cells)
 
# convert to NumPy arrays
forest = np.array(forest)
agriculture = np.array(agriculture)
other = np.array(other)
developed = np.array(developed)
 
print "forest: {mean} +- {std}".format(mean=np.mean(forest), std=np.std(forest))
print "agriculture: {mean} +- {std}".format(mean=np.mean(agriculture), std=np.std(agriculture))
print "other: {mean} +- {std}".format(mean=np.mean(other), std=np.std(other))
print "developed: {mean} +- {std}".format(mean=np.mean(developed), std=np.std(developed))

Python: Forest fragmentation

We will need the addon r.forestfrag. You can install it in the command console with:

g.extension r.forestfrag

import grass.script as gscript 

def forest_fragmentation(landuse):
    gscript.mapcalc("forest_2035 = if({landuse} >= 41 && {landuse} <= 43, 1, 0)".format(landuse=landuse))
    gscript.run_command('r.forestfrag', flags='a', input='forest_2035', output='fragment', window=25)
    result = gscript.parse_command('r.stats', input='fragment', flags='cn', parse=(gscript.parse_key_val, {'sep': ' ', 'val_type': int}))
    gscript.run_command('g.remove', type='raster', name=['forest_2035', 'fragment'], flags='f')
    return result

current_mapset = gscript.gisenv()['MAPSET']
landuse_maps = gscript.list_grouped('raster', pattern="landuse_2035_*")[current_mapset]
landuse_maps = landuse_maps[:5]
for landuse in landuse_maps:
    result = forest_fragmentation(landuse)
    print result