FUTURES tutorial

r.futures.* is an implementation of FUTure Urban-Regional Environment Simulation (FUTURES) which is a model for multilevel simulations of emerging urban-rural landscape structure. FUTURES produces regional projections of landscape patterns using coupled submodels that integrate nonstationary drivers of land change: per capita demand (DEMAND submodel), site suitability (POTENTIAL submodel), and the spatial structure of conversion events (PGA submodel).

This tutorial shows how to prepare data and run the model. You can download sample dataset from here (TODO) and simulate urban growth in The Triangle, region in the Piedmont of North Carolina in the United States with rapidly growing cities Raleigh, Durham and Chapel Hill.

Software

Required software includes:

• GRASS GIS 7
• r.futures addon
• R (≥ 3.0.2)
• SciPy (optional)

Input data used in this tutorial

• digital elevation model (NED)
• NLCD 2001, 2011
• NLCD 1992/2001 Retrofit Land Cover Change Product
• transportation network (TIGER)
• county boundaries (TIGER)
• protected areas (Secured Lands)
• cities above 20000 (USGS)
• county population past estimates and future projections (NC OSBM) per county (nonspatial)

Workflow

Initial data preparation

First we will set computational region of our analyses to an extent covering our study area and so that the cells are aligned with one of our landuse rasters:

g.region n=1648005 s=1480695 e=1624815 w=1462965 align=landuse_2001 -p

We will derive urbanized areas from NLCD dataset for year 2001 and 2011 by extracting categories category 21 - 24 into a new binary map where developed is 1, undeveloped 0 and NULL (no data) is area unsuitable for development (water, wetlands, protected areas). First we will convert protected areas from vector to raster and set NULLs to zeros (for simpler raster algebra expression in the next step):

v.to.rast input=protected output=protected use=val
r.null map=protected null=0

And then create rasters of developed/undeveloped areas using raster algebra:

r.mapcalc "urban_2001 = if(landuse_2001 >= 21 && landuse_2001 <= 24, 1, if(landuse_2001 == 11 || landuse_2001 >= 90 || protected, null(), 0))"
r.mapcalc "urban_2011 = if(landuse_2011 >= 21 && landuse_2011 <= 24, 1, if(landuse_2011 == 11 || landuse_2011 >= 90 || protected, null(), 0))"

We will convert vector counties to raster with the values of the FIPS attribute which links to population file:

v.to.rast input=subregions_triangle type=area use=attr attribute_column=FIPS output=selected_counties_id

If we need to run FUTURES only for smaller number of counties, we can extract them based on the identifier:

v.to.rast input=subregions_triangle type=area where="FIPS in (37183, 37101)" use=attr attribute_column=FIPS output=selected_counties_id

Potential submodel

Predictors

Slope

We will derive slope in degrees from DEM using r.slope.aspect module:

r.slope.aspect elevation=DEM slope=slope

Distance from lakes/rivers

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

Distance from protected areas

We will use raster protected of protected areas we already created. We compute the distance to protected areas with module r.grow.distance and set color table from green to red:

r.grow.distance input=protected distance=dist_to_protected
r.colors map=dist_to_protected color=gyr

Forests

We will smooth the transition between forest and other land use, see NLCD category Open Water:

r.mapcalc "forest = if(landuse_2011 >= 41 && landuse_2011 <= 43, 1, 0)"
r.neighbors -c input=forest output=forest_smooth size=15 method=average
r.colors map=forest_smooth color=ndvi

Travel time to cities

Here we will compute travel time to cities (> population 20000) as cumulative cost distance where cost is defined as travel time on roads. 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 based on the type of road:

g.copy vector=roads_2015,roads
v.db.addcolumn map=roads columns="speed double precision"
v.db.update map=roads column=speed value=50 where="MTFCC = 'S1400'"
v.db.update map=roads 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=roads 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_20000
r.colors map=travel_time_cities color=byr

Road density

We will rasterize roads and use moving window analysis (r.neighbors) to compute road density:

v.to.rast input=roads_2015 output=roads use=val
r.null map=roads null=0
r.neighbors -c input=roads output=road_dens size=25 method=average

Distance to interchanges

We will consider TIGER roads of type Ramp as interchanges, rasterize them and compute euclidean distance to them:

v.to.rast input=roads_2015 where="MTFCC = 'S1630'" output=interchanges use=val
r.grow.distance -m input=interchanges distance=dist_interchanges

Development pressure

Computing development pressure with r.futures.devpressure:

r.futures.devpressure -n input=urban_2011 output=devpressure_0_5 method=gravity size=20 gamma=0.5 scaling_factor=0.001

When gamma increases, development influence decreases more rapidly with distance. Size is half the size of the moving window. To explore the effects of these parameters, you can view the development pressure as a surface with draped urban_2011 raster over it.

TODO We should derive more layers with different gamma parameters for the potential statistical model.

Rescaling variables

First we will look at the ranges of our predictor variables by running a short Python code snippet in Python shell:

for name in ('slope', 'dist_to_water', 'dist_to_protected', 'forest_smooth', 'travel_time_cities', 'road_dens', 'dist_interchanges', 'devpressure_0_5'):
    minmax = grass.raster_info(name)
    print name, minmax['min'], minmax['max']

We will rescale some of our input variables:

r.mapcalc "dist_to_water_km = dist_to_water / 1000"
r.mapcalc "dist_to_protected_km = dist_to_protected / 1000"
r.mapcalc "dist_interchanges_km = dist_interchanges / 1000"
r.mapcalc "road_dens_perc = road_dens * 100"
r.mapcalc "forest_smooth_perc = forest_smooth * 100"

Sampling

To sample only in the analyzed counties, we will clip development layer:

r.mapcalc "urban_2011_clip = if(selected_counties_id, urban_2011)"

To estimate the number of sampling points, we can use r.report to report number of developed/undeveloped cells and their ratio.

r.report map=urban_2011 units=h,c,p

We will sample the predictors and the response variable with 10000 random points in undeveloped areas and 2000 points in developed area:

r.sample.category input=urban_2011_clip output=sampling sampled=selected_counties_id,devpressure_0_5,slope,road_dens_perc,forest_smooth_perc,dist_to_water_km,dist_to_protected_km,dist_interchanges_km,travel_time_cities npoints=10000,2000

The attribute table can be exported as CSV file (not necessary step):

v.db.select map=sampling columns=urban_2011_clip,selected_counties_id,devpressure_0_5,slope,road_dens_perc,forest_smooth_perc,dist_to_water_km,dist_to_protected_km,dist_interchanges_km,travel_time_cities separator=comma file=samples.csv

Development potential

Now we find best model for predicting urbanization using r.futures.potential which wraps an R script.

r.futures.potential input=sampling output=potential.csv columns=devpressure_0_5,slope,road_dens_perc,forest_smooth_perc,dist_to_water_km,dist_to_protected_km,dist_interchanges_km,travel_time_cities developed_column=urban_2011_clip subregions_column=selected_counties_id min_variables=4

Demand

TODO: check if we need to mask roads

v.to.rast input=roads type=line output=roads_mask use=val
r.mask roads_mask -i

We will use r.futures.demand which derives the population vs. development relation. The relation can be linear/logarithmic/exponential/exponential approach.

The format of the input population CSV files is described in the manual. It is important to have synchronized categories of subregions and the column headers of the CSV files (in our case FIPS number). How to simply generate the list of years (for which demand is computed) is described in r.futures.demand manual.

To simplify things, change your current working directory to where your input CSV files are - you can do that from menu Settings → GRASS working environment → Change working directory.

r.futures.demand development=urban_1992,urban_2001,urban_2011 subregions=selected_counties_id observed_population=population_trend.csv projected_population=population_projection.csv simulation_times=2011,2012,2013,2014,2015,2016,2017,2018,2019,2020,2021,2022,2023,2024,2025,2026,2027,2028,2029,2030,2031,2032,2033,2034,2035 plot=plot_demand.pdf demand=demand.csv separator=comma

In your current working directory, you should find files plot_demand.png and demand.csv.