# Talk:Workshop on urban growth modeling with FUTURES

(diff) ← Older revision | Latest revision (diff) | Newer revision → (diff)

## 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.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
```

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
```