November workshop

FUTURES workshop

Date: Monday, November 9 at 3:00 PM in 5119

Data will be provided in form of GRASS GIS location. Our case study is Triangle. Participants should bring their laptops with GRASS GIS 7 and R (≥ 3.0.2) installed.

Outline

Brief overview of r.futures components
Demand - r.futures.demand
Potential - sampling using r.sample.category and running R for selecting model
Calibration - r.futures.calib - only limited range of values to calibrate
PGA - r.futures.pga - for one county

Overview

See r.futures manual page to review all the components.

Input data for case study:

Spatial:

developed/undeveloped rasters for 1990 2000 and 2010 derived from NLCD
subregions - raster with counties
predictors

distance to protected areas
slope
road density
distance to rivers/lakes
forests
agriculture suitability (soils)
distance to urban centers (traveltime)

Non-spatial:

past population for 1990 2000 and 2010 for each county and projections

For this workshop, we will have data already prepared in the right formats.

Data preparation

First we will set computational region of our analyses:

g.region n=1595760 s=1494780 w=1509630 e=1608600 res=30 -pa

We will extract 2 selected counties (Wake and Johnston) and convert them to raster with the values equal the FIPS attribute which links to population file:

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

Then we change the categories of the raster so that they start with 1 (needed for FUTURES simulation):

r.clump input=selected_counties_id output=selected_counties_1

Demand

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

First we set the computational region to match our study area:

g.region raster=subregions

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 the 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_1990,urban_2000,urban_2010 subregions=selected_counties_id observed_population=demand_input_known_pop.csv projected_population=demand_input_predicted_pop.csv simulation_times=2010,2011,2012,2013,2014,2015,2016,2017,2018,2019,2020,2021,2022,2023,2024,2025,2026,2027,2028,2029,2030 method=logarithmic plot=plot demand=demand.csv separator=comma

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

Potential

Preprocessing of the predictor variables

Slope

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

r.slope.aspect elevation=DEM slope=slope

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=rdens_1km 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=euc_interchanges

Development pressure

Computing development pressure with r.futures.devpressure:

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

When gamma increases, development influence decreases more rapidly with distance. Size is 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_2010 raster over it.

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

Sampling

First we will rescale some of our input variables:

r.mapcalc expression="rdens_1km_100 = rdens_1km  * 100"
r.mapcalc expression="forest_1km_100 = forest_1km  * 100"
r.mapcalc expression="farm_1km_100 = farm_1km * 100"
r.mapcalc expression="euc_interchanges_100 = euc_interchanges / 100"

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

r.mapcalc "urban_2010_clip = if (selected_counties_1, urban_2010)"

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

r.sample.category input=urban_2010_clip output=sampling sampled=selected_counties_1,slope,devpressure_0_5,rdens_1km_100,water_dist,forest_1km_100,farm_1km_100,euc_protect,euc_interchanges_100 npoints=5000,1000

And export the attribute table as CSV file:

v.db.select map=sampling columns=urban_2010_clip,selected_counties_1,slope,devpressure_0_5,rdens_1km_100,water_dist,forest_1km_100,farm_1km_100,euc_protect,euc_interchanges_100 separator=comma file=samples.csv

By selecting only certain columns for export, we can specify the input predictors.

Then run R script for automatic model selection which outputs summary of the "best" model and exports a CSV file with the model coefficients:

Rscript potential.R -i samples.csv -l selected_counties_1  -r urban_2010_clip -m 3 -x 7 -o potential.csv

Note: We must edit the potential file and remove double quotes from the names of rows, otherwise PGA won't read it correctly!

Creating incentive table

We need to prepare incentive table, which changes probabilities on the landscape with a power function. It's easy to generate one using python commands in r.futures.pga module:

import numpy as np
power = 2
length = 1001
x = np.linspace(0, 1, length)
np.savetxt(fname="incentive_table.txt", fmt='%1.3f', delimiter=',', X=np.column_stack((x, np.power(x, power))), header=",{}".format(length), comments='')

Paste those commands for example to the Python Shell in GRASS (line by line). You can change the power to a number between 0.25 and 4 to test urban sprawl/infill scenarios.

Calibration

We will calibrate patches on Wake county only. We first extract it:

r.mapcalc "wake = if(selected_counties_1 == 1, 1, null())"
g.region zoom=wake -p

Module r.futures.calib calls r.futures.pga. So we have to have all the inputs prepared.

The order of predictors must match the order of columns in potential.csv.

r.futures.calib development_start=urban_2000 development_end=urban_2010 repeat=1 compactness_mean=0.4 compactness_range=0.1 discount_factor=0.6 patch_threshold=1800 patch_sizes=patches.txt calibration_results=calib.csv nprocs=1 development_pressure=devpressure_0_5 predictors=euc_interchanges_100,euc_protect,forest_1km_100,rdens_1km_100,water_dist n_dev_neighbourhood=10 devpot_params=potential.csv incentive_table=incentive_table.txt num_neighbors=4 seed_search=2 development_pressure_approach=gravity gamma=0.5 scaling_factor=1 num_regions=2 subregions=wake demand=demand.csv

PGA