Modeling the Growth of Cities Using Satellite Images
Please cite the following papers if you feel that your research benefits from the codes and ideas in this chapter:
Wu, L. & Wang, C. (2015). Tracing the Attention of Moving Citizens. Scientific reports, 6.
Satellite Images
Satellite images of earth at night are beautiful. They are not only appealing, but also provide data for analyzing and modeling the spatial distribution of human activities (e.g., wealth creation and energy consumption) in cities.
The center of Beijing
Reid Wiseman, an American astronaut, took the following photo of Beijing and put it on Twitter:
Can we find the downtown, or the Central Business District (CBD) of Beijing without any background knowledge, and only by analyzing the photo? It is reasonable to assume that downtown contains the central point of the city. But how to define the "central point" ?
Let's do some simple but fancy analysis using scipy.ndimage functions.
import numpy as np
import matplotlib.pyplot as plt
from scipy import ndimage
from numpy import linalg as LA
from scipy.ndimage import measurements
im = ndimage.imread('/Users/csid/Desktop/beijingnightmap.jpg',flatten=True)
mg = plt.imshow(im,cmap = plt.get_cmap('gray'))
y,x = ndimage.measurements.center_of_mass(im)
xlim = mg.axes.get_xlim()
ylim = mg.axes.get_ylim()
plt.xlim(xlim)
plt.ylim(ylim)
plt.imshow(im,cmap = plt.get_cmap('gray'))
plt.plot(x,y,'ro')
plt.show()
Using the above codes, I firstly convert the RGB image into a gray scale image, which is a 1000*1000 matrix whose elements vary from 0 (dark) to 1 (light), indicating the lightness of pixels.
Then I calculate the "center of mass" of pixels weighted by their lightness. If we assume that the lightness of pixels indicates the level of human activity in the corresponding location, then the "center of mass" is also the center point of the city.
The above figure shows that we can use this method to identify the center point (marked in red) of Beijing.
The downtown of Beijing
Now it would be interesting to analyze how human activity decays over distance from the center point of CBD in Beijing:
def OLSRegressFit(x,y):
xx = sm.add_constant(x, prepend=True)
res = sm.OLS(y,xx).fit()
constant, beta = res.params
r2 = res.rsquared
return [constant, beta, r2]
def linearBin(x,y):
a=[]
q=sorted(zip(map(int,x),y))
tag = ''
d=[]
for i in q:
x1,y1 = i
if x1 == tag:
d.append(y1)
else:
if tag:
a.append([tag,np.mean(d)])
tag = x1
d = []
a.append([tag,np.mean(d)])
nx,ny = np.array(a).T
return nx,ny
def add_subplot_axes(ax,rect,axisbg='w'):
fig = plt.gcf()
box = ax.get_position()
width = box.width
height = box.height
inax_position = ax.transAxes.transform(rect[0:2])
transFigure = fig.transFigure.inverted()
infig_position = transFigure.transform(inax_position)
x = infig_position[0]
y = infig_position[1]
width *= rect[2]
height *= rect[3] # <= Typo was here
subax = fig.add_axes([x,y,width,height],axisbg=axisbg)
x_labelsize = subax.get_xticklabels()[0].get_size()
y_labelsize = subax.get_yticklabels()[0].get_size()
x_labelsize *= rect[2]**0.5
y_labelsize *= rect[3]**0.5
subax.xaxis.set_tick_params(labelsize=x_labelsize)
subax.yaxis.set_tick_params(labelsize=y_labelsize)
return subax
n = 10000
indices = np.random.randint(0,len(im),size=(n,2))
gs = [im[x,y] for x,y in indices]
ds = [np.linalg.norm(i-np.array([x,y])) for i in indices]
nds,ngs=linearBin(ds,gs)
ml = im.max()
fig = plt.figure(figsize=(8, 8),facecolor='white')
ax = fig.add_subplot(1,1,1)
x,y = np.array([[x,y] for x,y in zip(nds,ngs/ml) if x>100 and y>0]).T
constant, beta, r2 = OLSRegressFit(np.log(x),y)
plt.plot(x,y,'b.')
plt.plot(x,constant+beta*np.log(x),'r-',label=r'$g\,=0.25-0.007log(r)$')
plt.xscale('log')
plt.xlabel(r'$r_i$',size=20)
plt.ylabel(r'$g_i$',size=20)
plt.legend(loc='lower left')
#
subax = add_subplot_axes(ax,[0.6,0.6,0.35,0.35])
subax.plot(ds,gs/ml,marker='.',linestyle='',color='RoyalBlue',alpha=0.03)
subax.plot(nds,ngs/ml,marker='.',linestyle='',color='orange',alpha=0.3)
subax.plot(x,constant+beta*np.log(x),'r-',label=r'$g\,=0.25-0.007log(r)$')
subax.set_xlim(1,2000)
subax.set_ylim(0,1)
subax.set_xscale('log')
plt.xlabel(r'$r_i$',size=16)
plt.ylabel(r'$g_i$',size=16)
#
plt.show()
It turns out that the decay of activity over distance follows a log function. There exists a "turning point" beyond which the decay dramatically accelerates. Before this point (x < 100), the level of human activity does not change a lot. This 100-pixel-radius area is the downtown of Beijing.
fig = plt.figure(figsize=(8, 8),facecolor='white')
ax = fig.add_subplot(1,1,1)
im = ndimage.imread('/Users/csid/Desktop/beijingnightmap.jpg',flatten=True)
mg = plt.imshow(im,cmap = plt.get_cmap('gray'))
y,x = ndimage.measurements.center_of_mass(im)
xlim = mg.axes.get_xlim()
ylim = mg.axes.get_ylim()
plt.xlim(xlim)
plt.ylim(ylim)
plt.imshow(im,cmap = plt.get_cmap('gray'))
plt.plot(x,y,'ro')
circle = plt.Circle((x, y), 100, color='gold',fill=True,alpha=0.5)
ax.add_artist(circle)
plt.show()
The above figure shows the downtown (marked in red) of Beijing.
The economics of scale of cities
Let's now turn to countries, which are systems containing many cities. From the NASA website we can download the following satellite image of the U.S.
And import it for analysis as:
im = ndimage.imread('.../Desktop/us.jpg',flatten=True)
im1 = im > 150
plt.imshow(im1,cmap = plt.get_cmap('gray'))
plt.show()
The above codes set a threshold equaling 150 (the maximum is 255) and turn the image into an adjacency matrix.
In the adjacency matrix as shown above, we can define a natural city area as a slice of the matrix in which all elements are connected. Using the following codes we detect 39,988 clusters from the U.S. image and plot them out.
lw, num = measurements.label(im1)
area = np.log(measurements.sum(im1, lw, index=np.arange(lw.max() + 1))+1)
areaImg = area[lw]
from mpl_toolkits.axes_grid1 import make_axes_locatable
plt.figure(figsize=(10, 8),facecolor='white')
ax = plt.gca()
ims = ax.imshow(areaImg, interpolation='nearest',cmap=plt.get_cmap('RdGy_r'))
# create an axes on the right side of ax. The width of cax will be 5%
# of ax and the padding between cax and ax will be fixed at 0.05 inch.
divider = make_axes_locatable(ax)
cax = divider.append_axes("right", size="2%", pad=0.3)
plt.colorbar(ims, cax=cax)
plt.show()
Obviously, the number of clusters depends on the threshold we choose. Therefore it is necessary to systematically investigate the relation between cluster numbers and the given threshold.
def nOfClusters(k,n):
nim = im > k
lw, num = measurements.label(nim)
area = measurements.sum(nim, lw, index=np.arange(lw.max() + 1))
return len([i for i in area if i>n])
a = [[i,nOfClusters(i,0)] for i in np.linspace(1,255,20)]
b = [[i,nOfClusters(i,3)] for i in np.linspace(1,255,20)]
c = [[i,nOfClusters(i,10)] for i in np.linspace(1,255,20)]
fig = plt.figure(figsize=(10, 8),facecolor='white')
x1,y1 = np.array(a).T
x2,y2 = np.array(b).T
x3,y3 = np.array(c).T
plt.plot(x1,y1,'bo-',label=r'$cluster\,size\,>0$')
plt.plot(x2,y2,'r^-',label=r'$cluster\,size\,>3$')
plt.plot(x3,y3,'gs-',label=r'$cluster\,size\,>10$')
plt.legend(loc='upper right',numpoints=1,fontsize=14)
plt.xlabel(r'$Threshold$',size=16)
plt.ylabel(r'$N \,of \,Clusters$',size=16)
plt.show()
It turns out that the number of clusters decays smoothly with the increase of threshold when threshold > 50. In other words, our choice of threshold would not change the result dramatically as long as it is greater than 50. In the current experiment, we use a threshold that equals 150.
We fin that the total activity within a city, measured by the sum of the lightness over all pixels within a cluster, increases super-linearly to city area, which is measured by the number of pixels in the cluster. This explains why city exists: when people get together, they create more interactions, wealth, and innovations than independently. A similar result was also discovered by Bettencourt et al. in 2006. They suggested that such a scaling function quantifies the economics of scale in cities.
Model city growth and human mobility
Inspired by the empirical observations from the satellite images, Zhang et al. proposed a "random geometric network" model to simulate the growth of cities and successfully replicated the super-linear scaling as well as the heterogeneous spatial distribution of activities. In another paper we extended this model and use it to explain the mobility of city residents.
This model is defined in a d-dimensional Euclidean space. Initially a seed node is placed at the center of the space. At each time step a new node is generated randomly within the space, added to the graph, and connected to all existing nodes within distance r. If there are no existing nodes within r, this new node will be removed.
We call this version “Model-all” and change the linking rule to develop two versions, including 1). the new nodes are connected to the highest-degree node within distance r (Model-max), and 2). the new nodes are connected to the lowest-degree node within distance r (Model-min). We find that Model-min and Model-max replicate the patterns observed in the mobility network and in the attention network, respectively.
We suggest that this is because both the browsing behavior and the physical movement of mobile phone users are governed by two underlying driven forces, the exploration of new sites and the preferential return to favorite sites.
The exploration tendency is much stronger in the virtual world than that in the physical space. With a high probability users hop between popular and non-popular websites, connecting high degree nodes with low degree ones.
In contrast, in the physical world users preferentially return to living and/or working places. To get back to the preferred place, one has to go through a sequence of consecutive base stations in the local clusters. Thus, whenever users are visiting a new base station that is seldom visited (a low degree node) during a long trip across the city, it is very likely that they will meet another less frequently visited base station.