Modeling the St. Louis Arch - Part 2
Posted on 2019-11-23 in how-to • 9 min read
Continuation¶
When we left off last time, we had gotten to a 2D representation of the centroid of the arch cross-sections. However, the centroid is used as “construction” geometry in this model - it doesn’t represent something physical. Rather, it is used for constructing, or laying out, portions of the model that will actually be built. Yes, this is a little confusing because we are talking about construction in terms of concrete and steel, but construction geometry as a term used in this context is common in 2D CADD and other modeling packages.
To get started, we’ll import the same packages, use the same constants, and utilize the same functions previously developed.
Note: Typically at this point in a data science project, this would be a great time to start factoring notebook cells out into scripts and modules. However, I’m going to make the call that I can live with some technical debt for the time being.
import math
import matplotlib.pyplot as plt
plt.rcParams['figure.figsize'] = [11, 8.5]
from mpl_toolkits.mplot3d import Axes3D
import numpy as np
import mpmath
import pandas as pd
fc = 625.0925 # Maximium height of centroid
Qb = 1262.6651 # Maximum cross-sectional area at arch base
Qt = 125.1406 # Minimum cross-sectional area at arch top
L = 299.2239 # Half of centroid at arch base
A = fc / ((Qb / Qt) - 1)
C = math.acosh(Qb / Qt)
def El(Y):
return fc - Y
def X(Y):
return (L / C) * (math.acosh(1 + Y/A))
Intrados and Extrados¶
The NPS web site has served us well thus far. However, we could certainly use some more detail on a methodology for calculating the inner and outer edges of the arch to fully develop a 2D elevation view. The detailed methodology for developing the arch geometry was laid out in the blueprints.
A Word About Symmetry¶
The equations developed for cross-sectional charactistics of the arch were developed for the right-hand side of the arch. Speaking from experience, bad things will happen if you try to apply the same relationships to the left (x < 0) side of the arch. Therefore, we’ll just calculate the right leg, then mirror about the X axis to get the other half of our model.
With that important caveat in mind, let’s build pythonic versions of these equations:
def Y(X):
ret = A * (math.cosh((C * X) / L) - 1)
if ret == 0:
ret += 1e-16
return ret
def slope(Y):
return (L / C) * (1 / math.sqrt((2 * A * Y + Y**2)))
def alpha(slope):
return math.atan(slope)
def Q(Y):
"""cross-sectional area at any Y"""
return ((Qb - Qt) / fc) * Y + Qt
def H(Q):
"""height of the equilateral triangular section"""
return math.sqrt(Q * mpmath.cot(mpmath.radians(30)))
def W2(H):
"""side of the equilateral triangular section"""
return 2 * H * math.tan(math.radians(30))
# test some known values
print(f"Value of Y at X=0 is {Y(0):.4f}")
print(f"Elevation at X=0 is {El(Y(0)):.4f}")
print(f"Maximum cross-sectional area is {Q(0)}, or Qb = {Qt}")
print(f"Maximum cross-sectional area is {Q(Y(L))}, or Qt = {Qb}")
print(f"Height of triangle at the top is {H(Qt):.4f}")
print(f"Length of triangle side at the top is {W2(H(Qt)):.4f}")
print(f"Height of triangle at the bottom is {H(Qb):.4f}")
print(f"Length of triangle side at the bottom is {W2(H(Qb)):.4f}")
print(f"Slope angle (alpha) near top of arch is {math.degrees(alpha(slope(0.000001))):.4f}")
Great! Everything is checking out so far.
Side note: ideally, we would be using a test suite such as pytest in lieu of print-mode debugging. However, in a notebook such as this the focus is more on telling a story than creating a piece of production software.
Data Structure + Function = Class¶
Initially, my next step was to start defining individual functions to calculate points on the interior and exterior of the arch at any given location along the centroid. But, that would have duplicated a lot of the calculation logic. So, let’s think about what we actually want to do: given any x coordinate on the centroid, calculate the y coordinate of the centroid, along with the coordinates of the intrados and extrados. This is the perfect time to write our first class. This will allow us to encapsulate the calculations and results into a single, tidy spot.
While we’re at it, let’s go ahead and make the leap to three dimensions. One wrinkle here is that our coordinate system defined thus far is consistent with a “Y-up” convention that is common in graphics. This can be confusing for those of us that are used to “Z-up” in the CADD world. For the arch, the centroid and intrados will be on the Z-plane (Z=0). So that means we only need to worry about calculating Z coordinates for the extrados. Since the arch doesn’t have any torsion (twist) by design, we actually only need to calculate a single Z offset, then assign the positive value to one extrados point and the negated value to the other extrados point. This Z offset will be one-half of W, the length of each side of the triangular cross-section.
class NormalSection():
"""calculates the triangular cross-section of the arch at any provided
value of x (offset from axis of symmetry)"""
def __init__(self, x):
max_x = 299.5458
if x < 0 or x > max_x:
raise ValueError(f"Invalid value for x (must be between 0 and {max_x})")
if x == 0:
x += 1e-6
y = Y(x)
a = math.atan(slope(y))
h = H(Q(y))
h1 = 2 * h / 3
h2 = h1 / 2
z = 0.5 * W2(h)
# intrados
x3 = x - h1 * math.cos(a)
y3 = y + h1 * math.sin(a)
# extrados
x1 = x + h2 * math.cos(a)
y1 = y - h2 * math.sin(a)
# counterclockwise winding
self._pt1 = [x3, El(y3), 0]
self._pt2 = [x1, El(y1), -z]
self._pt3 = [x1, El(y1), z]
@property
def vertices(self):
return np.array([self._pt1, self._pt2, self._pt3], dtype='float64')
An easy way to test this class is to confirm that the extrados point at the top of the arch is equal to the total height of the arch (630 feet). The distance between legs of the arch at the base should be 630 feet as well.
print(f"Maximum height of arch is {NormalSection(0).vertices[1][1]:.4f}.")
print(f"Total width of arch at base is {(2 * NormalSection(299.5458).vertices[1][0]):.4f}.")
Blueprint Stations¶
Now that we have a class defined to provide triangle vertices at any value of x on the centroid construction geometry, we need to know what values of x should actually be evaluated. The arch blueprints included a list of 72 stations with these values.
We can use pandas to parse this file so that we are creating our virtual model of the arch geometry exactly as intended per the blueprints.
Note: I renamed my local copy to have a *.txt extension
so that it didn’t wreak havoc with the static site generator tooling.
The read_html method returns a list of dataframes corresponding to
all of the tables in an html document.
From inspection of the input data, we know that we only want to keep the first one.
df = pd.read_html('../static/centroidxdat.txt', header=4, index_col=0)[0]
df.drop(df.columns[1:], axis="columns", inplace=True)
df.rename({df.columns[0]: "centroid_x"}, axis=1, inplace=True)
stations = df.loc[:, "centroid_x"]
y = np.array([El(Y(_)) for _ in stations])
sections = [NormalSection(_) for _ in stations]
x3 = np.array([_[0] for _ in [p.vertices[0] for p in sections]])
y3 = np.array([_[1] for _ in [p.vertices[0] for p in sections]])
x1w = np.array([_[0] for _ in [p.vertices[1] for p in sections]])
y1w = np.array([_[1] for _ in [p.vertices[1] for p in sections]])
x1e = np.array([_[0] for _ in [p.vertices[2] for p in sections]])
y1e = np.array([_[1] for _ in [p.vertices[2] for p in sections]])
triangles2d = np.float64([[[s.vertices[0, 0], s.vertices[0, 1]], [s.vertices[1, 0], s.vertices[1, 1]]] for s in sections])
x2d = triangles2d[:, :, 0]
y2d = triangles2d[:, :, 1]
axes = plt.gca()
axes.set_aspect('equal')
axes.set_ylim(0,650)
plt.plot(stations, y, label='centroid', color='green')
plt.plot(-stations, y, label='centroid - north', color='green')
plt.plot(x3, y3, label='intrados - south', color='red')
plt.plot(-x3, y3, label='intrados - north', color='red')
plt.plot(x1w, y1w, label='extrados - south', color='blue')
plt.plot(-x1w, y1w, label='extrados - north', color='blue')
for xs, ys in zip(x2d, y2d):
plt.plot(xs, ys, color='gray')
plt.plot(-xs, ys, color='gray')
plt.title("Fig. 1 - Elevation View of Extrados, Centroid, Intrados, and Sections")
plt.show()
Confirm Results¶
Ok, so that looks good.
However, I said earlier that we were going to go ahead and make the jump to 3D.
Our NormalSection class calculates 3D points for all of the triangular cross sections.
So, with a few more lines of code, we should be able to make a similar plot,
but with the inclusion of the Z-axis,
which is the direction coming in and out of the screen in Figure 1.
z3 = np.array([_[2] for _ in [p.vertices[0] for p in sections]])
z1w = np.array([_[2] for _ in [p.vertices[1] for p in sections]])
z1e = np.array([_[2] for _ in [p.vertices[2] for p in sections]])
# add the first point again at the end to make a closed shape
triangles3d = np.float64([[[s.vertices[0, 0], s.vertices[0, 1], s.vertices[0, 2]],
[s.vertices[1, 0], s.vertices[1, 1], s.vertices[1, 2]],
[s.vertices[2, 0], s.vertices[2, 1], s.vertices[2, 2]],
[s.vertices[0, 0], s.vertices[0, 1], s.vertices[0, 2]]] for s in sections])
x3d = triangles3d[:, :, 0]
y3d = triangles3d[:, :, 1]
z3d = triangles3d[:, :, 2]
fig= plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot(xs=stations, ys=np.zeros_like(stations), zs=y, color='orange')
ax.plot(xs=-stations, ys=np.zeros_like(-stations), zs=y, color='orange')
ax.plot(xs=x3, ys=z3, zs=y3, color='red')
ax.plot(xs=-x3, ys=z3, zs=y3, color='red')
ax.plot(xs=x1w, ys=z1w, zs=y1w, color='blue')
ax.plot(xs=-x1w, ys=z1w, zs=y1w, color='blue')
ax.plot(xs=x1e, ys=z1e, zs=y1e, color='green')
ax.plot(xs=-x1e, ys=z1e, zs=y1e, color='green')
for x3d, y3d, z3d in zip(x3d, y3d, z3d):
plt.plot(xs=x3d, ys=z3d, zs=y3d, color='gray')
plt.plot(xs=-x3d, ys=z3d, zs=y3d, color='gray')
ax.set_xlim(-350, 350)
ax.set_ylim(-350, 350)
ax.set_zlim(0, 700)
plt.title("Fig. 2 - 3D View of Extrados, Centroid, Intrados, and Sections")
Serialize output¶
Now that we’ve calculated all of the geometry, let’s serialize it (write coordinates to a text file) so that we can potentially re-generate this model in another platform such as Blender.
We’re going to reverse the order of the triangles (arch cross-sections) so that it starts at the bottom and progresses up to the top. That way we can build our model in the same order as the way that actual construction was carried out.
import json
with open("triangles.json", "w") as f:
json.dump(triangles3d.tolist()[::-1], f, indent=2)
Conclusion¶
Given the density of the sections, it’s hard to pick out the centroid in the 3d plot. This makes sense, given that the centroid is merely ‘construction’ geometry and does not represent anything physically present in the real world.
Matplotlib does not provide any interaction out of the box. This means we cannot zoom, pan, or otherwise adjust our view of the arch model. However, now that we have our geometry, we can use any number of external libraries for rendering the geometry dynamically in the browser. Similarly, we could write the geometry as a DXF polyface mesh and load it into a CAD package such as MicroStation. These next steps will be left for a (potential) future blog post.
Why all this code? Why not just model the arch directly from a spreadsheet or other data source?¶
There are three main benefits of utilizing computational geometry for 3D modeling such as I have done with the Gateway Arch:
1. Parametric Values for Maximum Flexibility and Iteration Speed¶
The entire process, from start to finish, is based upon parameters that define the overall arch geometry. If this were a proposed project still in the design phase, we could easily do multiple variations as part of the ‘what if?’ iterative nature of design. For example, we could change one or more of the constants to scale the arch to a greater height, or adjust another to vary the arch from the 1:1 aspect ratio.
2. Change management with source control tools such as git.¶
As any architect or engineer will tell you, change is a constant during the design process. Ideas will be discarded one day only to be resurrected days, weeks, or months later. Because we are developing all of the geometry with python code, the parametric relationships are all captured in plain-text files versus proprietary binary formats typical in CADD and BIM software. Utilizing a tool such as git means that every iteration can be checkpointed with a commit hash, timestamp, and optional tag. Reverting to a previous version of the design (or even carrying multiple variations forward at the same time) is as easy as checking out a specific commit or branch in the revision history.
3. Reducing human error¶
When geometry is being created or entered by a human operator, there is always the possibility of ‘fat-fingering’ the data via typos or other transpositional bugaboos. By allowing a computer to generate the coordinate geometry algorithmically, we eliminate this particular source of error. Yes, there is still the potential for bugs or other errors in the code, but careful checking along the way either manually (as I have done in this example) or via an automated test suite (best practice when using code in production), those potential bugs can be fleshed out early in the process.