The field of urban studies is a continuous pursuit for regularities that expand our capacity to describe and understand cities. From its origins, the field has been intimately related to ideas and methodologies in the field of statistical physics, and an overarching summary of the path and main ideas of the application of statistical physics to urban environments is presented in [1].

We cannot claim that we understand how cities evolve as long as we do not have an exact set of equations that relate every variable to each other allowing us to understand the effects of population growth. This is paramount for a large number of fields, including research, urbanism and political/economical science. The aim of this work is to create an analytical framework for the analysis of the most important variables in a city from a geometrical standpoint.

This paper presents a theory of the physical aspect of the city using a fractal framework. Cities are defined by their occupation of space, and I show over the next sections that cities increase their fractal dimension as they grow in population, which is a fundamental property, and which was already noted in [2]. In fact, the study of cities as fractals has a long tradition in the scientific literature [2–8]. This view is fundamental to understand cities, since it commands the occupation of space for a given city, and therefore it is the only valid way to extract its geometric framework.

The work presented in this paper studies the population, the road network, the occupation of space, the fractal dimension of the city, the average heights, the maximum heights and the interactions or GDP of a city. I show how all these variables are related between each other and how they were derived. This work is not the first to attempt to produce a set of equations that describe the main variables of the city. Some previous work include [9, 10] which show how congestion influence the growth of cities, the aforementioned [2] that presents a theory of growth of cities based on scaling theory or several theories on the growth of cities [11–13].

This work delves as well into how this fractal framework affects and influence the scaling of a number of variables. Scaling theory [14–17] in urban science [18–23] study how allometric relations appear between variables such as the length of roads, the number of gas stations, the cost of maintaining a city’s infrastructure, its GDP and many other quantities to the size of the city in terms of population. In previous work [24], it was shown that those scaling exponents could be derived from a simplified and approximated version of the equations that governed those same variables. The current work presents an improvement on the calculation of the scaling exponents going beyond what was presented in [24], since we now account with better and more precise equations to describe length of roads, population and GDP.

The work also shows an application of the set of formulas to develop an approximation of the internal distribution of a city, derived from its growth. In order to do so, it assumes that the same formulas that describe a city in its current state have to remain meaningful to describe a city at any specific instant of its history, meaning that this growth is ergodic.

Cities begin to form with the construction of a single house, slowly other houses join, percolating space, and soon the city’s fractal dimension starts increasing. New occupied sites are incorporated with a certain probability of occupation over the territory that surrounds the city while the existing urban fabric gets densified. This probability of occupation tends to a constant value because of a self-optimization pattern. Going any lower would break the city apart into different clusters and expand the city over a long area, increasing travel times and decreasing economies of scale. Going higher would increase traffic and other problems derived from density, leading to an optimal solution in which the probability of occupation is as low as possible while still keeping the city as one single cluster. This probability is very close to the critical probability of a percolation over a squared lattice in two dimensions because the topology of cities is in average similar to that lattice.

In its first stages the city densifies its road network, subdividing the occupied sites which increases the density of occupation. At some point in its growth, the density of people that can live in a planar city saturates and in order to keep growing, the city needs to extend into the third dimension increasing its height, eventually pushing the fractal dimension of the population above 2. From its initial state, the density of population keeps on ever increasing as the city first densifies and later grows further into the third dimension.

The current section constructs the main framework of this work, presenting the derivation of its equations for several geometrical variables in a city. In order to simplify the equations and reasoning, I will use along this section an idealised system, in which I will avoid talking about multipliers or characteristic scales that need to be fitted in order to obtain realistic values, I will show how to calculate those multipliers in the last section of the paper, where I will adapt the equations to work with a real system and give the value for the constants in the specific case of the United Kingdom which will serve as an example throughout the paper.

As a city grows, the pre-existing city does not disappear, meaning that it must maintain a minimum of the current probability of occupation of sites where buildings are constructed. Considering a squared city with a linear dimension L, with a planar fractal dimension d and a certain number of occupied sites n = L ^d (shown in Figure 1B) coming from a probability of occupation ρ over an area a = L ². Then we have ρa = n which in turn means that L ^d−2 = ρ so: d = 2 + ln ρ ln L (1)

This equation means that the fractal dimension of a city needs to increase as its linear dimension grows in order to avoid having its preexisting city disappear, this was already noted in [2]. Otherwise, if it were to remain constant (or decrease) we see that the only solution would be to decrease its probability of occupation, meaning that in order to occupy sites in its outskirts the city would need to vacate sites in its preexisting city. This explains the behaviour observed in real systems (Figure 1A) which shows a normal distribution of the error between the predicted and real value with mean −0.0142 and standard deviation of 0.0677.

Cities grow vertically above its two-dimensional footprint. In fact, population becomes a fractal volume, that starts below dimension 2 but as cities become larger it surpasses it. Furthermore, the population has always a larger fractal dimension than the footprint of a city. As it was shown on [24] the fractal dimension of the population (d _p ) can be obtained by adding a fractal vertical component η to the planar fractal dimension of sites d, that is d _p = d + η and therefore the population (shown in Figure 1C) can be expressed as p = L d + η (2)

Throughout its growth, a city starts densifying its street network increasing the quantity of people that can live within it, and at some threshold x _c the density of population per site in a planar city saturates and cannot longer continue growing through this process. In order to keep on growing above that critical threshold x _c it needs to start increasing its height. Therefore, the growth of people per site p n = L η is absorbed either by the increase in average length of roads ⟨ℓ⟩ in a site of the space of the city or by the average number of vertical levels (floors) of the city ⟨h⟩. This means that p n 〈 ℓ 〉 〈 h 〉 = k (3) is a constant.

Numerically, at x _c , the density of roads reaches its maximum value of 1 (understood as the probability of finding a road segment in a site) ⟨ℓ⟩ _c = 1 and the average number of levels of the city is also 1, ⟨h⟩ _c = 1, since it was 1 from the beginning of the growth of the city and it only starts increasing right after x _c . So we have that: 〈 ℓ 〉 c 〈 h 〉 c = 1 = p c k n c = 1 k L c η (4)

We also have that since ρa = n, p a 〈 ℓ 〉 〈 h 〉 = ρ k (5)

This last equation means that at x _c the population is proportional to the area since ⟨ℓ⟩ _c = 1 and ⟨h⟩ _c = 1 and k and ρ are constants. Therefore, the fractal dimension of the population at x _c is the same as the fractal dimension of the area ( d p c = d c + η = 2 ) . In other terms, at the point in which the population starts growing into the third dimension is the threshold in which its fractal dimension grows above 2, as expected. We can see that since at x _c , d _c = 2 − η and using Eq. 1 we have that η is η = − log ρ log L c (6) and since p _c = L ² = a _c and a/n = ρ ⁻¹ then L c η = p c / n c = ρ − 1 , therefore, we can see that our constant in Eq. 4 is k = ρ ⁻¹ which means that, p n 〈 ℓ 〉 〈 h 〉 = ρ − 1 (7) for any city size. This is a constant for any city, the limit of population density in each vertical floor belonging to a site per meter of road (how dense is the network, how subdivided is the system), the system cannot hold more people than this value per site, per level. Furthermore, Eq. 5 becomes p a 〈 ℓ 〉 〈 h 〉 = 1 (8) which means that at the critical threshold x _c , when saturation is reached and both the density of roads and the number of levels is 1, the population equals the area. Moreover, and as an indirect consequence of Eq. 6 we also have that the threshold x _c is reached when the linear dimension of the city is L _c = ρ ^−1/η .

A city cannot grow without limit, the equations portrayed in this work show that its density would go to impossible amounts, the heights of its buildings would reach levels that are physically unattainable and a large number of other issues such as congestion and competition for space would arise. As we see in Eq. 7, ρ ⁻¹ is the density limit for the night time population, how many people can live in each site per level, and this is a hard limit, no city outgrows this. As cities grow further than this threshold (x _c ), the day-time population will spread over its area, people will walk down the parks, the plazas, the avenues and will of course be present in buildings. Since the area of the city cannot contain several levels, it means that at some point in its growth, the population spread over its area will also reach this same limit p a = ρ − 1 and the city will be completely collapsed. This is expectedly a hard limit for growth, and cities will not surpass it. In fact, using our equations and the multiplying factors measured for the United Kingdom, this only happens when the population of a city reaches 36.1 million people, and the largest city on Earth has a very similar level of population, of course, this is a numerical result that is highly dependent on the approximated value of η and ρ obtained through a genetic algorithm as explained in the next section, and very slight modifications to those values change greatly this specific value. We call this maximum threshold x _m .

If at x _m we have that ρp _m = a _m and for all cities ρa = n, then p m n m = L m η = ρ − 2 which means that L _m = ρ ^−2/η and using Eq. 1, d _m = 2 − η/2.

Regarding ⟨ℓ⟩ and ⟨h⟩ given that we know that both are complementary, since at x _c both are 1, and one cannot exist in the numeric range of the other (one has to be less than 1 and the other more than 1) then we have that ⟨ h ⟩ ⟨ ℓ ⟩ = p a = ρ p n = ρ L η = ρ 1 + η d − 2 and given that ⟨ℓ⟩ is a probability and cannot go above the value 1: ⟨ ℓ ⟩ = min ρ 1 + η d − 2 , 1 = min L d + η − 2 , 1 (9)

which is shown in Figure 2D, and since ⟨ h ⟩ = ρ p n ⟨ ℓ ⟩ then: ⟨ h ⟩ = max 1 , ρ L η = max 1 , ρ 1 + η d − 2 = max 1 , L d + η − 2 (10) portrayed in Figure 2B.

We can also obtain from Eq. 9 the total length of roads (shown in Figure 1D) which is: ℓ = ⟨ ℓ ⟩ n = min 1 , ρ 1 + η d − 2 L d = min 1 , L d + η − 2 L d (11)

In order to calculate the maximum number of levels of the city (⌈h⌉) we use an approximation from [24] where we obtained that d p = d + η ∼ d + ln ⟨ h ⟩ ln ⌈ h ⌉ using box counting analytically. Since this was only an approximation and meant to work above the critical threshold, we need to consider ⟨h⟩ _x = ρL ^η (without limiting its lower bound to 1). Doing so we have that η ∼ ln ⟨ h ⟩ ln ⌈ h ⌉ (Figure 2C) which means that ⌈ h ⌉ ∼ exp ( ln ( ρ L n ) η ) ∼ k h L . Because of the approximated nature of the equation we can not obtain directly the value of k _h and to calculate the true value of this multiplying constant we have to see that at the minimum possible linear dimension, when the city is composed of a single house, the average height and the maximum height should have the same value. This happens when d = 0, and since L = ρ ^1/(d−2) then L _d=0 = ρ ^−1/2, the average number of levels is ⟨ h ⟩ d = 0 = ρ L d = 0 η = ρ 1 − η / 2 and equating it to ⌈h⌉ _d=0 = k _h L _d=0 = k _h ρ ^−1/2 we obtain that k _h = ρ ^1.5−η/2 therefore: ⌈ h ⌉ = ρ 1.5 − η / 2 L = L d 3 − η + η − 1 2 (12) which is shown in Figure 2A.

We can also obtain equations that describe the average number of population in a site projected to the floor (collapsing all levels) ⟨p _p ⟩ (Figure 2E), the average number of people per meter of road ⟨p _ℓ ⟩ and the average number of people per site and per level ⟨p _h ⟩ ⟨ p p ⟩ = p n = L η (13) ⟨ p l ⟩ = p l = max L η , L d − 2 + 2 η (14) 〈 p h 〉 = p n 〈 h 〉 = min L η , L d − 2 (15)

Over the next sections I will show how to adapt this framework to real data and how to obtain the value of η and ρ to be able to get the final values of our exponents.

Originally a theory derived in the field of biology, scaling theory studies the allometric scaling of variables in a city as it they relate to its population growth. Some of those variables scale sub-linearly with the size of the city, meaning that the larger a city gets the slower that variable grows, this is the case for variables where economies of scale arise, such as the length of roads needed to cover the city, the number of gas stations, etc. Other variables grow linearly, because they correspond to some fixed value per person, such as the amount of water consumed. Finally, some other variables grow super-linearly, meaning that they grow faster than the population, usually these arise through feedback effects and include elements like traffic congestion, criminality, interactions or the GDP of a city. In [24] we showed that in fact, this relation to size was due to the fractal nature of cities, and calculated the expected exponents from the fractal scaling of the population and road network.

In that previous work [24] we reasoned that since the length of roads was proportional to L ^d and the population was proportional to L d p , the scaling exponent should be equal to γ = d d p = d d + η . We now have an improved formulation that describe both quantities and no longer need to make any approximations we can just solve the scaling equation ℓ = p ^γ and calculate the exponent γ = ln ⁡ ℓ ln ⁡ p . This gives us that: γ = min d d + η , 1 + d − 2 d + η (16) this means that the previous reasoning still stands, but only for the largest cities (those above x _c ) but the small cities are better represented by a different γ. The resulting length of roads fits the data to a very high degree of precision (normal distribution of the differences between the logged real and predicted values with μ = 0.01466 and σ = 0.1618) as shown in Figure 3A where Figure 3C shows the value of γ.

In that same paper, we also reasoned that interactions occur when people go to the street and that therefore it should be proportional to the square of the quantity of people in the ground level, multiplied by the number of locations in which that were possible. In that work, the equations were approximated and we used ℓ ∼ n which gave us that i = (p/n) (p/n − 1)n, where i represents the total possible interactions, but since in this work we are distinguishing between the two values (ℓ and n), it is more precise to say that the people in the street interacts, and the number of possible locations is the length of the street network. therefore: i = p ℓ p ℓ − 1 ℓ = p p ℓ − 1 = max L d + 2 η , L 2 + η − L d + η (17)

In order to obtain an approximation of the super-linear exponent of interaction i ∼ p γ sup , we need to drop the exact equation and approximate i ∼ p ℓ 2 ℓ = p 2 ℓ = p 2 p γ = p 2 − γ . Therefore the approximate value of the super-linear exponent is: γ sup ∼ 2 − γ = max d + 2 η d + η , 2 + η d + η (18)

Similarly to what was done in [24], we assume that the GDP of a city is a direct consequence of the interactions between individuals, and use that quantity to showcase the validity of the formulation as shown in Figure 3B while the super-linear exponent is shown in Figure 3C.

The current framework represents an idealised system, it is unitless because everything is divided by an implicit characteristic scale that we will make explicit in this section and there are not any multiplying constants for the sake of simplicity. This section completes the framework, by including those factors and thus creating the final set of equations for the system.

From this point on forward we will use the subindex r to refer to real variables as measured from the data.

To determine the side of our real square (in meters) we use the area of the city. L r = a r 1 2 (19)

The characteristic scale for the length of the side of our squared area is called L ₀ and it is measured in meters. L = L r L 0 (20) where for the United Kingdom L ₀ = 538.924 m. This value was calculated through measuring the fractal dimension for all cities and their areas. An approximation can be obtained through performing those measurements for the largest city (max (d _r ), max (L _r )) and calculating our theoretical max(L) = exp (ln(ρ)/(max (d _r ) − 2)), to find L ₀ = max (L _r )/max(L). Of course, for this we need to determine ρ, this can be done either through directly measuring occupied space (buildings and roads) against open spaces in the city (parks, and plazas) or assuming our theoretical value for ρ = 0.5991 taken from the next section. As an example, London has an approximated 40% surface occupied by parks, which means that its ρ _r = 0.60. I use this value for L ₀ as a starting guess and perform a least square estimate of L ₀ using the measured fractal dimension and area for all cities (I assumed the theoretical ρ to be valid). Notice that we cannot use L _r and d to directly calculate ρ because. ρ r = ( L r / L 0 ) d − 2 ≠ L r d − 2

For completeness, we will show how to obtain the area as a function of the side of the square. a r = a 0 L r L 0 2 (21) where a 0 = L 0 2 m 2 d = 2 + ln ρ ln L r / L 0 (22) p r = p 0 L r L 0 d + η (23) where for the United Kingdom p ₀ = 1,207 people, which was measured by adjusting the theoretical population to the real data measured until their differences were minimised. n r = n 0 L r L 0 d (24) where for the United Kingdom n ₀ = 7.85 sites, measured against the real data. ℓ r = ℓ 0 min L r L 0 d , L r L 0 2 d + η − 2 (25) where for the United Kingdom ℓ ₀ = 7,700 m, measured against the real data. ⟨ ℓ ⟩ r = ℓ 0 n 0 min 1 , L r L 0 d + η − 2 (26) ⟨ h ⟩ r = h 0 max 1 , L r L 0 d + η − 2 (27) where for the United Kingdom h ₀ = 2 levels, measured against the real data. h m r = h 0 L r L 0 d 3 − η + η − 1 2 (28)

The constant that limits growth becomes: p r ℓ r ⟨ h ⟩ r = ρ − 1 p 0 ℓ 0 h 0 people m level (29)

and the densities: ⟨ p p r ⟩ = p 0 n 0 L r L 0 η (30) ⟨ p l r ⟩ = p 0 l 0 max L r L 0 η , L r L 0 d − 2 + 2 η (31) ⟨ p h r ⟩ = p 0 n 0 h 0 min L r L 0 η , L r L 0 d − 2 (32)

Regarding the scaling equations, we have that the real length of the road network as a function of the real population becomes: ℓ r = ℓ 0 p r p 0 γ (33) Notice that, the typical equation of a scaling ℓ r ∝ p r γ is in fact ℓ r = ( ℓ 0 / p 0 γ ) p r γ , and since γ changes with the population size, the multiplying factor is not constant, as the equation ℓ r ∝ p r γ assumes. The fact that both the exponent and the multiplying factor vary with size explains a lot of the problematic that exists around measuring precisely the scaling exponent, although much of the variability becomes negligible if we only use cities above x _c .

For the GDP we have: g r ∼ g 0 p r p 0 2 − γ (34) where for the United Kingdom g ₀ = 2.3 ⋅ 10⁷ euro, measured against the real data. This is an approximation, and from my perspective it is preferred to use the actual equation instead of an approximated scaling law, whenever possible, even though both equations look indistinguishable when presented against each other or the data. The equation for the interaction of population is: i r = i 0 max L r L 0 d + 2 η , L r L 0 2 + η − L r L 0 d + η (35) with a value i ₀ unknown, since there is no data to measure it. This factor i ₀ represents the probability that a potential interaction becomes a real one. Furthermore, if the assumption between proportionality of interactions and GDP stands: g r = g 0 max L r L 0 d + 2 η , L r L 0 2 + η − L r L 0 d + η (36)

One interesting side effect of this, is that given that L ₀, n ₀, p ₀, ℓ ₀ and h ₀ or even g ₀ are pure constants for a system of cities (the variability is absorbed through the rest of the equation), they are much better descriptors of a system, and when calculated for other systems, they will allow us to make comparisons with less noise between different countries.

To render this analytical approach useful we need to be able to obtain the values of our two constants η and ρ. My approach was to use a genetic algorithm, whose inputs were the area (a _r , from where we obtain L _r ), fractal dimension (d _r ) and population (p _r ) for each city and the parameters to be optimised are η _p and ρ _p . Using this, I apply a two steps approximation.

In the first step, in order to obtain the heuristic value for each individual, I calculate L ₀ using the parameter ρ _p given by the algorithm, and then after calculating the real density ρ r = ( L r / L 0 ) d r − 2 for each city, calculate an average density that mixes the parameter and the measured density, ρ ̂ = ( ρ p + ρ r ) / 2 . I then calculate a theoretical fractal dimension d = 2 + ln ( ρ ̂ ) / ln ( L r / L 0 ) and obtain the theoretical population for each city using it p = ( L r / L 0 ) d + η p , where η _p is the second parameter to be fitted by the algorithm. Then obtain the linear fit between the theoretical population p and the real one p _r = ap + b, discard b and return as the final heuristic the L1 norm between the logarithms of p _r and ap. In this phase we obtain the value of η and an approximation of the value of ρ.

In the second step, we fix η and only optimize ρ, allowing the search only in the neighborhood of the approximated value we obtained in the first step. The only difference between the two, is that we no longer use the real ρ _r and instead use directly the parameter ρ _p at every step (to calculate L ₀ and d = 2 + ln (ρ _p )/ln (L _r /L ₀)). From this second step we obtain our constant ρ.

The values obtained were: η = 0.2092 (37) ρ = 0.5991 (38) were less significant digits were discarded. These values mean that d _m = 2 − η/2 = 1.8954 and d _c = 2 − η = 1.7908. The measurements of η and ρ were obtained from approximated processes and these measurements could be improved in the future.

This is surprisingly close to the values for a site percolation in a 2d-lattice given in the literature, were η = 0.2083 is the exponent for the function that controls the probability of two sites belonging to the same cluster as a function of distance, p _c = 0.5927 is the critical probability, d _f = 1.8958 is the fractal dimension of the percolating cluster. Given that percolation has been tied in the literature [25] to the formation of cities, I expect that there exists a logical link between the two but the reasoning behind this numerical coincidence falls outside the scope of this paper and is left for future work.

Scaling studies have shown that different city systems across the world have very similar scaling exponents [18]. Following our derivation we see from Eq. 16, 18 that the scaling exponent depends on d and η. Since d is a function of the linear size of the city and ρ (Eq. 1) the scaling exponent is a function of ρ and η. If the scaling exponents are truly universal it would then mean that in fact ρ and η are universal and therefore these values should remain stable for different systems.

In the following section we show a possible application of the framework contained in this paper, in order to demonstrate its expressiveness.

We can apply the same reasoning presented above to obtain the internal distribution of a single city, since at each stage of its growth, the city has to follow the equations presented for fractal dimension, population, number of sites, and length of roads if we consider urban growth to be an ergodic process.

Upon growth, the city increases from a current linear length L to L + dL. In this change of linear size, it modifies its fractal dimension from d _L to d _L+dL , and its population change is d p = ( L + d L ) d L + d L + η − L d L + η .

This population change will be partitioned between the stripe of land added to the city and the existing urban tissue. I assumed a simple formula for this, that uses a weight to balance the two, w. We then consider a value δ that is the density of population added at each step, which multiplied by the respective areas gives us the increase of population in the new area and the preexisting one. The basic formula for the population at a certain stage of its growth is then: p L + d L = p L + δ w L 2 + 1 − w L + d L 2 − L 2 (39)

We can calculate p _L+dL and p _L using Eq. 1, 2 and their difference is the increment of population dp. So we can express δ as δ = d p L + d L 2 1 − w + L 2 2 w − 1 (40)

where w is adjusted to fit the real distributions, in the modeling process w has been made dependent on the step size, so variations on the step size would not influence the final distribution, the adjusted value was w = d L 25 , of course this can only be valid as long as w < 1, our step chosen for the model was dL = 1, the value used in the modeling was w = 1 25 .

Of course, as we add new population, each city stripe must remain under the maximum possible population. This maximum possible population can be calculated from Eq. 8, where max(p) = max (a⟨ℓ⟩⟨h⟩) = a max ⟨ℓ⟩ max ⟨h⟩ and since max ⟨ℓ⟩ = 1 then max(p) = a⌈h⌉. So each stripe must remain below its area multiplied by the maximum height of the city for the current linear dimension.

When deciding where to locate in the city, a new inhabitant only cares on the distance to the center, in order to simulate this extent when distributing the population (δwL ²) over the pre-existing city, we weight each strip by how many more people fit in it, divided by its perimeter, max ( p i ) − p i 4 L i and use this factor to distribute the population on the existing city.

We can repeat the same steps for the number of sites and the length of roads, obtaining the most important variables. For number of sites, we choose a maximum possible density of 0.9 (being 1 complete occupation), this value was obtained from the data observed, while length of roads is limited by the number of sites. From it we can calculate the heights of buildings expected and the density of sites per area or of people per site.

The height of buildings for the real data is a direct measurement taken from the LIDAR available at the Copernicus site [27] and no transformation was applied other than dividing it by 3 m, which is taken as an average floor height, this is shown in Figure 4E. In order to calculate the density of a site, we calculate how many occupied sites (there exists population in that element of the grid using data from the GHS [26]) are in the surrounding area of each site (with a radius of 6.250 km) and divided it by the maximum possible number of sites in that circle, as shown in Figure 4F. The last comparison (population per site) is more complicated, and we need to think how this data was created. The population data is obtained from the Global Human Settlement layer [26], this data has been produced by taking the population in censal sections, determining the building footprints from satellite data and interpolating the population with the perceived density of buildings, also, most probably, since we do not see any clear cuts from the censal sections, a spatial interpolation averaging large discontinuities was performed. Both interpolations (and even the data aggregated to a censal section) reduce the peaks of population, softening the overall distribution. Therefore, in order to create a fair comparison, we performed similar steps to our results. In Figure 4C both the real distribution obtained (dotted points) and the distribution obtained after a process of clustering and interpolation is shown.

The correspondence of the distributions obtained using the model with the real ones is fairly strong, indicating that this process could be a valid model for the internal growth of a city. However, and as we can notice in Figure 4B we can see that at the outskirts of the city there is a strange behaviour, where the height of buildings start growing again instead of decreasing, which means that there is still room for improvement. This problematic is created because the number of sites decreases faster than the population for that range.

The analytical derivations that give rise to the equations portrayed in this work, makes them exact functional forms of many aspects of the city’s physical environment. This is of extreme importance, since every derivation made from them, every operation will still represent what they are meant to convey. As it is often said, we stand on the shoulder of giants, and approximated equations of similar quantities have been portrayed before in the literature, and while these brought light to a lot of issues they are of limited applicability, because of their approximated nature.

I believe that following this text, new ideas will become easier to test and derive, aiding the process of solving the puzzle of cities.

This work portrays the equations for fractal dimension, population, area, length of roads, different densities of population, average and maximum heights (levels) for a city, and interactions (or GDP). Moreover, it shows how using this framework we can study the internal distributions of those same variables within the city.

The article uses population data from the Global Human Settlement Layer (GHS) [26], height data from the Copernicus satellite LIDAR data [27], height and road data from OpenStreetMap [28] and GDP data from Eurostat [29].