Here, we’ll start with no data and build a simple CA model with a single transition rule using a spatial Markov framework. That is, the probability of a unit transitioning into different uses is a probabilistic function of its current state, and the states of the units around it. Using data from NLCD, we can observe how often these transitions occurred in the past, including how those transitions are conditioned by different spatial contexts. Then we can use that model to simulate conditions into the future.
26.1 Data
Start by grabbing tract-level data for the Detroit-ish region
Make this Notebook Trusted to load map: File -> Trust Notebook
Now, create a regular hexagonal grid that covers the same area and convert the land cover data from NLCD from a raster into a vector, collapsing contiguous cells into regions
Now aggregate the pixels into the (coarser) hexgrid and include a column of colors to match the conventions in land use mapping
Code
def img_to_hex(hexes, img_path, year): hexes = hexes.copy()# take the LULC image and convert from raster to vector lulc = extract_raster_features(hexes, img_path)# convert to a string for easier representation as categorical variable lulc["value"] = lulc["value"].astype(int).astype(str)# pixel value 255 is NA lulc["value"] = lulc["value"].replace({"255": np.nan}) temp = gpd.sjoin(lulc.to_crs(hexes.crs), hexes) temp = ( temp.groupby("hex_id")["value"] .agg(lambda x: pd.Series.mode(x)[0]) .astype(str) ) # cheat and take the min when multiple modes are returned... hexes = hexes.join(temp) hexes["year"] = year hexes["value"] = hexes["value"].replace(lu_mapper).replace({"nan": np.nan}) hexes = hexes.dropna()return hexes
Code
nlcd_colors = pd.read_csv('nlcd_colors.csv',)
Code
# convert rgb to 0-1 rangenlcd_colors.color = nlcd_colors.color.apply(lambda x: np.array(x.split(",")).astype(float) /255)# convert rgb to hex (easier to hold a single value in a column instead of a listnlcd_colors.color = nlcd_colors.color.apply(lambda x: to_hex(x))
Using geosnap and PySAL, we can look at these transitions more formally, modeling the change over time as a spatial Markov process. Using this framework, we can ask whether certain transitions are more likely than others given their spatial context. For example, is a piece of undeveloped land more likely to become developed if it’s neighbors have already been developed?
Code
sm = transition(lulc_hexes, unit_index='hex_id', cluster_col='value')
/Users/knaaptime/Dropbox/projects/geosnap/geosnap/analyze/dynamics.py:123: FutureWarning: `use_index` defaults to False but will default to True in future. Set True/False directly to control this behavior and silence this warning
w = Ws[w_type].from_dataframe(gpd.GeoDataFrame(gdf_wide), **w_options)
/Users/knaaptime/miniforge3/envs/urban_analysis/lib/python3.12/site-packages/libpysal/weights/contiguity.py:61: UserWarning: The weights matrix is not fully connected:
There are 6 disconnected components.
There are 2 islands with ids: 16374, 38660.
W.__init__(self, neighbors, ids=ids, **kw)
('WARNING: ', 16374, ' is an island (no neighbors)')
('WARNING: ', 38660, ' is an island (no neighbors)')
/Users/knaaptime/miniforge3/envs/urban_analysis/lib/python3.12/site-packages/numpy/core/fromnumeric.py:3504: RuntimeWarning: Mean of empty slice.
return _methods._mean(a, axis=axis, dtype=dtype,
Code
type(sm)
giddy.markov.Spatial_Markov
Inference from the spatial Markov object shows the transitions are not random, but show different rates depending on the “neighborhood” of each observation
Code
sm.Q_p_value
0.0
Code
sm.LR_p_value
0.0
Plotting all the transition matrices can be daunting, but it provides a detailed look at how transition rates vary by context
We can see that barren land, for example, is much more likely to become developed when it’s neighbors are low-intensity developed (from top-left, see the second row, second column)
26.3 Simulating Future Land-Use
With these transition probabilities in hand (and the set of land uses in the last time period), we have all we need to simulate potential land uses in the future. First, we assume that the empirical transitions observed in prior time periods will remain static, and that only the prior time period influences the probability of transition between landuse types. Then, we can use the spatial context of each unit to select the proper transition matrix, and draw a new land-use type using the transition probabilities observed in the past.
Repeating this process for each unit in the study area gives a stochastic realization of land-use transitions; repeating this process for a range of years lets us look at a potential evolution of the region into the future. We could then change transition probabilities or manually edit land-uses in certain places (e.g. to simulate processes like urban development, sea-level rise, or a protective land-use policy) and see how those changes propagate through the system into the future.
Code
w = Queen.from_dataframe(lulc_hexes[lulc_hexes.year==2016].dropna().reset_index())
/var/folders/j8/5bgcw6hs7cqcbbz48d6bsftw0000gp/T/ipykernel_82155/4138579200.py:1: FutureWarning: `use_index` defaults to False but will default to True in future. Set True/False directly to control this behavior and silence this warning
w = Queen.from_dataframe(lulc_hexes[lulc_hexes.year==2016].dropna().reset_index())
/Users/knaaptime/miniforge3/envs/urban_analysis/lib/python3.12/site-packages/libpysal/weights/contiguity.py:347: UserWarning: The weights matrix is not fully connected:
There are 5 disconnected components.
There is 1 island with id: 34415.
W.__init__(self, neighbors, ids=ids, **kw)
/Users/knaaptime/Dropbox/projects/geosnap/geosnap/visualize/mapping.py:374: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`). Consider using `matplotlib.pyplot.close()`.
fig, ax = plt.subplots(figsize=figsize)
<Figure size 1152x1152 with 0 Axes>
<Figure size 1152x1152 with 0 Axes>
<Figure size 1152x1152 with 0 Axes>
<Figure size 1152x1152 with 0 Axes>
<Figure size 1152x1152 with 0 Axes>
<Figure size 1152x1152 with 0 Axes>
<Figure size 1152x1152 with 0 Axes>
<Figure size 1152x1152 with 0 Axes>
<Figure size 1152x1152 with 0 Axes>
<Figure size 1152x1152 with 0 Axes>
<Figure size 1152x1152 with 0 Axes>
<Figure size 1152x1152 with 0 Axes>
<Figure size 1152x1152 with 0 Axes>
<Figure size 1152x1152 with 0 Axes>
<Figure size 1152x1152 with 0 Axes>
<Figure size 1152x1152 with 0 Axes>
<Figure size 1152x1152 with 0 Axes>
<Figure size 1152x1152 with 0 Axes>
<Figure size 1152x1152 with 0 Axes>
<Figure size 1152x1152 with 0 Axes>
<Figure size 1152x1152 with 0 Axes>
<Figure size 1152x1152 with 0 Axes>
<Figure size 1152x1152 with 0 Axes>
<Figure size 1152x1152 with 0 Axes>
<Figure size 1152x1152 with 0 Axes>
<Figure size 1152x1152 with 0 Axes>
<Figure size 960x960 with 0 Axes>
Now we can plot the simulated land-use moving into the future. This model is pretty naive, but the results seem plausible, with development growing outward from existing cores. We also get lots of infill (but the model isnt smart enough to know some places can’t be developed, etc). It would be fairly easy to include some additional restrictions that some observations are fixed because of policy or topology
We can also quickly compute the share of land in each type and compare two time periods. Since the output is a regular grid, we just need to count the number of cells in each category without needing to actually look at the area
value
Barren Land (Rock/Sand/Clay) -0.001453
Cultivated Crops -0.007623
Deciduous Forest -0.024469
Developed High Intensity 0.006696
Developed, Low Intensity -0.002695
Developed, Medium Intensity 0.049971
Developed, Open Space 0.000758
Emergent Herbaceous Wetlands -0.002548
Evergreen Forest -0.000147
Grassland/Herbaceous -0.001369
Mixed Forest -0.002738
Pasture/Hay -0.006317
Woody Wetlands -0.008065
Name: proportion, dtype: float64
This would suggest a continued trend toward urbanization, with particularly large relative losses in deciduous forest. Low-intensity development is decreasing, but medium-intensity seems to be growing strongly. Following this analysis, some natural questions that arise include, what would it mean (in social, economic, environmental, etc.) terms to have a 5% increase in medium-intensity development and a 2.5% reduction in deciduous forest? What are the implications for climate change, anticipated housing growth, traffic congestion, change in tax revenue, and so on? And how might these conditions differ under a different scenario, such as one where we have alternative land regulations or urban growth constraints?
Brail, R. K. (Ed.). (2008). Planning support systems for cities and regions. Lincoln Institute of Land Policy.
Chakraborty, A., Kaza, N., Knaap, G.-J., & Deal, B. (2011). Robust Plans and Contingent Plans: Scenario Planning for an Uncertain World. Journal of the American Planning Association, 77(3), 251–266. https://doi.org/10.1080/01944363.2011.582394
Chakraborty, A., & McMillan, A. (2015). Scenario Planning for Urban Planners: Toward a Practitioner’s Guide. Journal of the American Planning Association, 81(1), 18–29. https://doi.org/10.1080/01944363.2015.1038576
Dietzel, C., & Clarke, K. (2006). The effect of disaggregating land use categories in cellular automata during model calibration and forecasting. Computers, Environment and Urban Systems, 30(1), 78–101. https://doi.org/10.1016/j.compenvurbsys.2005.04.001
Guan, Q., & Clarke, K. C. (2010). A general-purpose parallel raster processing programming library test application using a geographic cellular automata model. International Journal of Geographical Information Science, 24(5), 695–722. https://doi.org/10.1080/13658810902984228
Jantz, C. A., Goetz, S. J., Donato, D., & Claggett, P. (2010). Designing and implementing a regional urban modeling system using the SLEUTH cellular urban model. Computers, Environment and Urban Systems, 34(1), 1–16. https://doi.org/10.1016/j.compenvurbsys.2009.08.003
Knaap, G., Engelberg, D., Avin, U., Erdogan, S., Ducca, F., Welch, T. F., Finio, N., Moeckel, R., & Shahumyan, H. (2020). Modeling Sustainability Scenarios in the Baltimore–Washington (DC) Region. Journal of the American Planning Association, 0(0), 1–14. https://doi.org/10.1080/01944363.2019.1680311
