PSI Wave => Buildings

by Buddy McCloskey         02/16/2011     Contact

Philly 3D

    1. Physical Parameters:

If we let m=lb, Work=ft-lb, Horse & Watts Power=ft-lb/sec, (ft-lb/sec) x sec = KWH Work=(ft-lb) again.
Work <-> Energy = Force x Distance = mg x ft = lb x ft/sec^2 x ft = lb x ft^2/sec^2 = lb x (ft/sec)^2 = lb V^2.
Kinetic Energy mVt^2/2 as terminal velocity Vt is bought to 0, we have avg V = Vt/2 or m V^2 again, m=lb.
Pressure P = Force/Area = mg/ft^2 = m x (ft/sec^2)/ft^2 = m x (1/sec^2)/ft m/(ft sec^2).

Now if we apply a Press P for a distance (ft), the distance a building moves downwind before it stops, we have m x (ft/sec^2)/ft^2 x ft = m x (ft/sec)^2/ft^2 = m V^2/ft^2. Then times the face area (K ft^2) perpendicular to the wind, m V^2/ft^2 x K ft^2 = mK V^2.

            A. Facing the PSI:

An approach to this rough estimate of PSI is attempted here. When a wall of wind hits a standing wall, the air can't escape downwind but compresses to build up a standing pressure wave which takes the path of less air pressure resistance. This would be laterally or vertically to the sides and top of the building, but not through the ground. Consider a 1'^2 area with an adjacent other 1'^2 area also trying to have it's overpressure escape. So each air volume pushes against each other and a balance occurs. The other means of escape is upwind but that's against the steady flow against it.

Not only is a larger face area increases the total pressure on a building, for the flow to move laterally or vertically is inhibited as only 1/2 the say 1^2 area areas along the top and edges can escape, the PSI builds up more in the center. Also, with a wider face dimension, the support to hold the wall are further away at the sides. Thus, a larger pressure buildup in the center of the face wall that has the least support. However, the higher buildings increase the face area also but increases the torque by leverage on the lower wall section also. Thus by wind, it's center and downward section should blow out first.

We don't have such pressure gradient & torque factors to include here but know that it is somewhat proportional to the width & height increase for the PSI and the height for the torque. At least as the face area increases, both these effects increase as well as the total force against the face wall as it's area increases.

            B. Wind Velocity <-> PSI Estimate:

Let's assume that the air has no way to escape around the edges and top of the face wall. Therefore, we would have back pressure upwind until a balance occurs. This might be analogous to holding a large drinking glass into the wind. So, considering 150 mph, and now the fudge icing on the cake is to estimate how long until this pressure balance is achieved? Let's take this piece of fudge as 1 sec, but probably less. We then have a deceleration term as

d = (150mph-0)/sec = (150mi/hr)/sec = [(150mi x 5280'/mi)/(hr x 3600 sec/hr)]/sec = 220'/sec^2. Well the weight of air is 14.7 lbs/in^2 and air pressure is given by mg/in^2. Well the m is constant to both normal and horizontal so the ratio is the deceleration / acceleration as seen in (md)/(mg) = d/g. Using g as 32'/sec, we have 220/32 = 6.875 times as much pressure as the weight of the atmosphere itself. That's x 14.7 = 101.1 PSI. Overpressure (>14.7) = 101.1-14.7 = 86.6 PSI. Water = 62.5 lb/ft3, that's 62.5/144 = .434 psi for 1 ft of water. So this is equivalent pressure as diving down to 86.6/.434 = 199.5' of water.

Now we can see the destructive power of a class 5 hurricane or tornado. With just a little wind gust, umbrellas are reversed. However, sailboats trap to make good use of PSI's. In our building face, 2 PSI << 86.6 PSI or only 2.3% of 86.6 = 2 PSI meaning 97.7% of the overpressure would escape around the building acting as a fluid flow.

            C. Apply to WTC's:

For the WTC, it was supposed to withstand a wind velocity of 150 mph, if I remembered my research correctly. At that point, lets say that it crumbles because of the lateral shear stress. If it had a skinny side such as the United Nations Building, and the wind was perpendicular to it's length, then the thin width downwind doesn't support it as much so may topple. With the wind was perpendicular to it's width, the reverse is true.

      WTC      height=1365'       length=208'       width=208'

Now we substitute, bh = 1365, length = 208 and V=150 mph, and m = density or weight of air that is relatively constant, K = 283920 x 150^2 mph = 6.3882 x 10^9. If we have a higher height or a longer length with the width remaining constant, the building should topple. If we have a wider width backing the face, the building should stand.

Now we didn't change mph to '/sec but this would just make a bigger constant. We also have the width = 208' to stop the building from toppling. This could be expressed as 6.3882 x 10^9 x 208/bw. Here we have bw=208 so 208/bw = 1. So for bw < 208' backing, K V^2 > 6.3882 x 10^9 it topples. For bw > 208', K V^2 < 6.3882 x 10^9 it stands.

Not knowing all the details of stress because we don't know how far each part of the building would move until it collapses, we can't calculate the PSI to topple it. However, we know that K V^2 is proportional to PSI. What do we want to assume the PSI on the building face would bring it down? A review of 1 Megaton PSI Effects again, we chose the 2 psi. Noted that most higher buildings are gone, we have the height factor that increases the psi area so this could be the dividing value for regular buildings that also were designed with less elasticity and could not withstand 150 mph if constructed the same dimensions as the WTC's. The near 2 psi overpressure would correspond to WTC face of 1365' x 208' = 283920'^2. So now we have psi/2 x 283920'^2 as our toppling barrier. That is if psi=2, psi/2 x 283920'^2 = 1 x 283920'^2 = 283920'^2. To include the bw term, gives us psi/2 x 283920'^2 x 208/bw as our formula.

Let this be toppling power estimate be denoted by Pe and let Pc be the critical value of 283920'^2, so we have Pe = psi/2 x bl x bh x 208/bw = Pc = 283920'^2 until a better estimate comes forth from the cement construction industry. So we have psi by the logarithmic formula by distance from ground zero, and the given bl, bh and bw building dimensions. If Pe < Pc stand, Pe > Pc topple.

            D. More Evidence of PSI in Action:

The term PSI is overpressure above the normal atmosphere. This is usually 14.7 psi (lbs/in ^2). A standard mercurial barometer has a supply of mercury in it's cistern base connected to a vertical vacuum tube. Why? Mercury's high density allows it to rise in the tube only near the standard of 29.92" at mean sea level (MSL). High & low pressure systems increase or decrease this height. The highest I have observed is about 30.65" and the lowest about 28.80". I'm sure I may have missed more extremes but other locations probably experienced more. During the winter of 2011, there was a pressure reading of 31.10 on our US weather map.

Notice my difference is 1.85". Now hurricanes have had pressure readings near 27" within the eye. If a hurricane approaches and say we had a pressure of 30", the 3" difference is like the barometer cistern and tube but on a larger scale. The hurricane's vertical eye is like the tube and the area outside is the cistern of ocean water. The tube isn't a vacuum (0.00") but 30-27 or 3/30 = 1/10 or roughly 10% of one. This allows a welling up bulge of the ocean that travels with the eye causing a storm surge as it approaches over land.

Another case is the tornado. In it's eye, usually 1/2 mile or less, we have even higher winds but the same principle. We introduce a new term dp/dt, the change in pressure per unit of time. For a hurricane, we may have 3"/3 hr but for a tornado, it may be 5"/1 min. Over water, we see the upward surge like the hurricanes produce, but here, they are called waterspouts.

If air at regular pressure is trapped in a house when a tornado passes over, even though the violent winds could blow the house down, the sudden dp/dt blasts the house apart like an explosion. This is why they would want you to open the doors & windows first so the pressure inside the house would escape to the lower pressure outside to equalize the pressure difference. If the house were air tight, the house would explode outward. This is some relief considering that flying glass is directed away from inside the house. Of course, you have the violent winds again for a short time. On inspection over farm areas, they found cows killed from internal bloating as the trapped air inside could not adjust that rapidly to the ambient air pressure. So the worst thing one could do is to take a deep breath and hold it just before a tornado passes over.

            E. Slow or Fast PSI?:

In the examples above, we have some crude estimates of dp/dt. It was 3"/3 hrs or 3"/180 min to 5"/1 min. However, the energy of the hurricane eye may cover 10-30 miles while the tornado about 1/2 mile. However, the dp/dt ratio is (5/1)/(3/180)=300. So the faster dp/dt causes more extremes locally. Now consider the pressure wave from a blast. Here we have the dp term as 2.0 over the ambient air standard of 14.7 psi. This is the same range of variation we have observed over 50 years and approximately the same with an advancing normal hurricane. However, the blast's dt term is inversely proportional to it's travel speed.

This is a compression wave that travels at the speed of sound (mach1), 761.2 mph. Assume the wave length is 1/10 of a mile, probably less but at 600 mph=600 mi/60min=10 mi/min, D=rt, t=D/r = 1/10 mi/(10 mi/min) = 1/100 min. So dp/dt is 2"/(1/100 min) = 200"/min. The ratio to the tornado is 200"/min /(5"/min) or 40 times as great. Obviously, using the term minutes is poor because the pressure jump won't last that long or we would have a vacuum within 5 seconds. So seconds would have been a better choice but the 40 times as great would still be the same. Actually, it's 761.2/600 x 40 = 50.7 times as much. Thus the overpressure psi of 2" should probably bring down most buildings, especially higher ones.

            F. Lower 2 PSI Already?:

With the 1 Megaton PSI Effects again, the building heights of multi-story buildings is the factor. How shall we classify such? Probably a 10 - 20 story? As of right now, we settle on a 10 story with a 2 PSI or maybe a 1.4 psi @ 20 story? This might be a fair guess because given a constant width, the PSI x Height here is 20 & 28 respectively? The second choice seems good until we find a better estimate. However, what would be a standard width for buildings?

            G. Ramifications of a Lower Limit PSI:

                (1) Outer PSI Ring:

As we notice the last line of inflow.bug , it is the last psi ring to make it into the stem's radius updraft @ r=5.259 mi. So, the psi @ r=5.26 mi should be calculated. Actually, let's do it for a range of r's from .05 -> 10 mi for every .05mi. In review of rpsi.gif, we had a better fit with 29 instead of 28.5 with a corresponding -1.73 instead of -1.653. We now have psi=29mega*mi^-1.73. See the match in PSI <-> mi. So for r=5.25 mi, psi=1.646. Thus, if we choose psi=1.4 for a 20 story building, this would be at 5.75 mi. The building could fall but never be sucked into the Nuclear Cloud's stem since 5.75 > 5.26 mi. We could also adjust the 1.4 psi x 20 story = 28 to 28/1.646psi = 17.01 for a 17 story equivalent.

                (2) Inner City Building's PSI Value:

Since r=5.26 mi is our outer limit, we only have to consider buildings within this radius. Most cities are bounded by ellipses, circles or rectangles imaginary perimeters. Most all the raw building data that can be obtained is within 5 mi^2 except NYC is a rectangular 9 mi^2. Now the r=5.26 mi gives us an area of pi*r^2 = 3.14 * 5.26^2 = 86.9 mi ^2. This means that all buildings > 17 stories would fall that are tabulated within 5 mi^2. This is an equivalent circle of pi r^2=5, r=(5/pi)^1/2 = 1.26 mi which corresponds to nearly 19.713 psi by our calculated table seen in PSI <-> mi. again.

                (3) Represented Width:

We have obtained 50 building samples within Chicago and we consider 17 stories x 11'/story = 187'. The proportion of height to width is about 7.5 near 200'. So 187'/7.5 = 24.9'. For height to length, it was 5 so So 187'/5 = 37.4'.

                (4) New Pc Formula:

We had Pc = psi/2 x bl x bh x 208/bw for WTC but now replaced by Pc = (psi/1.646) x (bl/37.4) x (bh/187) x (24.9/bw) where bw is the down wave support equivalent calculated by the topics below. Note that the face area is bl x bh and the down wave supporting length is bw. So smaller width and larger face area = larger Pc. Simplified, Pc = psi (bl x bh)/bw x 24.9/(1.646 x 37.4 x 187) = (psi bl bh)/bw x 2.163e-3. A check shows with psi=1.646, bl=37.4', bw= 24.9' and bh=187', (psi bl bh)/bw x 2.163e-3 = (1.646 x 37.4 x 187)/24.9 x 2.163e-3 = 462.32 x 2.163e-3 = 1.00 since (1.646/1.646) x (37.4/37.4) x (187/187) x (24.9/24.9) = 1 x 1 x 1 x 1 = 1. So now, if Pe > Pc=1, building falls, Pc < 1, building stands. Since psi = 19.7 near the edge of most city limits, then 19.7/1.646 = 12.0, then a matching 17 story building could be reduced by 1/12 or 17/12 = 1.42 stories to fall. That is, some 2 story buildings may stand within center city but are there many? Most all are skyscrapers > 2 stories and even if any, their contribution to the total dust from skyscrapers are minimal.

This means, we would bring most all buildings down within typical center city radii. This topic is expanded to discover this below since, some low level really elongated buildings oriented with it's width perpendicular to the PSI traveling wave may still have to be considered as standing by use of Pe, but are there any? This would be a poor use of real estate within center cities. PSI face vs. support is developed in the following topics below.

    2. Calculating Building Dimensions to Remold:

Now to compute face areas of the psi wave perpendicular to the building, view Building - PSI Wave Orientation. The building outline is in red with the psi down wave support behind the face in blue. Thus, ground zero is 270 degrees from the building. As seen, the length (bl) components support against the psi pressure wave. The top triangle (3 white vertices that use x1,y1 & x2,y2 is congruent to the bottom triangle (3 white vertices that use x3,y3 & xo,yo). The cwy face (cross wave y) is perpendicular to the direction the wave front travels. These 2 cwy's overlap. The total face length perpendicular to the wave is the y between the extreme corners pairs y1-y3. With a different orientation, the extreme corners may be y2-yo or the y difference between diagonal vertices. The vertices pair will be the even or odd set, differing corner subscripts by 2. The psi pressure wave is all along the perpendicular component length y1-y3 but not all at once. It hits different parts of the building at different times as the wave front advances across the building.

That's enough. To do more is a bore, after we thought this through with a breakthrough. If you want to see the original grueling analysis to calculate the 3 triangle areas that is superseded below, view Building Area Components .

    3. Weighting Factors:

Now we can substitute these values and weigh the blue lines as areas, but other orientations and change of signs cause negative areas and overlaps cause other problems. Noticing that we are summing the 3 areas that the blue lines make, then the summing all areas between the blue line is nothing more than the total area of all 3 areas (2 triangles & 1 trapezoid) itself which is the complete area (ba) of it's rectangular roof or base.

Now just consider the length of all the blue lines and averaging would also be a good estimate. This would be 21 lines, including the line length=0 at the bottom vertice. All we need to do is find the length of the maximum vertical line (y1-y3) normal to the pressure wave direction and let this be a fictitious side (new length or width) and all the blue lines normal to it can be simply averaged across this cross wave dy (cwy) distance or simply called the cross wind width (cw). Each of the 21 blue lines are 1 pixel high and (y1-y3) is say approximately 400 pixels. Now suppose we double the 21 lines to 42, double again to 84, but then triple to 252 lines. Here we have over half the rectangle filled with blue lines. Now if we double again, we overlap some of the 504 lines as the complete rectangle is filled with blue. So, if we fill all 400 possible pixels with the 400 lines, then the sum of all the blue lines make up the base building area ba pixels. So the sum of these 400 blue lines by cw=(y1-y3) is simply the area (ba), then it's average length is ba/(cw=400) or called the downwind building length (dwbl). So with this one dimension length (y1-y3)=(cw), we have the other equivalent dimension or dwbl=ba/cw. Likewise, cw=ba/dwbl. Here we have the new downwind support term (dwbl) and vertical (exposed to wave face cw) components alignment to the pressure wave front direction of travel. The rectangular building has the same area with any orientation, but as far as the pressure wave acting on it is concerned, it's impact is on these new fictitious dimensions due to it's different orientations.

    4. Computer Graphics to the Rescue:

Well if a rectangular building's orientation and rotation isn't such a potential problem, how about non-rectangular buildings? This is a whole new ball game with dimensions, area and alignment calculations. We let the computer calculate the building areas by counting colored pixels. We also need a method to handle multiple vertices instead of the 4 corners which is a subset of multiple vertices.

        A. PSI - Building Alignments:

            (1) Building Parameters:

With the vertices read in, we can calculate the building's length, width, area, center and alignments to true north, direction and range to ground zero. Most important, the vertices must be read in order of rotation or else calculate the order. For a 4-corner building, we would have
dxl=x[2]-x[1], dyl=y[2]-y[1], bl=sqrt(dxl*dxl + dyl*dyl) where l denotes length. With w for width, we have
dxw=x[1]-x[0], dyw=y[1]-y[0], bw=sqrt(dxw*dxw + dyw*dyw). An added feature that we might need is the diagonal diag=sqrt(bw*bw + bl*bl). For the building alignment of the length to true north, we call a subroutine vect() which is deg=rad*atan(dy/dx) or the direction of longest adjacent 2 vertices for 4 corners. First, a correction for using only radian degrees, and a trigonometric to geographical direction switch by if(dx < 0) dir=90-deg & if(dx > 0) dir=270-deg.

However, for the multi vertice corners, we have to calculate a representative bl by comparing the distance for every 2 vertice combinations, not just the adjacent's. The maximum vertice separations with their distance & vertice # is saved to calculate the representative alignment by vect() again. The need for this alignment is not apparent yet but calculated now rather than all this over again later.

            (2) Building's New Coordinate System by Wave Direction:

What is needed, is the building's alignment to wave direction. Defining the coordinates of ground zero as xgz,ygz we use the middle of the building xmid,ymid. For the 4 corners, this is the average of the opposite corners - the mid value coordinate of each diagonal which intersect at this same point. For the multi vertice corners, this may not be the exact center but close enough considering that the building dimensions are in 100's of feet but a few miles to ground zero and considering the closer in buildings are all demolished anyway. Then with dx=xmid-xgz & dy=ymid-ygz, we use vect() again but add range as rng=sqrt(dx*dx + dy*dy) because we need to know this distance so we can calculate the psi at each building's face. We now transfer 'dir' to the new variable name wavedir.

The coordinate system must be relative to the direction of the pressure wave axis. For reference, but not apparently needed, we have the difference between the direction to ground zero and the building alignment wadir=(dira-wavedir). We apply the formulas of the grid rotation around xmid,ymid as follows:

            x[n]=(xn-xmid)*cosd + (y[n]-ymid)*sind
            y[n]=(y[n]-ymid)*cosd - (xn-xmid)*sind

for each vertice n. We have to hold the old x[n] as xn before changing x[n] as it's used for y[n]. The cosd & sind is for the angle ddir=270-wavedir or the angle to rotate back to the original 270 deg coordinate system. Applying radians again, cosd=cos(ddir/rad) and sind=sin(ddir/rad).

            (3) More Use of Wave's Coordinate System:

As we rotated all the vertices, we kept track of the rotated max & min x,y and their corresponding vertice number n. After the rotation is done, we have rdy=ymax-ymin - the rotated building cross-wave (cw) max dy and rdx=xmax-xmin - the rotated building along-wave advance direction max dx. Then we finally have the prize, the dwbl=ba/cw - the building backup 'length' supporting the front face cross-wave advance direction across max dy. Well we know the area (ba) for a building rectangle as ba=bw*bl but not for multi vertice non-rectangle building.

            (4) Computer Graphics Breakthrough:

To handle the odd shape buildings, let's first plot the building's outline then fill it's interior with the same color as the outline and let the computer count all these colored pixels. Since we want to be as accurate as possible, instead of plotting each building as a different size image by a fixed pixel to length spacing, let's maximize each image to fill the screen for more accuracy then adjust accordingly. We need to fix a proportional standard example to use. We define a dpdx0 as dp/dx and dpdy0 as dp/dy - the '0' signifying the initial image. We also need to keep a constant aspect ratio for pixel spacing as arp=(swi/px)/(shi/py). The shi = screen height inches or some units and swi = screen width. The px,py are the screen's resolution. Keeping the graphics simple and sufficient for most CPU's, we keep the basic px,py resolution as 640,480 pixels.

For our initial '0', we have ba0=16800'^2, dpdx0=2.18168, dpdy0=-2.27572 and counted the colored pixel area as pixc0=83927 pixels. All these are applied by ratios and explained in specific details with Building Areas by Pixels . This is to make sure the scaling concept and values are applied correctly to the denominators or numerators. With all this explanation, we come away with pixcs=pixc*(dpdx0/dpdx)*(dpdy0/dpdy), bap=pixcs/pixc0*ba0 and wbw=bap/dwbl.

                (a) Computer Graphics Examples:

We bring in our screen capture subroutine scap() to capture this action. Our initial capture is Building # 1- Rotation # 1 . The solid darker blue area is the initial (0) building of 83927 pixels (pixc0). The lighter blue or cyan is the same area but non rotated with the area equivalent width of wbw. That is, the dark blue areas outside the cyan box should fit into the black area inside the cyan box. Think of the dark blue area as silly putty re-shaped to the cyan colored box.

Note that we have added the red PSI over pressure wave. There will be a different PSI value for each different radius of curvature at each building's location from ground zero. A 5-degree more CC rotation is seen in Building # 1- Rotation # 2 . Finally, we have the non-rectangular multi vertice plot in Building # 2- Rotation # 1 . Here we see the calculated maximum length for building orientation (yellow line) that was explained above. We check out all programs with listing of parameters and variables values by a bug file. For these last 3 pics above, the detailed info is seen in BPSI.bug , only for the hardy at heart.

    5. A Great New Purpose:

Why do we need this update? As Google Earth was used to locate the highest buildings in urban areas and count the building's stories as shown at the top of this page, the raw data to make up these building displays could be further used as a more detailed dust inventory if we consider each building in a city. Instead of having a representative building height for each of the 70 PSI rings, we calculate the % dust and % particle distribution for each floor slab for each building that is assigned to it's appropriate dr=.1 mi psi ring as explained above. Here, instead of incrementing each ix psi ring and approximate a common building height with fixed dimension for each building then estimate the number of these buildings in this ix psi ring by a density guess, we calculate each ix psi ring by range to each building from ground zero as well as their real dimensions and even their orientation for psi impacts.

    6. Lidar - Light Detection & Ranging Data:

Above, we had obtained 50 building samples within Chicago mentioned above. Some of the data was questionable and we tried these methods just explained above on the sample. Also, their complete costly data did not extend far enough from center city. After poor results trying to obtain more data and for free, it appeared as a dead end. Finally, the USGS has Lidar Data for 12% of the country. After struggling with their header file and data in LAS1.1 format, it was discovered no buildings were classified. So this data does not indicates buildings but a 'bare earth' or terrain type as we first thought. Later below, we bring this all back.

        A. Original Lidar Data Source:

Our National Hurricane Center in cooperation with Florida International University had obtained some Lidar Data from 2002-2004 but in another format. No big header data but 5' square grids with heights. However, they also had the matching grided data for the 'bare earth' so the top of building height - bare earth heights = actual building height. However, a new coordinate system called the State Plane is used. This is Florida East tiles that are 1001x1001 5' cells. More details can be seen in Metadata_topsurface.

        B. LIDAR.c = LIDAR building heights :

Our problem is, we must group all the adjacent height points together. The minimum base height for a building was chosen at 25'. A 2-story house could be usually wooden frame. Any higher, was assumed an office building of brick, cinder block, stone or cement mason. The name of the 1001 x 1001 is given by the state plane coordinates of the SW corner in 1000's of feet. The file t920_515.flt means Top heights @ state plane's SW corner at 920,000 East & 515,000 North. First, we read the building height data of 4 bytes from t920_515.flt in IEEE format. We touch the conversion method here. The 4-bytes are 32 bits but reversed order from our normal number system (Big Endian) to (Little Endian). This allows the computer to access the lower values first as it always has to. The highest value to the right are varied so may be blank. This depends on your PC's CPU chips. The 32 bits are converted to a floating point scientific notation with base mantissa with only + exponents. The last bit is the +,- value of it's converted value. Next, read the ground height in MSL from b920_515.flt and subtract it from the building height - leaving us the building height above ground.

Now the fun begins by grouping all the adjacent height points into one building. Well, we search for the 1st point and try to find all the adjacent points on the same row that touch each other. Keeping track of the 1st & last point, we replace those found points with 0 so they won't be read again as another building. We make use of the file position (filp) by using the fixed record length (recl=24) of each building point for rewriting 0's after it is used. For each adjacent ix or iy, we sum the building height, ix's, iy's and the building count (nc) so we can calculate the average position and height when all points are tabulated. Also, from all the labor from the previous topics, we keep track of min & max ix,iy's to determine those 4-corners touching the x,y axes just explained above.

We set a rule that the next row's 1st valid point can be 1 less than the left and the last 1 more than the right of the previous line. For more details of this procedure, see IyIx.dia . Once all points of a building is complete, we originally wrote out the tally() of it to 920515.9bb file using 9 bytes describing the building's location & dimensions.

        C. Major Data Upgdate: 02/14/2012

As this date indicates, the date at the intro of the page indicates 1 year of activity, or in some cases inactivity due to personal distractions by house activity. A dead squirrel started my attic search till found that initiated an elapsed 30 year attic cleaning project and repacking activity. Also, a new bedroom floor caused by a bad rain gutter seeping water under the house from excessive rains. Thus a French Drain was installed taking 3 weeks to complete then a part of a new roof.

At least, this was in conjunction waiting for a response from Google Earth requesting the sharing of their 3D City Building Data that they advertise so well. Many letters explaining that I would only need to run their data once and not need it any further so as not to compete with any of their displays that would compete with theirs as this would violate their user's agreement contract. All they would say, if at all, was review their user's agreement. Actually, it did say in their agreement that one would need permission otherwise. However, the permission goes into one of those link loops that ends up where you started. Trying other personnel, I explained that I could run my program on their data at their site location and then forward the results. All the buildings would collapse into the dr=.1 mi radial rings so all the details of the buildings are washed out into cumulative tons of dust. Also explained, this is the best use of their 3D building data ever devised for GIS and is a great patriotic effort. If any one has a better connection with Google Earth than I, please contact them but I will try even harder for their data during 2012.

          (1) 4-Byte Original Lidar <--> 64-Byte Lidar LAS1.1:

Since we explained above that the newer & better Lidar LAS 1.1 data never classified buildings, we exploited the older 4-Byte Lidar. After successfully converting the newer 64 bit Lidar LAS 1.1 data, with large headers, the older 4-Byte seems a snap. The converted values looked good until I cross checked some data points with Miami's FIU that isolated a 5 mile radius around Miami and gave me a ftp site to download nearly 100 top and ground base tiles. Using a tile, I discovered that my conversion was off. They had run the data on their ESRI-ArcGis-ArcMap. After trying to resolve this 'bit' by 'bit', so to speak, no solution was found. My main contact went to Germany from Dec '11 - Feb '12 so could help me no more. In the meantime, I downloaded this gigantic ESRI-ArcGis-ArcMap10 for a 60 day trial program. Then I ran their toolbox for the data conversion Geo processing>ArcToolbox>From Raster>Raster to ASCII on the data. To my surprise, it didn't match FIU's values, or worse mine either. However, mine was closer, 1.5' less than FIU's.

                (a) Bite the Bits:

After all possible combinations of bytes and bits, see ieee.dia for more details, I could not reproduce either conversion, FIU's or my ESRI-ArcGis-ArcMap10. Talking to ESRI, I could not get unlimited help and whoever wrote that code was probably somewhere else and they don't use this data too much any more. However, I found a later version of Total Image Converter that I use to convert 100's of PCX files to GIF's by DOS command lines for Google Earth Displays as seen in previous pages. They converted this older 4-byte Lidar data (*.flt) but made a tile of about 36 pics, with blocky interpolation for every foot between the cell size of 5' and added shading features. This is obvious as seen in one of their generated subtiles T920_52005 . Obviously, this wasn't a simple file conversion but added bells and whistles for pic displays that only ESRI knows. The help desk of Total Image Converter is trying to contact the vendor but has been unsuccessful. Not knowing the correct value, we left off with the older 4-byte Lidar temporarily.

                (b) Return to 64-Bit Lidar LAS:

As I had originally ran this before, I noticed that the tile # run was not the real high rise center of Miami and the numbers really appeared as terrain. Some of the data was higher but was not classified as 'ground' but as 'unclassified'. So, after downloading a another tile (117542.las) and finding it's State Plane within the file (920E520N), we found many high rise values but they were still 'unclassified'. So, it appears that instead of doing all the classification work to reclassify points as buildings, it was easier for some end user to do this. However, they did classify ground & water. This means we could use this but sift out trees, shrubs & highways. A sample of 54 tiles with headers and min-max x,y,z values is seen in lidarlas.b54 . It is noted that the some minz & maxz values are out of range for this Miami region so filtering code is also required.

        D. Processing LidarLas1.1:

For a quick reference of all these consecutive programs, we can see the flow diagram in lidarlas.fld . We started with LIDAR.c above but now we modify it for the new reading format code and procedures as LIDARLAS.c . Notice, we have 2 choices of inputs for the unzipped & renamed *.las formats. We have the local C: drive but since these files are near 200 mb each, we have them stored on a portable external hard drive.

The las.dir file is the basic alphabetical *.las directory file listing. The geo.key input file is the tabled header info just seen above in the lidarlas.bug file. Since the State Plane coordinates are within the LAS file, we need to link Tile #, Las # and the SW corner (SwxSwn) of the tile file for later use. This is seen in the 54 tiles of Tile.spc . Using the SwxSwn, we have our 2 main output files 920E520N.grf & 920E520N.grg. The 'f' stands for # of floors in a building and the 'g' stands for ground elevation grid files. Instead of arranging to sequence all this data for x,y coordinates as in the older Lidar, we can put each point into it's proper grid position by a file position filpg=idy*1000+idx where idx & idy is the number of points E,N from the SW corner respectively.

However, noticing the # of data points (npr) in the bug file above, it is somewhat near 7 mega points. Allowing for one point as a pixel, we would have a 7 megapixel pix if no compression was used. Since we should have a grid file range of 5001' x 5001', we have about 25 megapixels. Even though the resolution space between points is sometimes less than 1', at other scan of rows or columns, it may be well over 1'. Obviously, a 25 megapixel area/7 megapixels points = 3.57 '^2/pixel. So 7 megapixels won't fill the grid and should leave some holes as seen in a cropped portion of the .grf converted to a pic by LAS1.1 Full Grid Cropped . Use the magnifier to zoom in to see the spaces. To resolve this issue so no gaps are present in order to have continuous building points, if we let 1 point represent 7^1/2 = 2.65', then we would exactly fill the 7 grid megapixel grid. However, since all points are not linearly spaced, a close examination of LAS1.1 Full Grid Cropped again shows some areas with at least a 2-pixel separation, so we need at least represent 1 pixel point = 2'x 2'. However, for a larger safety cushion to eliminate holes, we might as well go to 1 point = 5' by squeezing the grid down 1/5 in each dimension. Thus, idx=(x-minx)/5+.5 & idy=(y-miny)/5+.5 with x & y the number of feet E,N from the SW corner.

Now instead of holes in the building's grid, we are more likely to have overlap points. So, for building floors & terrain, we read the stored height value and take the highest of this or the newly overlapped floors or terrain point. For the next record, we jump over the header data and go directly to the x,y,z data by file position filp=offtpd+(nr-1)*pdrl and then it's classification by filp+15. Now pdrl is point data record length from the header file.

        E. Grids to Graphics:

            (1) TerPCXK.c: & BnfPCXK.c

Now to make sure that the data is converted, we want a graphic display to act as a Google Earth overlay to see if the pic tiles superimpose on top of Google Earth's 3D buildings. For building # of floors, we have BnfPCXK.c and for Terrain, we have TerPCXK.c . These programs are identical except for the names of the input & output files best seen in lidarlas.fld again. The value from the grid is a 1 character byte from 0 - 255 in value that takes the Building's Number of Floors (nf) in the 1st and elevation (z) for TERrain height in the 2nd programs respectively. The range of values are fine for all buildings in the US but 0 to 255' is for Miami terrain only. We may need less detail for higher terrain in other US cities. However, we restrict elevation to 0 instead of feet below sea level. These were small in ground areas around Miami but some others below sea level could be wave action and low tide marshes? We would not have the base of a solid building under water anyhow unless on stilts or piers.

Both programs have a file directory file, *.dir for input for the *.grf & *.grg files and the generic combo.kml file for Google Earth displays. Both have an optional output *.bug file, the *B.pcx & *T.pcx picture files as well as the *.kml files for Google Earth. Both have the same pcx.sp State Plane coordinates file that are the 4 corners of the tile. We pause the program to jump to the open window of Corpscon 6.0 to convert the 4 corners in pcx.sp to pcx.ll. It would be best to run both programs parallel so when each program continues, both can use the same pcx.ll, the 4 corners in Lat/Lon for Google Earth's *.kml files. The Total Image Converter must change the *B.pcx -> *B.gif files and the *T.pcx -> *T.gif files by their appropriate batch (*.bat) file listed.

            (2) Display Results:

Looking at the pic conversion before superimposing on Google Earth, let's examine our 920520B Tile. As the lighter reds are the higher buildings by # of floors, a zoom in shows no holes or black pixels within the buildings. However, we see the drop off (darker red) at the edges of the buildings indicating that those point areas used can be used as ground elevation in the corresponding Terrain Tile 920520T Tile. Here we see the inverse, the black areas indicating a building structure, although not officially classified as one, or water. Again, the lighter the green, the higher the green grass or ground area. Also notice back in the building pic, the bright red lines adjacent to some buildings. The only valid explanation is the higher steel construction cranes that have small but excellent reflection surfaces that were probably building onto those structures during that time period mentioned in their .xml file.

Some other interesting features are seen in 885525B Tile and 885525T Tile. Here we see a causeway in red built upon the plowed up ground mounds across the bay. In 885540B Tile and 885540T Tile, we see a built up high ground foundation in an otherwise residential neighborhood indicating a large building. Notice all the smaller blotched black areas indicating that this is a residential neighborhood of houses. We see this overlaid in Google Earth later.

        F. grtsbfpT.c & grtsbfpN.c = Grid To Single Building Floor Pics:

            (1) Corner the Buildings:

Now the difficult task - to capture each building depicted in the tile pics above. We use the *.grf grid file from LIDARLAS.c above to isolate each building and correct it's msl height to height above ground. We use a directory file list in bnf.dir to list the file inputs *.grf's and ter.dir for *.grg's. We have a file that links the Tile # with it's SW corner state plane in SwxSwy.til. Now each building will have it's own coordinates as state plane that can be used directly in the next program but also to convert to Lat/Lon for Google Earth overlay displays. This accumulated file 'pcxb.sp' is appended with each building's 4 coordinate extremes when extracted rather than convert each building's Lat/Lon when it occurs as we do for each tile. With possible 100's of building for each tile, we do them later as one big batch.

            (2) Procedures:

From the # of floors *.grf grid files, we construct the file which represents the delta Y and X Left & Right building limits. This file collects all adjacent building points along a line. It finds the 1st left building point then all the right consecutive points for each building on the same line. Then it does the same on next line and finds it's adjacent section there continuing as the previous line until no more adjacent building sections on subsequent lines. After is complete, we reread it and gather all building sections on one building at a time. Since we have the horizontal dimension on each line, we ignore all other building points on the line and find all the adjacent building structures for consecutive vertical lines down until we complete one building. We record a done=1 in the 1st column of the line so we won't read the same building again when we return to this first point on the first line of the just completed building by the initial file position filp0. Using filp0 saves CPU time by not rereading all the previously used points.

            (3) File -> Sub Files:

Viewing a partial of indicates we are collecting data to the last line, not shown, in the file to complete the 1st 3 buildings that were started on the first line. Notice the done=1 column after it was used. After completing all buildings, we should have done=1 in all lines of the first column. We carry over a 'nfr' term which is the total Number of Floors per that Row or the sum of nf's for each row. With each line of accepted as part of this isolated building, we copy this line without the done=1 & nfr to . Notice here, the 1st 3 buildings left to right consecutive solid points are aligned to the same building's line above it. The ixr-ixl indicated the color run for the *.pcx file while the ixr-ixl+1 also indicates the # of cells (nc) giving us an area count for each line and a total area for the building after completing all of the building's scanned lines. The comments and the : (continuing data) were added to the file for clarity. The ----- indicates a separate file for each building and 'eof' means end of file.

            (4) Building Parameters:

A seemingly trivial note about nc=ixr-ixl+1 rather than nc=ixr-ixl is expounded in EndPoints.dia . For the *.pcx picture file, we scale the red buildings from 0 to 254 @ 72 stories for Miami. Also, we need to correct the nf to nfa or the # of floors adjusted (nfa) by terrain. That is, having the terrain height z, we have rnf=nft/nct or the Real Number of Floors by the total area cell count (nct) divided into the total # of floors (nft). Then nfa=rnf-z/11+.5. To determine z, we had opened *.grg terrain files and matched the perimeter of the assumed rectangular building searching for a valid z. This worked well but didn't work for all cases which missed about 3%. So now we return to the old Lidar (32 bit) and assume that ESRI's translation mentioned above is correct. See Lidare.c in lidar.fld . In subroutine getz(), we read *.grg with the matching grid position (filpg) by averaging the min & max ix,iy's. This should be near the center of the building's base. We also condense the building parameters into 11 bytes that will be translated relative to ground zero in another program below. These 11 bytes are as follows:

    ixmin=2 bytes
    ixmax=2 bytes
    iyt=2 bytes
    iyb=2 bytes
    nct=2 bytes
    nfa=1 byte

            (5) Basic Ground Rule:

There are some basic rules for creating these building files. For reading the file explained above, if the next adjacent line below (iy+1) is from the same building's line above (iy), the mid position of that building run line above (iy) must fall between ixl & ixr of the next building's line (iy+1). Thus the code is ixmid=(ixr+ixl)/2 and if(ixmid >= ixl-1 && ixmid <= ixr+1) it is accepted. Notice, this allows for 90 degree corner walls sloping from the corner by extending ixl & ixr 1 data point.

            (6) Which of 2 Programs?:

In these 2 programs names, we have the last character an 'N' or 'T' which stands for Tiered or Non-Tiered buildings. In the tiered, we have a separate structure when we are > 1 floor from the central nf loop interval, We treat it as such but the tiered means we may really be connected to the same building. For the Non Tiered, we just have adjacent structures. Also, for the *.pcx files, we have them as TtileNfNb.pcx Tiered and NtileNbNfa for Non-Tiered. For identification and trace back, the 'tile' in both file names is the tile # (2 digits), the Nb is the Building # (3 digits) in tiered, (4 digits) for Non-Tiered. In the tiered, we have 'nf' = # of floors (middle floor) in the loop but 'nfa' is # of floors adjusted for the non tiered. To keep the files is a numeric order, for the tiered, we cycle through Nb with each increasing Nf so they are a consecutive numbering system. With Non-Tiered, we have no fixed Nf but Nfa which can be any number but the Nb is consecutive and can go to over 1000 buildings so this is a numeric order for the 1st 6 digits. For the tiered examples, see:

    1. Miami Tower
    2. Wachovia

This is all explained in detail by viewing grtsbfp.dia . Notice, we accept smaller areas for the tiered. Because of the larger number of building structures for the tiered, we will used the non-tiered and maybe adjust the volume upward by it's increase of 12%+. However, this may include trees so we dropped this 12%.

        G. PCXBKML.c = PCX Building pics to KML :

Using the pcxb.dir generated from the *.pcx files, and it's matching 4-corner pcxb.ll file from Corpscon 6.0 conversion of pcxb.sp from grtsbfpN.c or grtsbfpT.c just explained, we generate the tiered or non-tiered *.kml file for Google Earth displays here, such as seen in:

    3. Tiered
    4. Non-Tiered

By switching the displays in Google Earth, we can see if the composite individual buildings displayed here is superimposed on the tile pic 920520B Tile again by turning the 3D Google Earth Buildings on & off. However, as noted the black background around the individual non-tiered is larger because of the larger buildings. With the tiered, the tighter smaller background pics allow more of the complete building to show and not obscuring adjacent buildings.

Just for fun we modified PCXBKML.c to see what the old lidar from Lidare.c would look like since we are now using it for the *.ele terrain file for z to correct to actual floor heights (nfa). This is seen in 920520 Terrain Old Lidar Tile.

Well the first surprise was that originally, the picture was upside down compared to the LAS1.1. This means that the correct orientation in Google Earth showed which one was wrong. It shows that the old Lidar was correct. This means that the old Lidar, although named from the SW corner, the data starts from the NW corner. However, from the LAS1.1, it starts from the SW corner so going down the iy rows is really going up north. This is why the batch files for grtsbfpN.c, grtsbfpT.c, BnfPCXK.c & TerPCXK.c need it's images flipped vertically even though their Lat/Lon corners are correct. However, we are setting file position 'filpg' assuming 0 from the NW corner in these LAS1.1 files. This is OK as long as we flip it later but to determine z from this older Lidar, we redefine the filpg from the bottom of the file instead of the beginning top. See the explanation diagram in lidar Old-New.dia .

Another note from 920520 Terrain Old Lidar Tile is that the base of the buildings indicate a uniform terrain, which is actually done in the ground leveling before construction, but this scan didn't just occur at that time for every building. This means that for this early version, they just picked a fixed elevation across each building's base area but an unrealistic slightly higher height is assigned or else these bases wouldn't visually show up that well. We finally go to the National Elevation Data 'NEDS' base because of the adjusted heights under buildings (fudging) here. The link is at the bottom of this page.

        H. XB11B.c = Extract 11 Byte Buildings :

As nice as the pics are showing the wave, we are going to describe each building using only 4-corners of each building as we have access to better data now as just explained above. We will also show the adjustments for odd shaped buildings below. From now on, we will stick to the 4-corner rectangular buildings as we need to find the building's orientation's dimension by rotating the 4-points touching the axes until the maximum area reduces to the actual area. This is best explained in BoxedBox.dia . We extract the 11 bytes listed above, convert then list these final parameters as numerics described below.

            (1) Building Length, Width & Alignment:

We use the loop to find bl & bw but using substitution reduced the equations to find bl,bw exactly by
bl=x12*cosd where x12=x[2]-x[1] & bw=y10*cosd where y10=y[1]-y[0] after cosd=sqrt(area/(x12*y10)) was derived. Also, we have the building alignment by a=acos(cosd) that uses radians so dira=ddir=a*rad for degrees.

            (2) NESW Building Coordinate:

Using bl & bw, we found the NESW coordinates corners xr,yr's by

xr[0]=x[0]+bw*sin(a), yr[0]=y[0]
yr[1]=y[1]-bl*sin(a), xr[1]=x[1]
xr[2]=x[2]-bw*sin(a), yr[2]=y[2]
yr[3]=y[2]-bw*cos(a), xr[3]=x[2]

            (3) Ground Zero <--> Building Coordinate System:

One way is to now take the xr,yr's and rotate them more to the wave direction from ground zero and calculate the building area face towards ground zero. Another method is to ignore the NESW xr,yr's and set up a new coordinate system with ground zero's axis as xg --> & yg ^.

In order to see what is needed, we show Ground Zero East . Here, the pressure wave hits all the length side of the building from the east. Now we see in Ground Zero North , all the width is used. This is more resistant to the psi wave because of the down wave support from it's length as explained above. However, suppose that ground zero is not exactly east or north and the building alignment is Ground Zero 012 Deg . Ignore ground zero direction for now and just see a 340 Deg Building Alignment . Now we add a Ground Zero = 045 Deg . Assuming the arc of the pressure wave is nearly perpendicular to it's radius, especially at distances out to near 5.25 miles, the wave impacts the building first at radius R1 then leaves the building at R2. Now we take away R1 & R2 and show the Building & Ground Zero Orientation Angles then add the Building Corner #'s .

Now with a little imagination, we first pretend to have simple building coordinates by assuming that it was oriented W-E along it's length. We know this is probably a rare case but will correct this later. This is seen in the bottom diagram (below -----) of BoxedBox.dia again along with the rotation formulas. This is a match to Ground Zero East except the length & width are switched. Here, we actually have the building orientation angle = wave direction angle, assuming we are orienting by building width instead of length. If the wave direction and building orientation angle are both rotated to Ground Zero 012 Deg , the wave effect on the building is the same. However, if the building is rotated by a 340 Deg Building Alignment instead of 0 deg with the same Ground Zero=045 Deg , we have a longer cross wave face (cw) value indicated by the perpendiculars @ R1,R2 diagonals instead of just the building width.

            (4) Ground Zero or Building Coordinate System?:

The 4 corners of the building are analyzed relative to the ground zero point. Since the coordinate system is related to the wave direction axis, so then also are the building coordinates to this new axis of origin. Again, let's assume the 4 corners were set up as if the wave dir was 0 deg shown in Ground Zero North . Now if the real wave direction (WD) is an axis @ 065 deg, we see it as the red axis rotated about it's center o @ 65 deg in WD 65-119 . We could also say that the 4 corners are rotated -65 deg around it's rectangular center o in Ground Zero North as ground zero still remains north as shown in 4-Corner Rotate -65 Deg . Notice, the corners have replaced their +-bl/2,+-bw/2 with the new rotated xr,yr's coordinates that still reference true north. This relationship to the original north WD axis is the same as the previous +65 WD rotation. Regardless which we use, we have the rotation equations for the latter with this easier WD=0->3, 4 corner coordinates initialized by the 4 +-bl/2, +-bw/2 combinations as seen in Ground Zero North . This means we also can rotate these corners in Ground Zero North -119 deg if we want the new corner coordinates relative to the blue WD = 119 deg in WD 65-119.

Finally, we have the NESW view in Ground Zero 119 Deg , and just showed that if the whole system (picture) is rotated, no relationships between the ground zero axis and the 4-corner coordinates change. So if we rotate Ground Zero 119 Deg back CC 119 deg (-119), we have the appearance of Ground Zero 119 Deg Rotated -119 Deg. Here we see the ground zero axis pretending to be back to true north again with the 4-corners rotated -119 deg with new corresponding 4 xr,yr's coordinates.

In B340-W12 , we have extended vertical lines in order to depict the angles. Here angle 'a' is building orientation (BO) = 340 deg or -20 deg from 360 while 'w' is the wave direction (WD) to ground zero or +12 deg from north. Following the light blue line as the building's orientation intersects the wave direction, makes an angle 'd' as the needed direction difference to work with. We don't need to bother with a NESW coordinate system as the whole system could be rotated and the direction change result is seen in Rotate 35 deg , but the relationship between the building and the ground zero axis is still the same. So we have angle 'd' as the difference between WD -> BO or 12 - -20 = 32 deg. Recalling the BoxedBox.dia again, the xo,yo's were set by 1/2 signed bl,bw so regardless of the apparent rotation of the WD, it's the real rotation of xo,yo's that are rotated to xr,yr's in reference to the WV axis. As seen in the explanation at the bottom of BoxedBox.dia , we have

sind=sin(d); cosd=cos(d);
for (n=0; n<=3; n++) { xr[n]=xo[n]*cosd + yo[n]*sind; yr[n]=yo[n]*cosd - xo[n]*sind; }

            (5) Cross Wave Face - Wave Back Up:

We keep track of the max diagonals of y that is |yr[2]-yr[0]| or |yr[1]-yr[3]| for the final Cross Wave CW shown in purple. Then finally, we have the wave back up length (width here) wbw=area/cw, or dwbl as derived above.

We also have 2 other output files. The 1st is 920515.BCh where we output
BldCor     E-sp     N-sp     bh'     as seen as the record header file. This is a state plane 4-corner coordinate file to convert to Lat/Lon by Corpscon 6.0 to Google Earth. However, the pcxb.sp from grtsbfpT.c or grtsbfpN.c makes this earlier step unnecessary. The 2nd output file sample seen is 920515dr.lwh where we output
Bldg#     Dir     @ rng(mi)     bl'     bw'     bh'     area'^2     cw'     dwbl'     as seen as the first record header file. This is used for a Lidar city's .1 mi dr rings for the next program CDRCRD.c back at the Demolition Dust Model .

    NED's - National Elevation Data for Lidar