Creating a Minimal LULCC Scenario using a Bayesian Model
Ryan Thomas
May 9, 2018
Background
A core objective of urban and regional planning is to shape land uses by implementing zoning codes, land use restrictions, and incentives that guide land development. Sustainability has become a prominent subtext to land use policies as the contributions of land use change to climate change have been revealed through research. Complex relationships between human and non-human systems determine land use change, and understanding the integration of these systems remains a key challenge for planners attempting to shape land use patterns. The institutional impacts on patterns and rates of urbanization are variable in different locations, and it is difficult to identify best practices of land use regulation. Despite the lack of strong evidence in the efficacy of land use planning at large scales land use continues to be a focal point of planning practice. A wide array of stakeholders, or publics, can influence future land uses through their actions and decisions (Burby 2003). This is, in part, why Burby (2003) argues that including diverse groups in the planning process leads to greater plan implementation by local governments. However, even stakeholders who are not included in planning processes have an impact on how land is used, and land use regulations typically place intensity ceilings, rather than floors. For example, some land zoned for agriculuture also allows low-density residential uses and may in practice be dominated these residential use. As such, understanding stakeholder preferences can be a valuable addition to land use modeling and prediction of likely future uses. This is a justification for including stakeholder preferences in land use models, a common feature of Bayesian models. These land use changes will in turn have fiscal impacts on municipal boundaries.
One promising body of literature in this area is developing spatial models that include institutions, existing land use patterns and rates of urbanization, as well as geocoded interviews to determine where land use intensity is likely to increase. Urbanization models have long utilized remote sensing to monitor land use change because of its consistent measurement of Earth’s surface over time. In the last two decades, satellite data sources have been increasingly paired with institutional analysis and interviews at multiple levels political and social structures, including government and households.
Yet, the science of these impacts is underdeveloped, as is the craft of land use policy implementation that achieves desired results. It is unclear what land use patterns planners should plan for and, assuming that question is answered, the mechanisms that would enable achieving it. Existing research on the drivers of land use change paints a dim picture for the impact of planning on large-scale patterns of urban development. Studies looking at large scale urbanization rates show little or no impact of local policy on regional rates of growth. However, localized studies, especially those that balance institutional analysis with observed urbanization patterns, show that local politics, institutional dynamics, and policies do affect urban development patterns. The combination of these two bodies of research suggest that effective local land use policies may push development capital into more favorable policy environments. Yet, the field is in need of new theoretical and methodological frameworks or mechanisms for monitoring leakage and spillover effects of land use policy. This paper explores one method of measuring land use change over time along the Massachusetts coast in a reproducible way. By expanding the analytical methods described here (e.g. by including more covariates or by collecting survey data on land use preferences), urban planners could establish baselines, targets, and develop scenarios to examine the social and institutional determinants of land use change.
Brief literature review
Drivers of urbanization include macro factors and local factors. Macro factors are often aggregated to the national or regional/state levels, including national GDP levels, global exports, fiscal transfers to regional governments, national transportation costs, population change. These factors are usually present in large-scale global and regional studies of land use change. Except for transportation costs, which has a negative relationship with urbanization rates, increases in the values of all of these factors have been shown to increase urbanization rates (Seto et al. 2011, Zhang et al. (2014)). This connection poses a particular threat to global ecology, since increasing the values of these macro factors (again, except for transportation costs) has also been conceptually linked to moral mandate to improve quality of life and reduce poverty. This linkage between macro economic factors and urban expansion rates suggests that local land use planning, which is often used as a mechanism for economic development, may have an impact on urbanization rates. At the national level, fiscal transfers have also shown to increase rates of urban expansion. Transportation costs, generally proxied by fuel costs and assumed to be constant at the national level, have shown negative impact urban expansion rates – lower transportation costs correspond to higher urban expansion rates (Torrens 2000). The final macro factor that commonly found in studies of urbanization is population change, which positively impacts urbanization rates due to housing demand. Population change, though, can typically be measured at a local level in smaller scale studies and does not present a major challenge for estimating local models that could reveal differences in policy effectiveness.
Micro, or local, factors are typically studied through situated case studies that pay closer attention to power dynamics and institutional arrangements difficult to extrapolate to larger scales. Transportation infrastructure at this scale is typically measured by proximity to arterial roads and public transportation. Population dynamics at the local level include migration trends and the behavior of in-migrants and out-migrants that might develop particular areas of the city differentially. The most important difference between large-scale studies and smaller scale studies’ attention to drivers of urban land expansion is the presence of inter-jurisdictional relationships and land use policy in the latter. Interjursditcional dynamics are necessary to avoid local effects of leakage or spillover where restrictive policies in one municipality put added pressure on neighboring municipalities with the same levels of macro drivers. Coordination is needed for policy coherence at larger scales, and my previous research shows that multi-scalar governance can contribute to sound regional development patterns. Finally, local level models often include suitability analysis such as geological characteristics and critical ecological areas.
Research examining large areas has shown that, at the regional level, local land use plans have little influence on the rates of urbanization (Schmidt, Siedentop, and Fina 2017). However, local-level analyses show that urban politics can undermine strong land use policies, while multi-city and multi-level networks may improve coherence in land use policy implementation (Pierce et al. 2005). It is possible that different levels of analysis would show different results with respect to the rates of urban development in an ecological fallacy.
Numerous models of land use change have been piloted at regional scales with a range of the above drivers included as determinants. Preliminary research revealed that a BBN approach allows the incorporation of survey data indicating expert opinion, perceptions of land use decision-makers (Celio, Koellner, and Grêt-Regamey 2014), or firm-level decisions (Hashiguchi and Tanaka 2015) into the model. While survey data will not be included in the paper, other estimation approaches will be used to estimate a likelihood of conversion similar to ‘expert opinion’. Other land use change models could be used to estimate how likely are ares to be converted. Two common models for esimating and predicting land use conversion are CLUE-S and SLEUTH. The SLEUTH model is a cellular automata model that is based on the slope, land use, exclusion, urban extent, transportation, and hillshade for each cell (Silva and Clarke 2002, K. C. Clarke (2008), Hua et al. (2014), Park et al. (2017)). The initial plan for this study was to compare multiple models of land use change, including the BBN method with the SLEUTH model with the CLUE-S model. However, the data requirements proved too onerous given the time frame.
Area of interest
The area chosen to pilot this method is the Massachusetts coast, including all the towns with coastlines and an arbitrary expansion of the neighboring towns inland. The study area is delineated by towns so that future analysis could include town-level data, such as population growth, fiscal health, and zoning restrictions. The area is bounded by the area within 1974758.27198, 2290559.08483, 2138987.38199, 2496151.08937. The figure below depicts the most recent land use and land cover data over the study area
lab <- figs(name="studyArea","Land Use Classifications on the Massachusetts Coast in 1992")
levelplot(nlcd.2011, main=lab, col.regions=color.theme)
Methods
The approach has two steps. First, I develop a minimal Bayesian Network (BN) model of land use the study area, using only the data from 1992, 2001, and 2006. Land use transition tables between periods for each land use classification are constructed using these observations. This model is used to predict the land use classification in 2011, which is then compared to the observed data from 2011. This first step establishes the accuracy of the BN model in replicating observed data. The overall accuracy is greater than 98% for all land uses. I also apply the comparison method proposed by Pontius, Shusas, and McEachern (???) in a novel way to assess the accuracy of the model to predict gain, loss, swap, and net change for each of the land use classifications.
The second step is to complicate the BN model by introducing scenarios. To do this, I incorporate all the available data from 1992, 2001, 2006, and 2011, and generate three scenarios - development, conservation, and business-as-usual (BAU). Each of the scenarios are developed by applying a multiplier to the transition table from 2006-2011 to generate fictitious transition tables. Applying the three different scenarios, the BN model produces three possible outcomes for the 2016 time period, for which data are not yet available from the same data source (described below).
The following section describes the data used to develop the model and the necessary transformations to format the data for consumption by the modeling software.
Description of spatial data needed
While the literature review introduced a number of potential covariates, these have been omitted from this minimal model. All data come from a single source, the National Land use and Land Cover Database, which is described below.
USGS National Landuse Classification Dataset
The National Landuse Classification Dataset (NLCD) is constructed by the United States Geological Survey (USGS) using a classification algorithm applied to Landsat data. The table below shows the classifications in the processed data as downloaded from the website (Homer et al. 2015). To reduce the complexity of the model and coding procedures needed to produce a minimal example, the classifications were reduced to five classes - water, developed, undeveloped, agricultural, and wetland. The NLCD classifications in 1992 are slightly different from the classifications from subsequent years, which causes some classifications to be inconsistent. However, these thematic classifications are closely comparable and assumed insignificant to the model. In the below table, the original (1992) classifications are shown next to their update (2001-2011) classifications, along with their descriptions. The two right-most columns show how these classifications are reclassified into the model. Further below is a description and code example of how to read in and reclassify the data. The data for this paper were simplified in GRASS GIS, though it could also be done in R using a similar procedure to that shown further below.
The data were reprojected into universal trans mercator projection using zone 18 for MAssachusetts. The reprojection was also performed outside of R due to the personal preference of the author, but could be performed in R. Once read into R using the raster package, the following routine was performed for each of the NLCD data sets. For simplification, only the 1992 layer is shown in the code chunk below. The routine includes one step that is not necessar for the other data sets. This was done to illustrate how to reclassify in R using a simple example. The 1992 data layer contained some “transitional” classifications that were not present in other years’ data. These were converted to “undeveloped” using the third line of code below.
# Read in the land use files and change the variables to factors.
nlcd.1992 <- raster("lu_1992.tif")
# To avoid projecting the model over the entire bounding box, mask the areas out of range of the categories. The land use categories range from 10-90, so nothing over 100 should be displayed.
nlcd.1992 <- mask(nlcd.1992, nlcd.1992<100)
# Reclassify 30 (transitional) as 40 (undeveloped).
nlcd.1992[nlcd.1992==30] <- 40
# Create a raster attribute table (RAT) with the `ratify` function.
nlcd.1992 <- ratify(nlcd.1992)
# Modify the RAT with the labels
l <- levels(nlcd.1992)[[1]]
l$code <- c("Water.1992","Developed.1992","Undeveloped.1992","Agricultural.1992","Wetland.1992")
# Assign the modified RAT to the raster with the levels function.
levels(nlcd.1992) <- l
nlcd.stack <- stack(nlcd.1992, nlcd.2001, nlcd.2006, nlcd.2011)A simple table of the percent changes in land cover for each land classification is shown in the table below. Only developed land uses increase in each time period, though the rate appears to be slowing from 2001-2006 interval to the 2006-2011 interval. Note that the time series intervals are not consistent, with a 9 year gap followed by two 5 year gaps.
lab <- tbls(name="luChngRates", caption="Landuse Change Rates between Periods")
knitr::kable(lu.chng.rates, caption=lab) %>%
kable_styling(latex_options = c("striped", "hold_position"),
full_width = T)| Chng.1992.to.2001 | Chng.2001.to.2006 | Chng.2006.to.2011 | |
|---|---|---|---|
| Water | -0.12 | -0.05 | 0.01 |
| Developed | 6.66 | 1.66 | 0.75 |
| Undeveloped | -12.26 | -1.15 | -0.52 |
| Agricultural | 2.93 | -0.25 | -0.20 |
| Wetland | 2.80 | -0.23 | -0.04 |
Building a Minimal Model
The building a Bayesian belief network model requires constructing multiple intermediate parts, and two R packages have been helpful for doing this. In all, three key packages will organize this section, as a step-by-step guide for implementing a minimal model, lulcc, gRain, and bnspatial.
lulcc Package
The lulcc package is designed to analyze land use and land cover change, which give it its name. In this example, only one key feature is used, the crossTabulate function. However, the lulcc package also offers functions for generalizes linear model, random forest modeling of land use changes over time. Other functions in the lulcc package are designed to format a series of raster in the same projection and extent into observation data and covariates. Observations must have the same classifications. To estimate the probabilities for each of the relationships in the Bayes net, the existing data provide the transition probabilities. A simple cross table for each land use over each interval can suffice to estimate the ‘likelihood’ that any individual land use will transition to any other. The dataframes created here are the transition tables for ending year of the variable name. Note that the conditional probability tables must be converted to matrices in R before they can be input to the bnspatial package as described later in the section. The final variable, ctable.df.total contains the total land use change from 1992 to 2011. In the code below, the cross tables count the number of pixels that transitioned between states from one time period to the next. Some of the values had zero transitions, which causes errors in the bnspatial package. The work-around shown here is to add a small value (in this case, 1) to each value in the cross table. This will minimally change the ratios, but the values are large enough in this case that the overall accuracy of the model was not affected.
# LULCC method
# create a list of the RasterLayers that are covariates. This is not done here, but could be easily added by adding layers beginning with the charaters "ef" to the nlcd.stack.
#ef <- ExpVarRasterList(x = nlcd.stack, pattern="ef")
# create a list of the RasterLayers that are observations
obs <- ObsLulcRasterStack(
nlcd.stack[[1:4]],
categories=c(10,20,40,80,90),
labels=c("Water","Developed","Undeveloped","Agricultural","Wetland"),
t=c(0,9,14,19))
# Not used, but notable functions can be run as follows. See the help for lulcc package for more.
#part <- partition(x=obs[[1]], size=0.1, spatial=TRUE)
#train.data <- getPredictiveModelInputData(obs=obs, ef=ef, cells=part[["train"]], t=0)
#
#forms <- list(Developed~ef.002,
# Undeveloped~ef.002,
# Water~ef.002,
# Wetland~ef.002,
# Agricultural~ef.002)
#glm.models <- glmModels(formula=forms, family=binomial, data=train.data, obs=obs)
#rpart.models <- rpartModels(formula=forms,data=train.data,obs=obs)
#rf.models <- randomForestModels(formula=forms,data=train.data,obs=obs)According to Pontius, Shusas, and McEachern (2004), landscape changes should be measured in terms of the gain and loss each category, as well as with the swap and net change. The function to calculate swap is \(swap = 2 * min(Gain, Loss)\), or the amount of the land class for which loss and gain cancel out (Pontius, Shusas, and McEachern 2004). The table below shows the gain, loss, and between-group transitions analysis for the entire period.
The cross tables used for conditional probability are shown below. The conditional probabilities are only the rows and columns that correspond to specific land classifications. The total, loss, and gain values are not part of the conditional probability tables, but they are useful in getting a general idea of the rates of change in each period for each classification. The values in the first table show the percentage of the land area that transitioned from the classification in the row header in 1992 to the value in the column header in 2001. For example, 0.07% of developed land in 1992 was classified as Water in 2001. The diagonal is the percent of the land scape that persisted. For example, 27.13% out of a total of 38.88% of the 1992 developed land remained developed in 2001. Out of the total 11.752524% of the developed land that changed (persistent - total), 11.752524% was gained and 6.2911501% was lost.
The table below shows the gain, loss, swap, and net change over the entire period.
gRain Package
Now that we have seen the capabilities of the lulcc package, I turn to demonstrate the process for building a Bayesian network using the gRain package. The network used in this table is actually linear, in that each year is a node that only depends on the previous year. A the values= argument takes a matrix or vector that is of equal length to the number of classes in the predicted node (signified with the ~ symbol) and multipied by the number classes in each predictor node (which follow the | symbol). The example further below shows the format for a node with multiple dependencies, in which conditional probabilities for the predictor listed last should be first among the values= argument. See below for an example.
# create a vector with the land use classes IN ORDER
lu <- c("Water","Developed","Undeveloped","Agricultural","Wetland")
# use ctable matrices as conditional probability values
lu_1992 <- cptable(~lu_1992, values=c(lu.area.1992[,2]), levels=lu)
lu_2001 <- cptable(~lu_2001|lu_1992, values=c(as.matrix(ctable.df.2001[1:5,1:5])), levels=lu)
lu_2006 <- cptable(~lu_2006|lu_2001, values=c(as.matrix(ctable.df.2006[1:5,1:5])), levels=lu)
lu_2011 <- cptable(~lu_2011|lu_2006, values=c(as.matrix(ctable.df.2011[1:5,1:5])), levels=lu)Once the nodes are established, the compileCPT function is used to create a conditional probability table that considers all of the nodes. In order to graph this as a network and load it into the bnspatial package, the conditional probability table has to be loaded as anetwork into the environment with the loadNetwork() function. After this is done, we can visually check the network using the Rgraphviz::plotI() function.
# compile the nodes into a network that can be loaded into the environment
lunet1 <- compileCPT(list(lu_2011, lu_2006, lu_2001, lu_1992))
# load into the environment the network as a grain object
lunNet.gr <- loadNetwork(grain(lunet1), target="lu_2011") # load the network as an object in the environment
network <- lunNet.gr
# plot the network
lab <- figs(name="lunNet.gr", caption="Simple Bayesian Network")
plot(lunNet.gr, main=lab)
bnspatial Package
The network is clearly very simple, with each node dependent only on the previous node. However, as we will see, this is still a very good predictor of the 2011 land use, which we can test by mapping the most likely classes for 2011 based on the previous three observations and then comparing the results. This can be done with the bnspatial package by creating a bnspatial object with lu_2011 as the target. The code chunk below shows the process for creating a bnspatial object.
# set the target
target <- 'lu_2011'
# make a vector of the data without the target (2011 lan use layer)
spatialData <- c(nlcd.stack$lu_1992,
nlcd.stack$lu_2001,
nlcd.stack$lu_2006)
# nlcd.stack$lu_2011) Why don't we remove this, since it is predicted?
#paste0('lu_2006.tif'),
#paste0('lu_2011.tif'))
lookup <- paste0("luClasses.txt")
# An alternative way to do this is through the setClasses function as below.
#setClasses(c('lu_1992','lu_2001', 'lu_2006','lu_2011'),
# list(c("Water","Developed", "Undeveloped","Agriculture","Wetland"),
# c("Water","Developed", "Undeveloped","Agriculture","Wetland"),
# c("Water","Developed", "Undeveloped","Agriculture","Wetland"),
# c("Water","Developed","Undeveloped","Agriculture","Wetland")),
# list(c(lu.area.1992[,2]), c(lu.area.2001[,2])))
## Run a spatial query on the Bayesian network
bn <- bnspatial(network, target, spatialData , lookup, msk=(nlcd.stack$lu_2011<100))##
## "lu_1992" points to:
## -> lu_1992
##
## With states:
## Water Developed Undeveloped Agricultural Wetland
##
## Represented by integer values:
## 10 <-> 20 <-> 40 <-> 80 <-> 90
## ----------------------------------
##
## "lu_2001" points to:
## -> lu_2001
##
## With states:
## Water Developed Undeveloped Agricultural Wetland
##
## Represented by integer values:
## 10 <-> 20 <-> 40 <-> 80 <-> 90
## ----------------------------------
##
## "lu_2006" points to:
## -> lu_2006
##
## With states:
## Water Developed Undeveloped Agricultural Wetland
##
## Represented by integer values:
## 10 <-> 20 <-> 40 <-> 80 <-> 90
## ----------------------------------
## Lookup table to interpret "lu_2011" values:
## lu_2011 cell_ID
## 1 Water 1
## 2 Developed 2
## 3 Undeveloped 3
## 4 Agricultural 4
## 5 Wetland 5By default, the output of the bnspatial() function shows all of the nodes and their possible states. This is a good time to check that the intended model has been achieved. To compare the results with the known classifications of the 2011 layer. As shown in the table below, the highest probability classifications produced through the Bayesian Belief Network is over 98% accurate at predicting all five classes of land use. The diagonal shows the accurate prediction percentages for each class. The cells off-diagonal show inaccurate highed probabilities and the errant classifications. The column headers refer to the predicted classifications and the row headers refer to the actual 2011 classifications.
Results
The table below shows the accuracy of the predicted 2011 land cover map compared to the observed land uses. The table shows the predicted classes in the columns and the actual classes in the rows in percentage of the class total. The diagonal are accurate predictions. As mentioned above, this test is artificially high, since most pixels are likely to persist in their previous state (Pontius, Shusas, and McEachern 2004).
lab <- tbls(name="accuracy.pct", caption="Accuracy of BN model (rows) to Observed 2011 Land Uses (columns)")
knitr::kable(accuracy.pct, caption=lab)| Water | Developed | Undeveloped | Agricultural | Wetland | NA | |
|---|---|---|---|---|---|---|
| Water | 99.91 | 0.00 | 0.00 | 0 | 0 | 0.00 |
| Developed | 0.00 | 99.98 | 0.00 | 0 | 0 | 0.00 |
| Undeveloped | 0.00 | 0.00 | 99.86 | 0 | 0 | 0.03 |
| Agricultural | 0.00 | 0.00 | 0.00 | 100 | 0 | 0.00 |
| Wetland | 0.00 | 0.00 | 0.00 | 0 | 100 | 0.00 |
| NA | 0.09 | 0.02 | 0.14 | 0 | 0 | 99.97 |
| Total | 100.00 | 100.00 | 100.00 | 100 | 100 | 100.00 |
Using this assessment alone, a model predicting no change would be highly accurate. However, if we examine the gain, loss, swap, and net change, a more detailed picture emerges. Figure 3 below allows comparison between the observed transitions and predicted transitions, which reveals that the model systematically underpredicts transitions.
The prediction is more conservative with respect to transitions. There are fewer transitions inthe predicted model than in the actual land use changes, particularly for developed land. However, fairly limited change occured over the period, and the model reflects this.
Furthermore, the spatial distribution of the false classifications, shown as white pixels in the right panel below, appear to be random. The map series below shows the highest probability, the actual NLCD classifications, and the locations of accurate predictions.
Another feature of the bnsptial package is to generate a probability for each state of the target node, which can be useful to display maps of the various probabilities. This is done by adding the what="probability" argument, and the resulting bnspatial object will then contain a separate layer for the probabilities of each state, in this case the land use classification.
Developing Scenarios
Now that a minimal example has been shown and evaluated, this section will demonstrate the power of the bnspatial package to develop scenarios. First, we have to create two new nodes that contain the probabilty estimates for the future land use and under different scenarios. We also must create a new network using these new probability tables. For simplicity, the probability of transition for the future land uses, here called lu_current, are the same as the 2006-2011 probabilities. Three scenarios are developed - BAU (business as usual), Development, and Conservation. BAU transition probabilities are the same as the probabilities for 2006-2011. Development scenario transition probabilities are higher for all land uses to transition to developed or agricultural land, and lower for land to transition from development to any other classifications. Conservation scenario entails that land is less likely than in the 2006-2011 time period to transition to developed, and more likely that land will be reclaimed from developed to undeveloped. Values in the diagonal of each scenario are set to 1, since there is not a higher likelihood that, for example, agricultural land will remain agrcultural land. The matrices below, cons.matrix and dev.matrix show the multipliers that are applied to the cross table generated using crossTabulate as shown above. Note that this is not using matrix algebra, but simply multiplying the corresponding values in each position in both matrices, the scenario matrices and the probability matrices calculated above. This format makes it easier to calibrate the scenarios without calculating new probability tables.
# create a matrix for each scenario. This will be multiplied by the cross table to get the scenario transition probabilities.
cons.matrix <- matrix(nrow = nrow(ctable.df.total), ncol = ncol(ctable.df.total),
data = c(1.00, 0.80, 1.00, 0.80, 1.00,
0.80, 1.00, 0.95, 0.95, 1.00,
1.00, 1.10, 1.00, 0.90, 1.00,
0.80, 1.10, 1.10, 1.00, 1.00,
1.00, 1.00, 1.00, 1.00, 1.00))
dev.matrix <- matrix(nrow = nrow(ctable.df.total), ncol = ncol(ctable.df.total),
data = c(1.0, 1.2, 1.0, 1.1, 1.0,
1.2, 1.0, 1.2, 1.2, 1.2,
1.0, 1.2, 1.0, 1.1, 1.0,
1.1, 1.2, 1.1, 1.0, 1.1,
1.0, 1.2, 1.0, 1.1, 1.0))
scenarios <- cptable(~scenarios, values = c(3,3,3), levels = c("BAU","Conservation","Development"))
lu_current <- cptable(~lu_current|scenarios:lu_2011,
values=c(as.matrix(ctable.df.2011[1:5,1:5]),
as.matrix(ctable.df.2011[1:5,1:5]*cons.matrix),
as.matrix(ctable.df.2011[1:5,1:5]*dev.matrix)), # Wetland ctable.df.2011[5,]),
levels=lu)As the values have changed, the next step is to calculate another network and bnspatial object to compare the scenarios.
To compare the most likely class for each scenario, we can run the bnspatial() function for each of the scenarios and maps them separately.
## Run a spatial query on the Bayesian network
bn.s.dev <- bnspatial(network, target, spatialData , lookup, msk=(nlcd.stack$lu_2011<100), verbose = F,
scenarios="Development")## Lookup table to interpret "lu_current" values:
## lu_current cell_ID
## 1 Water 1
## 2 Developed 2
## 3 Undeveloped 3
## 4 Agricultural 4
## 5 Wetland 5bn.s.cons <- bnspatial(network, target, spatialData , lookup, msk=(nlcd.stack$lu_2011<100), verbose = F,
scenarios="Conservation")## Lookup table to interpret "lu_current" values:
## lu_current cell_ID
## 1 Water 1
## 2 Developed 2
## 3 Undeveloped 3
## 4 Agricultural 4
## 5 Wetland 5bn.s.bau <- bnspatial(network, target, spatialData , lookup, msk=(nlcd.stack$lu_2011<100), verbose = F,
scenarios="BAU")## Lookup table to interpret "lu_current" values:
## lu_current cell_ID
## 1 Water 1
## 2 Developed 2
## 3 Undeveloped 3
## 4 Agricultural 4
## 5 Wetland 5bn.s.dev.prob <- bnspatial(network, target, spatialData , lookup, msk=(nlcd.stack$lu_2011<100),
what="probability", verbose = F,
targetState = c("Developed","Agricultural","Undeveloped"),
scenarios="Development" )par(mfrow=c(1,3))
levelplot(bn.s.dev$Class, main="Development Scenario", margin=F, col.regions=color.theme)
levelplot(bn.s.cons$Class, main="Conservation Scenario", margin=F, col.regions=color.theme)
levelplot(bn.s.bau$Class, main="BAU Scenario", col.regions=color.theme)
The scenarios as first estimated cause an unrealistic number of cells to transition from development to water. This shows that the estimated land use transition multipliers should be much smaller. Given the format of the data, this can be easily done by replicating the code that produced the cons.matrix and modifying the numbers so they are smaller. This is done in a hidden chunk for ease of reading. The results mapped below show a more realistic projection.
# create a matrix for each scenario. This will be multiplied by the cross table to get the scenario transition probabilities.
cons.matrix <- matrix(nrow = nrow(ctable.df.total), ncol = ncol(ctable.df.total),
#Wat #Dev #Undev #Ag #Wet
data = c(1.0, 1.0, 1, 1, 1.00, #Wat
1.0, 1.0, 1.00001, 1, 1, #Dev
1.0, 0.999999, 1, 1, 1.00, #Undev
1.0, 0.999999, 1, 1.00, 1.00, #Ag
1.0, 0.999999, 1.00, 1.00, 1.00))#Wet
dev.matrix <- matrix(nrow = nrow(ctable.df.total), ncol = ncol(ctable.df.total),
#Wat #Dev #Unde#Ag #Wet
data = c(1.0, 1.0000001, 1.0, 1.0, 1.0, #Wat
1.0, 1.0000001, 1.0, 1.0, 1.0, #Dev
1.0, 1.0000001, 1.0, 1.0, 1.0, #Undev
1.0, 1.0000001, 1.0, 1.0, 1.0, #Ag
1.0, 1.0000001, 1.0, 1.0, 1.0))#Wet
scenarios <- cptable(~scenarios, values = c(3,3,3), levels = c("BAU","Conservation","Development"))
lu_current <- cptable(~lu_current|scenarios:lu_2011,
values=c(as.matrix(ctable.df.2011[1:5,1:5]),
as.matrix(ctable.df.2011[1:5,1:5]*cons.matrix),
as.matrix(ctable.df.2011[1:5,1:5]*dev.matrix)), # Wetland ctable.df.2011[5,]),
levels=lu)
# compile the nodes into a network that can be loaded into the environment
lunet2 <- compileCPT(list(lu_current,lu_2011, lu_2006, lu_2001, lu_1992, scenarios))
# load into the environment the network as a grain object
lunNet.gr <- loadNetwork(grain(lunet2), target="lu_current") # load the network as an object in the environment
network <- lunNet.gr
# plot the network
#plot(lunNet.gr)
# set the target
target <- 'lu_current'
# make a vector of the data without the target (2011 lan use layer)
spatialData <- c(nlcd.stack$lu_1992,
nlcd.stack$lu_2001,
nlcd.stack$lu_2006,
nlcd.stack$lu_2011) # This is now added to the spatialData object.
lookup <- paste0("luClasses2.txt")In the plots below, the Development scenario converts much of the developed land to water, and leaves much undeveloped land undeveloped. This claerly needs calibration and perhaps adding a more complex model would contribute to that.
par(mfrow=c(1,3))
levelplot(bn.s.dev$Class, main="Development Scenario", margin=F)
levelplot(bn.s.cons$Class, main="Conservation Scenario", margin=F)
levelplot(bn.s.bau$Class, main="BAU Scenario", margin=F)
The maps below show the probability for land to be developed, undeveloped, and agricultural in the development scenario. Clearly agricultural land is the likely to have the biggest gains under this scenario. Again, this result is not the intended result, and further calibration of the model is needed.
bn.s.dev.prob <- bnspatial(network, target, spatialData , lookup, msk=(nlcd.stack$lu_2011<100), verbose = F,
what="probability", targetState=c("Developed", "Undeveloped","Agricultural"),
scenarios="Development")
par(mfrow=c(1,3))
plot(bn.s.dev.prob$Probability$Developed, main="Developed")
plot(bn.s.dev.prob$Probability$Undeveloped, main="Undeveloped")
plot(bn.s.dev.prob$Probability$Agricultural, main="Agricultural")
Even with the modified scenario tables, the development scenario is unrealistic, with the majority of the developed land transitioning to water. It is not clear why the development scenario converts most of the developed land to water. This will have to be resolved later.
References
Burby, Raymond J. 2003. “Making Plans That Matter: Citizen Involvement and Government Action.” Journal of the American Planning Association 69 (1): 33–49. doi:10.1080/01944360308976292.
Celio, Enrico, Thomas Koellner, and Adrienne Grêt-Regamey. 2014. “Modeling Land Use Decisions with Bayesian Networks: Spatially Explicit Analysis of Driving Forces on Land Use Change.” Environmental Modelling & Software 52 (February): 222–33. doi:10.1016/j.envsoft.2013.10.014.
Clarke, Keith C. 2008. “Mapping and Modelling Land Use Change: An Application of the SLEUTH Model.” In Landscape Analysis and Visualisation, 353–66. Lecture Notes in Geoinformation and Cartography. Springer, Berlin, Heidelberg. doi:10.1007/978-3-540-69168-6_17.
Hashiguchi, Yoshihiro, and Kiyoyasu Tanaka. 2015. “Agglomeration and Firm-Level Productivity: A Bayesian Spatial Approach: Agglomeration and Firm-Level Productivity.” Papers in Regional Science 94 (November): S95–S114. doi:10.1111/pirs.12147.
Homer, Colin, Jon Dewitz, Yang Limin, Suming Jin, Patrick Danielson, George Xian, John Coulston, Nathaniel Herold, James Wickham, and Kevin Megown. 2015. “Completion of the 2011 National Land Cover Database for the Conterminous United States – Representing a Decade of Land Cover Change Information.” Photogrammetric Engineering & Remote Sensing 81 (5): 345–54. https://cfpub.epa.gov/si/si_public_record_report.cfm?dirEntryId=309950.
Hua, Lizhong, Lina Tang, Shenghui Cui, and Kai Yin. 2014. “Simulating Urban Growth Using the SLEUTH Model in a Coastal Peri-Urban District in China.” Sustainability 6 (6): 3899–3914. doi:10.3390/su6063899.
Park, Soyoung, Keith C. Clarke, Chuluong Choi, and Jinsoo Kim. 2017. “Simulating Land Use Change in the Seoul Metropolitan Area After Greenbelt Elimination Using the SLEUTH Model.” Journal of Sensors, 4012929. doi:10.1155/2017/4012929.
Pierce, S. M., R. M. Cowling, A. T. Knight, A. T. Lombard, M. Rouget, and T. Wolf. 2005. “Systematic Conservation Planning Products for Land-Use Planning: Interpretation for Implementation.” Biological Conservation 125 (4): 441–58. doi:10.1016/j.biocon.2005.04.019.
Pontius, Robert G., Emily Shusas, and Menzie McEachern. 2004. “Detecting Important Categorical Land Changes While Accounting for Persistence.” Agriculture, Ecosystems & Environment 101 (2-3): 251–68. doi:10.1016/j.agee.2003.09.008.
Schmidt, Stephan, Stefan Siedentop, and Stefan Fina. 2017. “How Effective Are Regions in Determining Urban Spatial Patterns? Evidence from Germany.” Journal of Urban Affairs 0 (0): 1–18. doi:10.1080/07352166.2017.1360741.
Seto, Karen C., Michail Fragkias, Burak Güneralp, and Michael K. Reilly. 2011. “A Meta-Analysis of Global Urban Land Expansion.” PLOS ONE 6 (8): e23777. doi:10.1371/journal.pone.0023777.
Silva, E.A, and K.C Clarke. 2002. “Calibration of the SLEUTH Urban Growth Model for Lisbon and Porto, Portugal.” Computers, Environment and Urban Systems 26 (6): 525–52. doi:10.1016/S0198-9715(01)00014-X.
Torrens, Paul. 2000. “How Land-Use-Transportation Models Work.” Paper 20. London: University College London.
Zhang, Qian, Jeremy Wallace, Xiangzheng Deng, and Karen C. Seto. 2014. “Central Versus Local States: Which Matters More in Affecting China’s Urban Growth?” Land Use Policy 38 (May): 487–96. doi:10.1016/j.landusepol.2013.12.015.