Talk:Workshop on urban growth modeling with FUTURES

Solutions to independent tasks

Distance from lakes/rivers

Distance to lakes or rivers may contribute to an amenity draw where people would like to build a house in a location where they can view water. First we will extract category Open Water from 2011 NLCD dataset:

r.mapcalc "water = if(landuse_2011 == 11, 1, null())"

Then we compute the distance to water with module r.grow.distance and set color table to shades of blue:

r.grow.distance input=water distance=dist_to_water
r.colors -n map=dist_to_water color=blues

Travel time to cities

First we specify the speed on different types of roads. We copy the roads raster into our mapset so that we can change it by adding a new attribute field speed. Then we assign speed values (km/h) based on the type of road:

g.copy vector=roads,myroads
v.db.addcolumn map=myroads columns="speed double precision"
v.db.update map=myroads column=speed value=50 where="MTFCC = 'S1400'"
v.db.update map=myroads column=speed value=100 where="MTFCC IN ('S1100', 'S1200')"

Now we rasterize the selected road types using the speed values from the attribute table as raster values.

v.to.rast input=myroads type=line where="MTFCC IN ('S1100', 'S1200', 'S1400')" output=roads_speed use=attr attribute_column=speed

We set the rest of the area to low speed and recompute the speed as time to travel through a 30m cell in minutes:

r.null map=roads_speed null=5
r.mapcalc "roads_travel_time = 1.8 / roads_speed"

Finally we compute the travel time to larger cities using r.cost:

r.cost input=roads_travel_time output=travel_time_cities start_points=cities
r.colors map=travel_time_cities color=byr

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="futures*")[current_mapset]
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)
    # 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

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