NED's - National Elevation Data for Lidar

by     Buddy McCloskey         Contact

    1. The Need for NED's: 04/16/2013

As examined in PSI Wave => Buildings: , Lidar used to extract ground elevation with buildings and trees, must be filtered to determine ground levels. It was thought to examine the more established elevation methods by our National Elevation Data base.

            A. DEM - Digital Elevation Model:

As our first approach, the USGS said that DEM's are part of the NED's system, however, found out that it is not the same as the NED's format. Instead of expounding on a system we didn't use, you can see the details in DEM.FLD . Well after all the procedures, at least there was a SDTS2DEM.exe file to download that translate SDTS to DEM. However, that's after you unzip the GZ then the TAR with another downloaded software Free Opener. I had done two of these and failed to unzip others. Working with the two first, after finding the *.dem file, we could use another downloaded program 'DEM2XYZN.exe' that converts the *.dem to xyz triplets - which we need. However, in ASCII and in a profile, x is held constant and y is profiled - so we have a repeated x for different y & z making too large a file for just one tile.

                (1) DEMPXKML.c = DEM PcX pics & KML :

After reading 1DEM0897.pdf , we gave it a try. Re diagrammed in DEM.DIA for clarity, we see that instead of overlapping data points, they had variable # of y data points until you go deeper into the x's. After programming fun, we finally got a pic of one tile as seen in 5838n.gif . It's OK but needs to be flipped south to north. To find an x,y point in 5838n.grd then match with an LAS file for z elevation, not knowing if all tiles will convert and not knowing what the numbering tile location (n5838) naming convention without plotting it or extracting corners in the file, we abandoned DEM.

            B. NED to the Rescue:

Even though this format is a far cry from the original neds from the 90's, it's worth cracking the format because of the size reduction. First, to download the data, we refer to NED.FLD . NED has the unzipped file and the sub folder name (30062092) but that unknown name has a sub folder NED file of n30w091 - meaning Lat 30N, Lon 91W. This time, we know the location by the file name. In these folders, we also have the matching header file n30w090.hdr as these 7 parameters are shown and explained in NED.FLD again. That byte order has been known to shift from MSBFIRST to LSBFIRST. Therefore, we have to solve the *.flt 32bit error we seem to have earlier.

                (1) More Bites on Bits:

After reviewing ieee.dia again, we see that a character sti[1]=.= 20. To convert this, a 'C' command is itoa(sti[i],bits8,2); called Integer to Ascii where 2 is the base or radix converted of sti[i] or character [.] is shown as bits8. We would expect bits8 = 00010100 = 16 + 4 = 20 in binary. However, a few minor adjustments caused the problems. First of all, we only are shown 10100. Well we can just add the zero's to the left? No, to make matters worse, 10100 is really left justified and could be taken as 128 + 32 = 160 if ignoring this. We wouldn't have to bother with those 0's to the right? We expected 00010100 but we have 10100000 instead. So we first have to shift all bits 3 to the right then add the 3 zeros. So the length of bits8=10100 or 5 characters so we are 3 short of 8. 10100 ---> xxx10100 --> 00010100. We could and should have written our own itoa().

The other feature that snagged me is that Little Endian having the least significant digit first, opposite to our numbering system (Big Endian) was always to be applied? This would mean, that it applies to bits too that we should now reverse the bit pattern to bits8r=00101000? The answer is no. Only the order of the 4-Bytes are reversed. Finally, we have the corrected program IEEE.c = IEEE format showing all 3 methods of extraction. However, *.flt can now be both Big or Little Endian depending on it's byte order.

            C. NED.c = National Elevation Dataset :

Before we use these ground elevations, we want to make sure we extracted the elevations and coordinate corners correctly so we have our stand alone program NED.c as seen in NED.FLD again. We brought over subroutine ieee32() from IEEE.c as we convert the 4-bytes into 32 bits after we write out the 4-byte order denoted by the 'byte order'. The SW corner (our xg[0],yg[0]) is xllcorner,yllcorner where ll = lower left from the header file. The northwest corner is in the file name n30w090.flt and extracted as 30 deg N & -090 deg W for our xg[1],yg[1] corner. The data starts from the SW corner. The intervals of x,y (ix,iy) goes from 0 -> 10,812 cells whose size (cs) is given by degrees from within the header file n30w090.hdr.

Since we need x,y in meters, we convert the SW & NW corner to utm's by our internal subroutine llutm(), which was supplied by Chuck Gantz from Austin,TX University from his notes "Equations from USGS Bulletin 1532". We call the sw corner utm's as swx,swy. Then x in meters (xm)=swx+ix*csm and y in meters (ym)=swy+iy*csm if we convert cs -> csm with csm = meters @ lat of the NED's iy row. It was not used but just referenced as coordinates if needed. For a cross check, we also used geographical coordinates xgix=xg[0]-ix*arc3 and ygix=yg[0]+iy*arc3, where arc3 = cs, also the calculated cell size (1/3 arc sec Lat), then converted to utm by llutm() again. This method should be slightly more accurate than the former. Actually, we ran the difference between the methods for a fixed ix=1000 for iy=0-> 1000 in NED.b0 .

            D. NEDFGPK.c = NED Flt -> Grid Pcx & Kml :

After a grid file (*.grd) is made from NEDFGPK.c above, let's convert the smaller grid file to a picture that will finally be placed on Google Earth for a test match in the next program below. Originally, a *.pcx file, that must be converted to be displayed here. Also, it has a 10,812^2 pixels or a 117 mb pic file. However, here each 1 deg x 1 deg grid file is scaled to 14 colors by the zmin -> z max range. I take an 11" x 17" @ 400 dpi section of these files which is 61.9% of the full size to sell as T-shirts & place mats when I do arts & craft fairs under the name of Geo-Graphics . Not sure of the max display allowed here on this site, which also depends on your set screen resolution, I scaled down the pics that will appear in less (10%) detail below.

                (1) Sample NED Pics :

Hopefully, since we are working with our 1st standard data set, we should show New Orleans-West . For a section near my area, we have DelMarva-North , DelMarva-South , S. New Jersey , NW Chesapeake Bay and New York City-Long Island . Even the background of this web page that you see here in blue & white, is not clouds in a blue sky but also a crop of upstate PA NED terrain colors. We develop a colored height legend for NED's & buildings for Pittsburgh below.

    2. LAS Needs NED's:

With more examination of the LAS data coming in for our US major cities, the various coordinate systems and naming file conventions must have been determined by state levels. Also, the classification of ground and non ground data points vary also. At first, Miami & New Orleans, it seemed trivial because of low elevation and even for center city Philadelphia, it appears insignificant. However, expanding out to the 5.25 miles needed for the dust rings, the terrain changed rapidly. For Pittsburgh, with terrain near 1000', this is crucial. The tile size area coverage also varied as well as their units and data density.

            A. The Resolution Revolution :

What was also noticed in the LAS tile files, is that some had more data points than others. Obviously, less data appears in remote area which is good reasoning and concentrate on the populated building areas is what we need. However, when we tried to use the 'ground' data classification, we found many 'sink holes'. We use this term as missing data, not as the physical. There are 2 ways to fill these sink holes. We can fill them with adjacent data or shrink the ground until they are swallowed up. As an example, we used the 'ground' classification for Pittsburgh's southern suburbs and show this in PIT S Ground Tile .

Note the series of blue buildings. These points are classified other than ground, which we will pick up later. If we were to use the same center building data point as a ground elevation in this file, it is missing because it is in the blue building area. We had a routine that scanned along the perimeter of the building that should pick up some ground points and if none, we scanned +,- data points from building center out to 700' looking for ground. This takes too much CPU time and even then we missed some ground. We can see the scattered blue botches off the building which could also be trees, not ground. Off to the right side, we see the terrain rising by the red -> pink color transition.

                (1) LAS Sink Holes:

Before we shrink the number of sink holes, we should determine how many we have and should allow. Each of the Phl & Pit tile represents 10,000' x 10,000'. If each data point represented 1 ft, we should have 10^4 x 10^4 = 10^8 = 100 x 10^6 or 100 mega pixels (data points). Now the data points in one sample tile indicate that we have 3793795 data points or 3.793795 x 10^6. Note, this is only 3.79/100 = 3.79% data coverage. Also noting, that the data points are split up between ground and non ground, as seen in the corollary to PIT S Tile Ground S by PIT S Tile Building S , we probably only have near 2% buildings and even less when subtracting trees. Not matching these pics exactly, it may be difficult to collect all the building's data points, now red, replacing the blue. Off to the right in the higher terrain, we note the isolated pinkish-purple building clearly seen. These colored tiles indicate building height + terrain.

                (2) Closing Sink Holes:

                    (a) Scale Down:

As an example, let's consider we have a grid of 100 x 100 = 1.e4 or 10,000 possible points with say 4000 red and 6000 blue points. If we shrink the grid by 1/10, we would have 1,000 possible grid points with 400 reds & 600 blues. This means we would lose 4000-400=3600 red points. In our case, we have data & non data points where we would call this data 'red' and missing data 'blue'. To keep the original 4000 red points, we want to scale down from 10,000 points to 4000 points and block any over prints of blue pixels or points by a 1 line code: if (pix == blue) continue. Similar to this, we initialize the 4000 grid points with null's or 0's (blank or black screen) and any data fills the initialized data. We would call this a 1:1 ratio or 1 data point for each grid point. With equally spaced distribution, the entire grid would be filled. However, this not being the case, we have overprints of at least 2 red data points meaning 1 missing red data point from each just missed adjacent grid point.

                    (b) Now Scale Up?:

We closed up these sink holes so we can have fewer consecutive missing data points that can be filled. Well, roughly, if we have a cell size near 30', then a street width including both side sidewalks and grass areas of 60' should truly separate buildings across a street. After experimenting, we decided we could fill in 1 or 2 missing data points by interpolation between 2 known data points without compromising the data except on rare occasions. Using the PHL center city tile, we find tile ranges dxm=maxx-minx & dym=maxy-miny which are 10,000' each in PA's LAS. For New Orleans, we have a different range and units (UTM) so dxm is not equal dym. We define cell size cs=(dxm*dym)/npr in ft^2/point. For holes being filled, we first explain this approach for New Orleans in Cellsize.dia . With equally distributed points, which they aren't, lx=dxm/cs & ly=dym/cs, or the # of ix's & iy's horizontal & vertical cells for an assumed 1:1 ratio. Here we have for 1:1, lx=379 ly=379 cs=26.359' with dxm=10,000' & dym=10,000'. That is, different data points for every 26.359 ft if equally distributed and no overlapping points. Then if we allow 2 missing points for 1 data point, a new cell size now covers 3 data points or 9 x original # of grid points @ 1:1. So with a new cs=cs/3, we now have 3:1 lx=1138 ly=1138 cs=8.786' from lx=dxm/cs & ly=dym/cs rounded up.

                    (c) Pittsburgh Scale Examples:

Since most all this scaling, terrain and new filters are first applied to Pittsburgh in the fullest, here we have a summary of a few tile scale factors based on the data points per tile (npr) from the program PITgs.c .

npr=3187058
1:1 lx= 319 ly= 319 cs=31.377' dxm=10000 dym=10000
3:1 lx= 956 ly= 956 cs=10.459'

npr=3424025
1:1 lx= 342 ly= 342 cs=29.205' dxm=10000 dym=10000
3:1 lx=1027 ly=1027 cs= 9.735'

                    (d) Picture Examples:

To notice this effect visually, we examine 1:1 in Center City Phila #1 . Here the holes are filled and we have nice solid colors. Earlier, we mentioned we used sf = 10' for PHL which is seen in Center City Phila #2 along with it's Center City Phila #2 Zoomed . We still have the solid colors but sort of uneven sides and edges but hardly any holes (black dots). Notice this especially towards the big red building (Older Reading Terminal Train Station) NE of city hall (purple hollowed rectangle mid left). Finally, with the 3:1, we have the final at Center City Phila #3 . If you need to zoom up, you will see that even though some holes were still visible, it appears that more points are concentrated over the building areas and we can use NED's to fill in the missing ground points. This creates straighter building sides and corners for more building dimension accuracy.

Even a more striking example of this is the New Orleans data. Zoom up on your PC monitor to notice New Orleans Super Dome . Boosting the scale factor of the pic on left, we see the results from using the 3:1 ratio on the right as 3:1, lx=1398 ly=1594 cs= 14.574'. The smoother edges of the Super Dome and especially the highways make the 3:1 on the right more accurate and a more appealing pic despite more holes from less resolution than that of PHL's cs=8.786'? Notice that it appears that no 3 consecutive black dots are present to break the dome into more sections by horizontal separation. However, we may inadvertently do this due to increasing height of the dome by more filtering rules later. As mentioned previously, even if a building becomes broken into segments, the total dust is nearly the same but it's base rubble piles for rubbed dust acts independently from each other and not as a single larger cumulative pile with more surface area for rubbing, thus making this slightly less by only a few %.

            B. Linking LAS to NED's :

                (1) New Orleans - UTM's:

For buildings, we have to convert the building's grid files from the linear State Plane's or UTM's into xg,yg's for input to create NED's *.grd files. Also, we must keep track of the N -> S order of *.LAS to that of the NED's. We applied the boundary limits by UTM's for finding the appropriate NED's file and call the external utmllf.exe subroutine to convert to xg,yg for the NED's, pcx's & kml in NEWb.c . We see most of this explained in LAS_NED.dia as applied mostly to New Orleans.

The New Orleans tile sizes are larger and supposed to match the Topo Quad sizes. Here are some pics that would completely fill 1 Deg Lat/Lon NED blocks. At roughly the same Lat, we have Daytona Beach,FL . Up in my area of the NE, we have Bedminster,PA . A brief summary at E.PA & NJ . Notice the 8 x 8 quads/per degree but I show 2 deg x 1 deg (2 x 64 topo quads).

                (2) PHL, PIT ..... State Planes:

At New Orleans, we converted the linear UTM's to xg,yg by an external program UTMLLf.c . However, using the State Plane, we had external State Plane large files (PHL.sp, PIT.sp, etc...) that Corpscon 6.0.1 would convert after closing PHLb.c or PITb.c programs. Since this is not a command line program, it is done through Windows with input & output files. However, we now try to do this internally with subroutine SPtLL(). The first approximation is to do this linearly, that is no converging of the longitudes northward. The second step is to account for these converging longitudes working from the meridian through the center of the tiles to building corners. This is explained best with diagrams in SP_LL.DIA .

Since this should work with building plots, we can also use the mid building coordinate xgmid,ygmid to access the NED's elevation for that building's mid point by ixz=(west-xgmid)/csn + 6 + .5 where 'csn' is in deg as the NED's cell size (1/3 arc sec explained above). Also, iyz=(north-ygmid)/csn+.5 where .5 is rounded up. Note that the +6 data points is from the west (whole deg) but not the north (whole deg) in the file name. The file position is then filpg=iyz*10812+ixz where 10812 is the west -> east # of data points in all NED's files.

    3. LAS-NED's Full Application - Pittsburgh, PA:

With PIT having the most complete Lidar data to date, we will use these methods as a standard for all the other incoming cities and will redo PHL, NEW & MIA as well. It will appear best to show this as a flow chart depicted in PIT.fld .

            A. *.XML Types 1 & 2:

For the PHL Lidar, we used the PASDA web site listed earlier but they did not go as far west as PIT. The main USGS Lidar web site has PIT. Even though the same PIT region has the same Lidar data or *.LAS, they have different file names by two different coordinate classifications as well as different *.xml formats depending if we download from 1-Java or 2-Earth View. Originally using Java, it appeared easier to find specific tiles by Earth View. The only real difference is that Type 2 gives us the 4 corners and mid-geographical xg,yg coordinates of the tile whereas Type 1 gives us the 4 side boundaries xg,yg's. You can view these details in Type 1 42134.xml and Type 2 in 3846.xml .

            B. PITXKML?.c, ?=1,2:

We set up these tiles to be plotted in Google Earth by using the directory list of the *.xml files Type 1 & 2. Even though there are no tile pics to plot, we can set up to plot when they are available. However, the main advantage of doing this beforehand is making the *.kml file to plot now because Google Earth, on the 32 bit PC's, will plot their big red X (missing *.gif tile pics) that at least shows what tiles are downloaded around a major city ( PIT Grid #1 ). This way, using the GE scale legend, we can determine if all tiles are found with the 5.25 mi radius from center city and to show what tiles loaded are not needed outside this radius. A summary of these are displayed in Grid PIT #1 and Grid PIT #2 .

            C: PITgs.c: -> PITall & PITn29 Grid Files :

Notice, there are the corresponding 2 types of *.kml files and also 2 sub types pitall?.kml & pitn29?.kml. The PASDA & USGS regional PHL.las files are identical. This also means that the # of different classified categories with data points are the same in PHL & PIT. See the accumulation of data points in different classes of a sample tile in PIT Classifications . However, they still haven't classified 6 = Building. They did classify 2 = Ground & 9 = Water. This means that all others besides 2 & 9 are some kind of structures as either man made or trees. So this explains that Non 2 & 9 is n29 or structures. All the classifications includes 'all' types. The heights from 'all' these include the added ground or water elevations below these structures. The tile's scale factor for plotting was just mentioned above.

            D. PITgnl.c = PITgnl Grid Files :

As mentioned above, we can now have a new tile grid file composed of NED's ground elevations (pre *.gnl). However, since the LAS Lidar data is supposed to be more accurate, but short on coverage, we can then modify the NED-LAS Tile grid file points to lower values by the *.all & *.n29 grid files if applicable. The trick is to match the LAS tile with the NED's coordinates. Some insight is given here in LAS-NED.dia but applies mostly to New Orleans. Here, we use the subroutine SPtLL(), as explained above and in SP_LL.DIA again. However, we add another subroutine nedz() that converts tile coordinates and cell size 'cs' to NED files coordinates and cell size 'csn' = 1/3 Arc Sec = 9.25926e-5 deg. Then we have to adjust for the 6 pix overlaps in the south & west side of the appropriate 1 deg NED. The main core of this is the code:

ixz=(west[nn]+xgxl)/csn + 6 +.5 csn in deg for NED's xg[0]= + 6 pixel overlap to NED
iyz=(north[nn]-ygy)/csn+.5 csn in deg for NED's no north overlap
if(ixz >= cols) { ixz-=cols; ixz+=12; nn=1; } LAS to Eastern NED's

Here ixz,iyz are NED's coordinates and the xgx,ygx is the linear converted State Plane's x,y to Lat/Lon in SPtLL(). With the converging Lon's, we use 'xgxl' as the final xg @ x-Lon in SPtLL(). The nn=0 is the western NED's (81 deg), nn=1 the eastern NED's (80 deg). Notice the required 6 pixel overlap twice = 12 for ixz. Cols=10812 is the number of x-cells (csn) in NED's.

            E. PITgtp.c = PIT Tile Pics:

The first problem is the 6 pixel overlap of the NED's files just mentioned above. This happens right through Pittsburgh as seen in a section of PIT tile 42134 without correction. The + 6 pixel correction image is seen in PIT tile 42134 + 6 pix .

                    (1) NED's Color-Height Scale:

Now ready to proceed with NED's elevation, we need to reference a color-MSL height scale. For PIT, we could use LAS tiles zmin & zmax but sometimes they are in error. Also, if you want to view adjacent tile elevations, we have to know their Z max & mins also. After looking at the PIT region in Google Earth, we decided on zmin= 700' & zmax =1260'.

Since we also must consider the intervals between min -> max, we have another but shorter flow diagram in Z-Palette.fld . From zPAL.bug , we have the 14 height levels to match 14 colors. Next, we then match the corresponding 14 intervals of 256 color numbers for zn (nz). Finally, we interpolate between the 14 nz values for 256 color values (palette) and interpolate the last 3 columns (red,green,blue) values as they increase and decrease from 0 <--> 255.

This last 256 palette line section part becomes the *.pal processed by zPALhtm.c
and converted to a *.htm color code then inserted by cut & paste into PIT NED's Color-Ground MSL' Scale , an html file. If the height values aren't aligned in columns, use your mouse to lengthen or shorten the screen size on the right window boundary. If we run out of colors, we may cycle through the levels again continuing above z max. We see this in the PIT regional view of pure NED's in PIT Regional NED.

We now rewrite any lower values from the *.all tile grid to each matching grid points of the initialized *.gnl tile file from pure NED's. Then we apply the same procedures with the structures (*.n29). Actually, this should not occur but some showed up. Probably a low bush or grasslands with it's lower ground elevation + bush is still lower than the NED's ground elevation. In any event, the speckled or motley lower values superimposed on the NED's are seen in LAS-NED's PIT tile 42134gnl . Another example of a LAS modified NED's PIT regional tile is seen in LAS-NED's PIT tile 42133gnl . In Google Earth, we used a few PIT regional tile structures (*.n29) and did a screen capture for LAS-NED's PIT tiles n29. This is showing mostly building heights added onto terrain heights which still uses the PIT NED's Color-Ground MSL' Scale . A more regional zoom out is seen in PIT regional # 1.n29 and even more zoom out is seen in PIT regional #2.n29. Now we zoom back in with better screen capture details with PIT Structures #1, PIT Structures #2 and PIT Structures #3.

            F. PITb.c = PIT's Buildings:

This is the program that captures single buildings using structures (*.n29) and ground (*.gnl) files. A series of filters were explained above ( Single Buildings Extractions.DIA ) for PHL which still stands. However, with such terrain here, we decided it best to match every building's grid point (*.n29) with the corresponding elevation grid point in the *.gnl tile file and add a few more filters to eliminate trees.

                    (1) Two Building Examples:

Before applying filters, we want each building's grid point to be actual height above ground by subtracting off the matching grid point in the *.gnl tile file. The first example is the US Steel Building . Shown is a horizontal slice across the USS Building at iy=752 in the Tile File 42134.n29. We have the building + terrain height @ bh, then we subtract off the ground height (z0) from it's corresponding tile file 42134.gnl to leaving us with the building height above ground (bh0). The 'mpb' is the missing building points (2 allowed) that fail the filter (.6x856=513.6, 1.4x856=1198.4) where last bh0=856. Finally it is seen that at the last iy= 752 ixl= 329 ixr= 349 zl= 744 zr= 750, l=left & r=right. Note that 'zl' is the 1st z0 @ ix=329 and 'zr' is the last z0 @ ixr= 349 which are also written to our new dyxlrz.bh file. From this file, the next ix set is filter tested at iy=753.

The 2nd building is the Pit Dome . This uses the same procedure but across the apex of the dome. Notice the peaking (169' @ ix=444) of bh0 in the middle slice of the dome and the few missing data points filled in by interpolatio

                    (2) Pennsylvania Trees = Penn's Woods:

We can't forget that before the buildings, the trees dominate especially near 3 rivers and the associated hills. Treating the trees as buildings, before applying tree filters, we have the following GE screen captures tree strips along the river banks in PIT Tree # 0. A zoom out south of this pic, a progression of filters produced this sequence: PIT Tree # 1, PIT Tree # 2, PIT Tree # 3, and PIT Tree # 4.

In Google Earth, you can bring the little man down in the upper right corner for a street view. After placing him in the right location, exiting from street view, zooming in & out and rotating the pic, you can see the hill supporting the trees, being picked up as buildings, before more filters are applied. View PIT Tree Bank # 1, PIT Tree Bank # 2, PIT Tree Bank # 3 and finally PIT Tree Bank # 4.

                    (3) Buildings on Banks? Don't Bank on it:

Now we have to estimate the minds of construction engineers as how much ground will they plow to create a level base for buildings. Combining our *.gnl tile with the Google Earth's 3D buildings, we can see that buildings in PIT Building Bank # 1 are not located on steep hills - indicated by our PIT NED's Color-Ground MSL' Scale , where colors change rapidly (height gradient) depicts side of hills. However, as soon as a plateau is reached, 'let's build' as PIT Building Bank # 2 clearly indicates. Now zoom out to a regional view across the river in PIT Terrain-3D Buildings.

                    (4) Buildings or Trees?:

Now we present our tree-buildings filters in Building Slopes.dia . Let's apply these new filters to a lone building jammed into the side of a bank as depicted in Lone Bank Building .

tile=02,nb=0680,nf=14 ixmin= 90 ixmax= 161 iyt= 209 iyb= 209 bhc( 11879')/nc( 72) = abh(165.0')
Keep: nb= 680 ixmin= 90 ixmax= 161 iyt= 209 iyb= 209 abh=165.0 zmax( 795)-zmin( 732)= 63.0 dis= 748.2 slope=0.084
abh=165.0 < 235, dis= 748.2 < 1200, dz/2=63/2=31.5 < abh/4=41.25 and slope=0.084 < .1.

We test as a tree because the average (building?) height is less than 235'. Since shorter length than 1200', we can't reject it. Half the hill is less than 1/4th the building height and the ground slope is just under .1 so we must keep it. These filters helped eliminate the persistent red tree streaks, and most of the lower green & yellow tree streaks seen in Last Tree Streaks.

The choice of 1200' length also eliminates some elevated Interstate Highways. Imagine having a nearly 15 story building running for miles. This would add too much dust if we considered it a collapsed building. It's supporting columns are spaced and several feet thick, giving it much backup support against a PSI Wave. The columns stand and the concrete road can be 75' against the PSI Wave also and that's only across the width as a minimum, not it' length. There is no 15 story face for the PSI Wave as most of it passes underneath between the support columns.

                    (5) Final PIT Regional Building's Plot:

Finally, with most trees eliminated, we have our final PIT's Single Building Plot .

                    (6) Other Outputs from PITb.c :

The captured GE screen plots were from subroutines sptll(), pcx() and kml(). We usually just do 1 or 2 tiles to determine if sptll() coordinates are plotted in Google Earth correctly. The pcx() converts the *.gnl tile file grid into a *.pcx file while kml() makes a folder of all the buildings with their building corner coordinates and mid building coordinates for Google Earth plotting (PITb.kml). The *.pcx file is converted to *.gif for the kml() by Total Image Converter. All this was seen in PIT's Single Building Plot above.

The Google Earth *.kml files reference the *.gif pic file names as follows:
Alpha Characters         PHL Building .GIF ID
1 -> 2                             Tile # 1 -> 99
3 -> 6                             Building/Section # 1-9999
7 -> 8                             # of Floors 2 - 99

Our main output for the RADFO Model is the 19-byte building parameters for each building in PIT.bb by tally(). Since we are using a linear state plane, we have the scan of the building by ix,iy's. A sample of PIT.bb is seen in this 22 line screen capture PIT 19 Building Bytes.png . Note, there are some unprintable characters that show up as blank spaces. This is a binary file, and if opened in notepad as other text files are here, a different character set is seen. So we used an old DOS Edit program screen capture that displays most of the original characters. This new 19-byte PIT.bb file replaces xb11b. Each byte character represents 1-256. The 19 bytes are broken down as follows:

byte 1: Tile # 1-40
byte 2: Tile SW-x corner/1.e4
byte 3: Tile NW-y corner/1.e4
byte 4 & 5: Tile cell size * 1.e3 (cs3) in feet
byte 6 & 7: Building Sequence # 1-9999
byte 8 & 9: Building ix min
byte 10 & 11: Building ix max
byte 12 & 13: Building iy min
byte 14 & 15: Building iy max
byte 16 & 17: Building's # of cells (nc) for area
byte 18: Building height by # of floors (nf)
byte 19: Carriage Return and line feed \n

bytes 1, 6&7 & 18 are the tile #, building # and # of floors for the *.gif file name in pitb?.kml, ?=1,2.
Example: 01301524 = 01 3015 24, tile=1, nb=3015, nf=24.

            G. XB19B.c = eXtract pit Building's 19 Bytes:

After all the building parameters in PIT.bb produced above are extracted, direction and range from ground zero (main building in cityl.dat), are calculated by rotation of the building until the calculated area matches the real area [nc*sf)^2] of the building. We also have the cross width (cw) where we can then calculate dwbl=area/cw, the down PSI wave backup length supporting the cross width (cw) face. A simple plot of this is seen in Building PSI Orientation . The other one is Building - PSI Wave Orientation. A review of this complete explanation with diagrams is found in 2. Calculating Building Dimensions to Remold: A sample of the top section of the main big output file is seen in drlwh.PIT . The header line parameters are part of the file and read before the data.

            H. CDRCRD.c & PD%R.c = Super Duper Looper II:

At this stage on, all the cities' dust is derived from these last 2 programs, being the same for each city. This has been all explained in Super Duper Looper II earlier. For the final dust particle distribution (1 -> 449u) for each of the 53 dr=.1 mi concentric rings, we have ctydst1.PIT .

    4. Philadelphia, PA Redone:

Rather than go through and repeat all the procedure details again, we will just show the amendments and results. The new flow diagram, following the standard used for PIT above, is seen in PHL.fld . Our 1. PITXKML?.c, ?=1,2 is the same except we only use ?=1 (Java xml's) at the USGS Lidar Site. The 2. PHLgs.c is the same as PIT's explained above as well as 3. PHLgnl.c = PHLgnl Grid Files .

            A. 4. PHLgtp.c = PHL Tile Pics DEM:

The color palette has been modified to show more height details with the accompanying color - height scale for PHL terrain. See PHL GNL Color Height Scale for PHL Ground NED-Lidar . Notice the end of the PA tiles cross into NJ but stop just south & east of Camden. We may need that section as the distance to there is close to our 5.25 mi needed to capture all the possible collapsed building into the nuclear cloud. We return to this New Jersey south section after completing the rest of PHL. We used this PHL's GNL background to zoom in and insert some of Google Earth's 3D buildings as seen looking west in PHL Ground NED-Lidar & 3D Buildings .

For the PHL region, we have improved the details as seen in PHL Ground NED-Lidar & 3D Buildings PHL , then for downtown area, we zoom into PHL Ground NED-Lidar & 3D Buildings Zoom 1 , then more zoom in PHL Ground NED-Lidar & 3D Buildings Zoom 2 , and a final zoom into center city by PHL Ground NED-Lidar & 3D Buildings Zoom 3.

            B. 5. PHLb.c = PHL Building's:

The filters to eliminate trees in PIT was more effective there than in PHL because of less trees and ground slopes in PHL. It appears that the rejection of trees by the 10% 'slope' was about 3:1 (PIT:PHL) but the 'bank' rejection was over 5:1 in PIT compared to PHL. However, a surprise in PHL was that 3 tiles sent the number of buildings over 9999 which made the *.GIF name to push over the 8-character prefix allowed. We only used the 2 center city tiles for GE pics that both were < 10,000 buildings so no problem. However, a solution for that is for >= 10000 buildings, we switch to an A-Z alpha-character identification by the subroutine alphac(). A sample is seen in alphac() . We have a base number of 26 (the alphabet) instead of 10 digits. We add another 17,576 onto the 10,000 or a possible total 27,576 buildings.

            C. 6. PHLXB19B.c = eXtract Building's 19 Bytes:

Applying the PIT's PITXB19B.c seemed fine except one major difference. The state plane east value for PIT is less than PHL's 2690000E. Above, for PIT, the x (E) is byte 2. However, for PHL, we have 269 x 1.e4. Since we only have 256 characters, when > 256, it cycles back to zero. So for a 269, we would extract a 269-256=13. Rather than making a major change for the eastern byte to bytes, we just add the lost 256 here. Since the range of the eastern state plane near PHL will have no 13 x 1.e4, we have no conflict.

            D. 7. CDRCRD.c = City Delta Rings Crush-Rub Dust:

There is no modification of this program but we are now at a crossroad. For the input file 'drlwh.PHL', we now need to know if the buildings south east of Camden, mentioned and shown above, are to be included. This brings us to a whole new state tile file again.

    5. New Jersey South:

Again, we try to keep the standard from PIT-PHL programs and so we refer to another flow diagram again by NJ.fld . Finally finding that missing tile in PHL Ground NED-Lidar , we found one tile that covers this and more.

Before I run any .LAS file in my programs, I run a sure established program from LAS Tools by Martin Isenburg. See his summary of 000102.inf . Notice the range of X is 334404.81-316721.39 = 17684 while Y is 402763.27-379881.70 = 22881. PHL was 10,000 x 10,000 so (17684 x 22881)/1.e8 = 4.0463 e8/e8 = 4.0463 as much area covered as PHL tiles. With over 4 times the area, we should have about 4 times the # of points. However, we have 46,280,629 points which is at least 10 times the # of points of a PHL single tile. Another big surprise was the fact that all the points are classified as either Low Vegetation (3) or 3 times as much High Vegetation (5). This means no buildings?

            A. 2. NJgs.c = NJ Grid & Scale:

Before we do program # 1. NJSXKML.c , we better do program # 2. NJgs.c = NJ Grid & Scale before we know what type of pic tiles we may have. What we first noticed is the scale factor:

1:1 lx=2023 ly=2617 cs= 8.743' dxm=17683 dym=22882
3:1 lx=6068 ly=7852 cs= 2.914'

Since most all of PIT & PHL, cs was usually near 10'. To use 2.914' is an overkill. Using the 1:1 cs= 8.743' is fine. As for the data classification, it claims just all vegetation. Therefore, to have a non 3 or 5 is useless so we will have a type '*.all' with NJs. For NJSXKML.c , we only have 1 tile.

            B. 4. NJgtp.c = NJ Grid Tiles to Pcx:

We may want to see what image features (trees & shrubs) are available with the .all 3&5 classification before we apply # 3. NJgnl.c . To our surprise, we see the plot in NJ S. All depicting definite buildings in the NE corner and western mid section near the Delaware River. Also, a north-south changing elevated highway is clearly seen. The accompanying color-height table is given by NJ S. Building & Terrain Color . Notice, we dropped the scale east of the river. We could have made it the same as PHL across the river, but then fewer height details would show. We wanted to see trees & possible buildings. For PA, we had ground & water separation, for NJs, we have all veggies. Considering NJ as the Garden State, it was a good classification for most but the important note is that they didn't classify buildings just like their western state neighbor. Therefore, by pic deciphering, all buildings must be there as in PA's.

            C. 3. NJgnl.c = NJ Grid Ground by Ned & LAS:

In NED.FLD , we plotted terrain by quads using NEDGPKQ.c = NED Grid -> Pcx & Kml by Quads. Here, we found and plotted the appropriate NJ S. Quad 18 . Recalling with samples pics above, we have 8 x 8 topo quads/ 1 Deg x 1 Deq. Now it definitely seems as if NJs is using the mentioned quarter quads. We superimposed the NJ S. tile onto the NJ S. Quad 18 to depict the NJs 000102 tile as a quarter quad tile. See NJ S. Quarter Quad Tile .

For the NED-Lidar of the 000102 tile, we have NJ S. GNL Tile . The top NE section of this is shown as pure NED in NJ S. NED's cropped then adjusted for vegetation or land below the NED's terrain is clearly seen in NJ S. NED's-Lidar cropped . Notice the building footprints must be representing ground, not vegetation.

            D. 5. NJb.c = NJ Building's pcx, kml & bytes:

Here we have some major adjustments to make as we list them in order.

                (1) Linking 2 Different State Planes:

Well, what is common between NJ & PHL? Not state plane but Lat/Lon and UTM's. However, we need to work with linear distances first, as what we are used to. The first approach is to convert NJ to UTM's. Unfortunately, each lies in different UTM Zone's - west & east of 75 deg Lon. Since the origin of assumed ground zero is Liberty Places in PHL, we force NJ to this UTM zone by use of llutm4.c. From 000102.xml, , we copy & paste the NW corner boundary Lat=39.9375002304473, Long=-75.125000228232. Rather excessive precision here. The output is 000102.utm . We read the converted NJs to PHL UTM's as xm[],ym[]'s with xm[0],ym[1] represents the NW corner of the 000102 tile as shown njs-phl.dia . Then Liberty Place to the NW corner of NJs tile is converted to feet by nwx=(xm[0]-utme)*fpm and nwy=(ym[0]-utmn)*fpm where fpm = '/meter. So now instead of adhering to either state plane, we now have a Liberty Place - Tile 000102 Plane. Note that xmid=(swx+sex)/2 & ymid=(swy+nwy)/2 for this new plane so sptll() from PIT & PHL does not have to be modified.

However, with this large tile, we have # of buildings well over 10,000 and the alpha-character subroutine alphac() to handle 27,576 is also exceeded. Actual count was 32,316 building sections. Since we only have one tile, we can drop 1 digit from the 2-digit tile but add 1 digit to the building number and include the 2-digit # of floors (nf) to complete the 8-character *.gif file name for the GE NJb.kml file as shown below:

Alpha Characters         NJs Building .GIF ID
1 -> 1                             Tile # 1 -> 9
2 -> 6                             Building/Section # 1-99999
7 -> 8                             # of Floors 2 - 99

Trying to display all these buildings in Google Earth took awhile but viewed in NJs Buildings and a zoomed view in NJs Buildings Zoomed . This verifies that buildings were coded but not classified. However, because NJ had no significant slopes, many trees slipped by.

            E. 6. NJXB16B.c = NJ eXtract Building's 16 Bytes:

For the NJs directory, we only have 1 tile no we can drop the 1 tile character. Also, the NJs state plane is replaced by the PHL Liberty Place - Tile 000102 Plane, we can reconstruct that link here and eliminate SW-x & NW-y state plane corners. Since we now can use less bytes for tally() in NJb.c just explained above, we can now reuse the equations in njs-phl.dia again. This new 16 byte NJ.bb output file from NJXB16B.c replaces PHLXB19B.c . Each byte character represents 1-256. The 16 bytes are broken down as follows:

byte 1 & 2: Tile cell size * 1.e3 (cs3) in feet
byte 3 & 4: Building Sequence # 1-65535
byte 5 & 6: Building ix min
byte 7 & 8: Building ix max
byte 9 & 10: Building iy min
byte 11 & 12: Building iy max
byte 13 & 14: Building's # of cells (nc) for area
byte 15: Building height by # of floors (nf)
byte 16: Carriage Return and line feed \n

            F. 7. CDRCRD.c = City Delta Rings Crush-Rub Dust:

Because there were so many buildings, we decided to reject those buildings > 5.25 mi here rather than process excessive buildings in the next program PD%R.c . If we examine NJ.fld again, it is noted, we never use the programs below the dotted line so we never used this program. However, to show the closest buildings, we allowed a bit over 5.25 mi just to see how close they are in NJs drlwh under the rng(mi) column. Even though less than .6 mi from our intake limit, we didn't use any buildings from this tile.

Why did we bother with any of NJ south after seeing the final results? Well, we shortened up the *.bb file and can use this method for other cities. Another good reason is that we see how to handle NJ tiles and when we go north to capture NYC western suburbs like Newark, we hope we can apply these same programs.

            G. 7. CDRCRD.c for PHL:

Now we return back to PHL where we left off at D. 7. CDRCRD.c = City Delta Rings Crush-Rub Dust. We have no NJs drlwh with rng <= 5.25 mi but for PHL, we have extracted a few lines from drlwh.PHL .

            H. 8. PD%R.c = Super Duper Looper II:

Our final output file for PHL is ctydst1.PHL .

    6. San Francisco,CA:

Now we go for the gold rush out west. Again, from the USGS website for Lidar, we started with type 1-Java and later tried type 2-Earth View again. Finally realizing that type 2 can call up your city map by name and their map is better to outline any chosen area. If we download one tile at a time, it is immediate. So we have type 1 & 2. Both need to be plotted to make sure there are no overlaps as that would result in counting buildings twice. One problem discovered is that there are 3 different group of names identifying the lidar tiles. That is, including long directory names after they are unzipped. These different sample tile names are given below:

ARRA-CA_GoldenGate_2010_000779\ARRA-CA_GoldenGate_2010_000779.xml
ARRA-CA_SanFranCoast_2010_000112\ARRA-CA_SanFranCoast_2010_000112.xml
CA_SanFranBayFEMA_2004_000305\CA_SanFranBayFEMA_2004_000305.xml

After making a program to extract the file name's 6-digits without the long directories, we discovered that we had duplicates as seen below:

ARRA-CA_SanFranCoast_2010_000096\ARRA-CA_SanFranCoast_2010_000096.xml
CA_SanFranBayFEMA_2004_000096\CA_SanFranBayFEMA_2004_000096.xml

This would cause all sort of problems as their locations differ. These 3 agencies made no attempt to cross check their numbering system although they have 6 digits to work with. We resolved this conflict by checking the directory name 'Coast' and 'Bay' in the 2 directories and with the range in 'Coast", we would have 096.xml as a final name while 'Bay' would have 0096.xml. With this settled, we can outline the flow chart in SFO.fld . Here, it follows much the same pattern as PIT.

            A. SFOXKML?.c, ?=1,2:

We have a GE view of partial downloads in Type #1 and Type #2 . If we view Type 1 again, we have some duplicates and different sizes hiding under the larger tiles on the western coast. This is indicated by the whiter backgrounds with the overlap and the fainter red X underneath. This is another reason we went to type # 2 as most are of uniform sizes with no overlaps yet. Now we have Type 1 & 2 . It took quite a few days to unravel this mess. The duplicates 1 & 2 had to be resolved as well as the overlaps of the 3 different types mentioned above. Finally, with no duplicates or overlaps, Type 1 & 2 is all the SFO tiles needed. We can now generate grids & pics below.

            B: SFOgs.c & SFOgtp.c = SFO Grid Files, Scale & Pics: :

                (1) Another Resolution Revolution:

Here we discovered a shocker that set me back a few hours. After running LAS Tools by Martin Isenburg mentioned above, his distribution of classifications are seen in 0912-Homes and 0982-Center City . We first notice this gigantic increase in the number of point records: 7530385 for Homes and 9408030 for Center City. Also notice, that instead of this spread over an area of 10,000' x 10,000', it covers an area 1500m x 1500m or 4921' x 4921' or 24.2% of PHL's & PIT's tile areas, which ran near 3x10^6 points/tile at times. This could have result in only .726 x 10^6 points instead of 7.53 x 10^6 here. For SFO Center City, we had even more points @ 9408030 ~ 9.41x10^6 points.

What effect did this have on the scale factor? We see this in 0912-Homes 1:1 first. We wanted near 10' but notice at the 1:1, we have .98' ~ 1' instead of 10'. Now notice that we usually reduce the scale factor down by 1/3 to bring us near 10' in PIT & PHL which triples the size of the grid's & pic's lx & ly and would result in a 15,060 x 15,060 grid. We don't want to do that with SFO. Even with 1:1 here, this is still too much resolution already. For BAL, we had good resolution to use 1:1. Even with a large 5020 x 5020 grid, let's see what a pic would look like. This size pic would be rejected on the web to date. Using the 'all' classification, for the time being and using a crude color height scale, we notice a cropped section in 0912ac0 1:1 . If you look closely, you can see the outline of a few homes in darker blue. However, in another crop, we notice a strange appearance in 0912ac2 1:1 . Notice the change in point resolution near the lower left (SW corner). It appears as bands of more lidar points with a different flight route? Even there, we still have more holes than expected.

For the calculation of a scale factor explained above, we assumed that the points are somewhat equally distributed . The points should appear like a 'screen' with consistent weaving. It appears that we have many overlapping points that are not classified as such. Within the adjacent scan, we have new but repeated points but the gaps between still remain. Examining the bug file 0912-sfogs 1:1 , the 'idx' & 'idy' are the rounded off grid coordinates from the x & y differences (dxm & dxy) from tile origin. Notice that x & y only changed about <= .04m which results in the same idx(293) & idy(5019) for 12 consecutive points. Then 3 consecutive but then 10 more consecutive points. This is an overkill on a given point area but missed some gaps between adjacent points.

So instead of expanding to 1:3 as before, let's condense the gaps, shrink the sink holes, till they almost disappear. First, we use 3:1 instead of 1:3 which reduces the grid to 1673 x 1673 and increase the scale factor to almost 3' as seen. With a better color height scale, we have a new crop @ 0912ac3 3:1 . If you can zoom in, you will see many holes filled but a few buck shots remain. We go to the next level with a 5:1 @ 0912-Homes 5:1 with the associated crop pic @ 0912ac3 5:1 by raising the scale factor to 4.901'. Thus, it appears as a nice smooth plot with only a few scattered sink holes and we have reduced the grid size down to a more familiar 1004 x 1004 range and doubled the point resolution.

                (2) How Does Your Garden Grow?:

At the bottom of 0912-Homes and 0982-Center City again, we notice 4174050 Medium Vegetation (4) and 5966505 Medium Vegetation (4) points respectively. This is the majority of points, even over ground. After seeing the NJs classification for low & high veggies, I would bet that SFO did the same and hid the buildings in with the medium veggies. The first pic from SFOgtp.c show us all structures (non 2) in Homes-0912n2.gif . OK, now we go into the garden of veggies with Homes-0912veg2.gif . Yes, look at the homes that sprung up after being planted in the garden? I guess you plant the 'home' seed and within the 30-year mortgage, the house sprouts into a ripe full fledged dwelling. Now looking at Center City structures in Center City-0982n2.gif , could high rises also be grown within the garden at Center City-0982veg.gif ? What a great seed company but how long does it take to grow high rises and highways?

                (3) Let's Impose a Superimpose:

Noticing the buildings are in both n2 & veg, are there any differences between the two? Depending on what graphic software opens up your pics, you will really see the differences. For Windows 7+, if you save your pic to your computer, open up the one Home (n2), then minimize it. Now in a new tab, open up the (veg) of the matching Home and minimize that too. You will see 2 overlapping minimized icons. Click on the top one, line up your mouse with the underneath tab listing but don't click yet. With your eyes focusing on the 1st pic, click your mouse. You will see the change. Do it again and go back to the 1st pic. Repeat this and you have a minor animation. If a picture is worth 1000 words, how many words does animation produce with this added dimension?

Now, if while viewing the pic in the window, maximize the window with both pics and repeat this. For those that could not do this, the 'n2' pic shows hazy parallel images of strips running from NW to SE. These are more points and cross into the streets with higher elevations. Since these (n2's) are supposed to be structures > ground, why do they show up higher in the streets? Since the veggies are a better indicator of structures only than 'n2', we will use 'veg' instead to extract SFO building heights.

Even a better way that works on Windows XP and above, save the 2 pics in a unique directory you make up. Go through [Computer] to find that directory and the only 2 pics you just saved to it. Now on one of them, right click and [Open With] > [Windows Picture and Fax Viewer] only. This overrides any other graphics arts programs that may try to hog it all. With the blue |< & >| icons, switch back and forth as it should be scaled to full screen mode.

            C. SFOgnl.c = SFOgnl Grid Files :

Now we look at the matching extracted tile in the downloaded NED's 1 x 1 deg file for SFO tiles 0912ned.gif-raw and 0982 ned.gif-raw . Notice it is a little 'blocky' because we have such a wide range of elevation this time covering an area of 4921' x 4921' or 24.2% of 10,000' x 10,000' or zoomed up by a factor of 4.13. Since the tile x-range is 1500 m and the NED sf ~ 10m, then we only display 150 data points or pixels along a row or column. For a full screen display of at least an 800 x-pixels/monitor screen, that's only a valid data point for every 5.33 pixels, where the LAS tile has all 800 pixels displayed with valid data points, but not all ground.

                (1) Smooth the Blocks:

For the reasons mentioned above, instead of using the closest data point in the NED's *.grd files, we might as well interpolate the missing ground elevation points as we do for the building heights in tiles. For the SFO tile 0982, we see that we are using a data point for every 3.924'. For our matching .gnl tile from NED's, we are using roughly 33.7' S->N and 26.6' W->E at SFO's latitude for every data point. For nearly flat terrain or small hills, this is no problem. On the NED's converted *.grd file in SFOgnl.c , we add a double interpolation to the subroutine nedz(). The details are explained in NED Smooth 2.dia . Using SFOgtp.c again, we have new pure NED's tiles in 0912ned.gif and 0982ned.gif . Do a superimpose of the new with the original unsmoothed raw 0912 ned.gif-raw and 0982 ned.gif-raw to notice the change to the smoothed.

Now we adjust this pure NED's Center City 0982ned.gif to this same Lidar tile's ground, veggies or even low structures that are lower elevation than the NED's for a better ground base as seen in 0982gnl.gif . Note, if we use the superimpose method again, you will see the uniform building base footprint. Now I realize, they must have used similar methods as we did earlier. Comparing heights, when they jump because of a building, they must take the surrounding terrain by their lidar, average them across the base of the building because they can't penetrate the building to obtain the plowed level ground, which probably has a basement anyway.

                (2) Recycle Colors?:

With the topography varying rapidly, if we choose a range of 256 colors, we may loose details of the coastal buildings and features while trying to capture the higher elevations. If we choose lower height increments for detail, we reach the max height too quickly and don't show higher features. So, if we use the 256 colors again, we really have 512 colors but you have to have an eye with imagination. Does the green represent 100-200' or say 1100-1200'? Knowing the real terrain somewhat, it is possible. A rule must apply here - use linear increments. If we have 1-100' by 20's, and 800 to 1000 by 100', when we reach over 1000', then we would have too much detail for 1000-1100'. As an example, we cycled through the colors twice in 0912 Homes-2 Cycle All . This is just an 'all' structures & terrain LAS tile. Notice the terrain features are somewhat ragged because of the structures added as compared to 0912ned.gif just smoothed above with the 2-cycle colors.

Note that the new cycle starts from white to black to green. Knowing which colors are higher, then it becomes easy. Is the yellow-green crescent in the lower-left (SW) corner a valley or a ridge? Since red > yellow, it's a valley. Reverse the colors, it's a ridge. Notice the green as the peak in the whitish high ridge at the top left. It is a steep ridge as the colors cycle rapidly. The associated color height scale is given in SFO Color-Ground MSL' Scale .

Now for the 'all' for SFO, we have SFO all - 2 cycle. With not many features, we now bring out the structures (.veg) in SFO Veg 1 - 2 cycle, SFO Veg 2 - 2 cycle and SFO Veg 3 - 2 cycle. Use the magnifier or full screen view for details. For the smoothed .gnl (Ground - Neds - Lidar), we have SFO gnl 1 - 2 cycle and SFO gnl 2 - 2 cycle. The Google Earth 3D buildings with the .gnl base is SFO gnl 3 - 2 cycle. A feature showing the surface over the ocean is SFO Bay Waves seen faintly and SFO Bay Waves Expanded.

            D. SFOb.c = SFO Building's :

The only remaining adjustment is for the UTM's coordinate system. We used the approach we had used in NJb.c = NJ Building's above by making a relative ground 0 (G0) to building plane. All buildings > 5.25 mi will be eliminated here in order to carry over those coordinates within 2 bytes. This eliminates the swx4 & swx4 2-bytes defining the NW or SW tile corners and hope this allows us to use 17 instead of 19 bytes for the SFO.bb file to be extracted.

Using the distance x & y to G0 seems nice but can we standardize this? At first try, it looks feasible but we should set a standard now for all remaining and possible revamping of previously Lidar within 5.25 miles of G0 for the .bb file. A detailed examination & diagram for UTM's is explained in UTM Building Bytes.dia . Based on this revived method, we hope we will have a final 21 bytes extract file given below:

byte 1: Tile # 1-255
byte 2 & 3: Tile cell size * 1.e3 (cs3) in feet
byte 4 & 5: Building Sequence (nb) # 1-65535
byte 6 & 7: NW Tile Cor X m (xm[0]-500000)
byte 8 & 9: NW Tile Cor Y m (ym[1]-4170000)
byte 10 & 11: Building ix min
byte 12 & 13: Building ix max
byte 14 & 15: Building iy min
byte 16 & 17: Building iy max
byte 18 & 19: Building's # of cells (nc) for area
byte 20: Building height by # of floors (nf)
byte 21: Carriage Return with line feed \n

To verify the compression in meters, then uncompressing the binary character and using feet, is seen in the 2 program outputs SFOb.c & SFOxb21b.c
. The cells 8, 26, 64 to feet 156.0', 48.0 & 384.0' to rng = 4084.4'= 0.774 mi match even with that NW tile corner & G0 UTM reduction correction factors. This should be the standard for all remaining city LAS files and for even redoing MIA & NEW.

Google Earth made some updates and a rotating series of dots to indicate loading pics. Actually, with them appearing slowly anyway, I didn't have to be reminded. However, it slowed the displayed layers over double. It took 1 hr & 20 min to load the 2 SFO business district 2 tiles. Then, better be awake to make a screen capture soon or your thrown off line and lose the pic. Here is the 2 individual building tiles pic in SFO Buildings - 2 tiles. The greater details show the coordinate match and local UTM - Lat/Lon subroutine utmll() is valid.

            E. Final 3 CDRCRD.c, PD%R.c & CTY3dd.c :

                (1) Shifting Ground Zero Building:

While running CDRCRD.c , we discovered that we were missing type 1 tiles on the west coast of the peninsula. The Ground Zero highest building Transamerica Pyramid Zoomed In and Transamerica Pyramid Zoomed Out was in the NE section called the financial district. Also having a small peak, we shifted the G0 SW 1.4 miles as seen by the yellow line in Ground Zero Shift which now includes some of the Type 1 buildings. We also had 2 other reasons for this shift SW. The Ground Zero is now the 400' building, half the height but bigger area just south of the Civic Center seen in the front left of center of Down Town Hub . We will have this same dilemma when considering Manhattan Island, with the 2 building hub areas.

We finished off with PD%R.c & CTY3dd.c which produced the final SFO input file ctydst1.SFO . CTY3dd.c is reserved for later to show the finale between all cities.

    7. San Diego,CA: 05/12/2013

While we are on the west coast, it might be better to complete it with the available Lidar data. Now we have come to the most difficult Lidar data to date. By the first approach, although we could have most by capturing Type 1 data, not all goes right with different PC's with different internet connections. Once we attempt to download Type 1, if we miss any, we usually go to Type 2 where we can almost pin point tiles to fill in the missing tiles.

            A. SAN Type 1 XML:

Here we see the captured tiles for SAN Type #1 . Well at least we have more hits than misses. The scale in the lower left is 4.91 mi as we need 5.25 mi. You can see we overreached to the north quite a bit but we can let the programs sift out those > 5.25 mi. A sample Type 1 xml file is seen in SAN 280841.xml. Discovering that the x & y range of these tiles is 2048' by their State Plane is only .39 mi, we have our hands full.

            B. SAN Type 2 XML :

Now we attempted to capture the missing tiles by Type 2. This yielded SAN Type #2 . Notice, the scale was made the same to compare to Type 1. A sample of their Type 2 xml file is SAN 000680.xml. Notice the range is still hopefully the same @ (2048'), as seen by the same small red x sized box, but not all x's. Now notice, we have the bigger red X's - is this good or bad? Well, the borders don't coincide with the smaller red x's and they overlap. It's best to use the smaller and then if any areas are missing, use the larger.

Notice the convenience of smaller xml file size and 2 more corner readings, but we haven't been told the elevation units. However, experience shows it's the same as Type 1. These 2 different types have the same *.las file but different *.xml files describing it. However, the statement LinearUnits 'METER' puts a scare into it all. Any difference in x,y,z units will show up vividly when we plot the tile data pics with Type 1 & 2 together.

            C. SAN Type 1&2 XML:

Now we have the whole party convention seen in SAN Type 1&2 . As before, GE screen captures shows the whiter background and the brighter x's with overlaps. Notice in Coronado, it appears we have 2 larger tiles almost overlapping exactly. Actually, the brightness of these larger tiles hint that there are also a few that do. Well, the first task here is to eliminate the duplicate Types 1 & 2. We decided to scratch the Type 2 overlaps. This takes quite some time even with the GE expanded list of tiles names to find the duplicates, write them down, scratch them from the old directory list then run it again. The numbering convention does help a bit. As seen, we might have to include a few more tiles to the east - by Type 2.

            D. Corrected Composite SAN Type 1&2 XML:

After a few days, to say the least, holding at bay the Large X's, we have the smaller tiles in Corrected SAN Type 1&2 . We could say we threw out those that had counterfeit tickets to this event and restricted those by size who could enter. Also, with no holes barred, we invited a few more tile friends.

            E. CPU Crunch Time:

After a look at these massive number of tiles, it's time to grind the data. Well as always, here is a preview of the tasks in the flow diagram SAN.fld . Yes, at first glance, it appears more complex. We have already seen the results of SAN?XKML.c, ?=1,2 above by GE plots of the data without the pictures (Red X's).

            F: SANgs.c: -> SAN Grid Files & Scale: :

Well, here we have only All and High Veg (5) classes again. Since we have much smaller tiles, the # of point records should also be down but by how much?

                (1) Cell Size Scale Factor:

We had the case with other cities that we have to zoom up the scale by allowing 2 missing data points. We went from a 1:1 to 3:1. Then in SFO, we had to zoom down because of too many data points. Finally here, we resolve this issue by testing the 1:1 scale factor first. That is, if the cell size in feet (csf) < 4', too much detail so we zoom back to csf=6'. For example, @ 1:1, if csf=3', then 6/csf=2 so we have a 2:1. If the 1:1 scale factor (csf) > 15', we allow 2 missing points so the 1:1 goes to 3:1. The csf range is 4' < csf < 15'. Most 1:1's were 5.5' < csf < 12'. A few were > 15'.

            G. SANgnl.c = SANgnl Grid Files :

Only one NED's file was needed for SAN - the NW corner of -118W Lon, 33N Lat. We incorporated the double interpolation in the subroutine nedz() used and explained above in SFO, and probably redo all previous cities to see if there were major changes with their new building heights.

            H: SANgtp.c = SAN Grid Tiles to Pics :

To combine all the .all & .veg tiles, we would miss details so with adjacent north-south tiles of each in SAN 2723 .all and SAN 2723 .veg , we see again that the .veg tiles represent structures better than the .all and it eliminated most all the elevated highway running northeast-southwest. The complete composite of *.gnl tiles with less detail is shown in SAN gnl - 2 cycle color . These colors range from 0' to 432' seen by SAN Color-Ground MSL' Scale .

            I. SANb.c = SAN Building's:

                (1) Adjusted State Plane Coordinate System:

Well we could have been fine with using the State Plane but the zone info I had from a SAN website put SAN in Zone 6. However, using Corpscon 6.0.1 reveals nothing close to what I had in any of the tiles. This means an adjustment. We could convert Lat/Lon to UTM's but in case we cross UTM zones, which we don't here, let's apply a method that we could use anywhere that imitates a State Plane. Before the gruesome details indicated by State Plane Building Bytes.dia , examine the real pics SAN North to G0 and SAN East to G0 application of the first diagram in State Plane Building Bytes.dia .

The G0 is the One America Plaza on the left center showing the yellow, N-S line from the top of the tile it is located in. This is easier to see when the zoomed 2nd pic line runs to G0 from it's west side of it's tile. What might appear is that this building isn't high being just in the green range. However, note that it is in the 2nd of the 2-cycle color range being enclosed by the 1st cycle's blue -> cyan -> grey -> to the 2nd cycle's black -> green -> yellow -> orange peak.

After reading State Plane Building Bytes.dia , we see that the first diagram indicates that dxg0 = spe-xf[0] but since spe,spn we bad here, we used the GE tools to determine the dxg0,dyg0. However, any city LAS file that is State Plane, we could use this method to make sure the G0's state place (spe,spn) is nearly correct and then use dxg0,dyg0 = spe-xf[0],spn-yf[1] where xf[0],yf[1] is the state plane's NW corner of the tile given in the tile's .xml file. However, as explained, the range of the tile here was 2048'. More of a range may cause some larger tiles with building's near 5.25 mile, to be over the limit for the *.bb file - but they are eliminated before then with the range test anyway.

                (2) Building Adjustment & CC Plot:

As a test, we always want to make sure the NED's-.gnl ground level are being subtracted to correct for building heights. We see a sample run in SANb.bug . Then the real test is the GE plot of 2 adjacent West-East Tiles shown in SAN Center City Buildings . You can see by the shadows to the north of the building have the shape and length of it's corresponding building and height as well as some tiny images of some underneath sides by GE's real photo images.

                (3) Tally():

After seeing State Plane Building Bytes.dia and noting we have 820 tile Type 1's and almost 300 tile Type 2's, we need 2 bytes to define tile # since only 1 byte brings us to 256 tiles. We have our new 18 bytes for san.bb as follows:

byte 1 & 2: Tile # 1-65535
byte 3 & 4: Tile cell size * 1.e3 (cs3) in meters 1-65535
byte 5 & 6: Building Sequence (nb) # 1-65535

! Building Corners from xg0,yg0 1-65535 !
byte 7 & 8:       (x[0]-dxg0+.5 + tx0 + 30000)) ixmin' to G0
byte 9 & 10:     (x[2]-dxg0+.5 + tx0 + 30000)) ixmax' to G0
byte 11 & 12:   (y[0]-dyg0+.5 + ty0 + 30000) iyt or iymin' to G0
byte 13 & 14:   (y[1]-dyg0+.5 + ty0 + 30000)) iyb or iymax' to G0

byte 15 & 16:   Building's # of cells (nc) for area 1-65535
byte 17:           Building height by # of floors (nf) 1-256
byte 18:           Carriage Return with line feed \n

            J: SANxb18b.c = SAN eXtract Building 18 Bytes :

Now it appears we have reached the end of the line. It's time now to consider the large red X's back in SAN Type #2 . It's almost a whole new city so now we almost treat this as one below. Thinking we were almost done, let's quote that famous Chicago song, Only the Beginning.

    8. San Diego,CA - Part II - Scripps:

As always, but especially here, we need to see SAN.fld again to see what we need. At the very top, we show some *.bat files but never explained them until now.

            A. Saving, Unzipping and Renaming .LAS & .XML Files:

For both Types 1 & 2, we save each in their directory SAN1 or SAN2. After this, we unzip 1 file at a time. The unzipped file has a file number that is only a USGS download number. It chooses the SAN1 or SAN2 directory or we give it another directory name. However, it has quite a name of it's own subdirectories. With only a few tile files, we usually wipe out these subdirectories and rename the file manually to only 8 prefix characters, the old standard I must use in all the programs. Here, with 820 Type 1's and 300 Type 2's, it is quicker to make a special program as seen below:

                (1) SANbat.c :

This creates a batch file that runs from a command line. After making a basic & alphabetical (/b /on ) directory list of the *.zip files, along with /s for subdirectories in the command prompt, we input this and a sample output is seen in SAN 2 Batch File . Knowing some DOS commands, you may notice that I found the *.las file name in a subdirectory with the subdirectory as part of the file name. Eliminating that and copying only to the main directory (SAN2), we have a 6 prefix name. We could have used the 'move' command instead of 'xcopy', but one typo error and who knows where the files went after hours of downloading Type 2 files one at a time. Safer to overload the hard drive then delete the *.zip and subdirectories later.

                (2) Who got involved?:

Notice ChulaVista - it was a duplicate of Type 1's so eliminated. Most of the SanDiego_2005 were used except where it duplicated Type 1. Then at the bottom, we have Scripps continual month and year updates. Scripps Institution of Oceanography (sometimes referred to as SIO, Scripps Oceanography, or Scripps) in La Jolla, California, founded in 1903, is one of the oldest and largest centers for ocean and Earth science research, public service, undergraduate and graduate training in the world.

Well I guess a change of personnel over the years created the problems I encountered. I had some tiles with the same name, was 000012. Although they named the file by having the date as part of the file name, I had to change to A00012,B00012,C00012 instead. These had the same location so I finally chose the, best or latest date and used 000012. Others had the different numbers but the same location with the different dates. Still others were duplicate area dimensions as seen above, but didn't cover the same data points. Without going through the gruesome details on what and what not to choose, the list finally came down to 4 tiles.

Let's just use one for example as seen in SAN2 000006.xml. Now here, they use 'LinearUnits' METER and they mean it. That is x,y,z in meters.

            B. SAN2gss.c -> SAN2 Grid Files Scripps :

After we make adjustments to convert csm<-->csf as cell size for meters & feet, we have the matching bug file SAN2gss.bug . Notice that our scale factor rule 4' < csf < 15' is applied over and over as the scales run amuck. Also notice the initial wide range of dxm=5900m dym=3486m at times as 1609m = 1 mi. Another surprise is the fact that all are Unclassified(1) so we only use *.all.

            C. SANgnl.c & SANgtp.c :

A cropped .gnl showing building ground footprints is seen in Scripps' Cropped 00032g . Everything is standard as before here but if the tile files are not displayed as pics, we have no idea what actual area it covers and what are the details. The first 000006.las was the same area as 000012, 000028 & 000029 but had better features. We tried to see the full tile plots such as Scripps' Scaled SAN2 00029 . It was scaled to fit here because the bug file above showed 1980 x 2774 that is too large to display here.

What to notice is that there are shadows? When does Lidar detect shadows or with my scale, lower height values? This must mean some data tinkering or somebody dug trenches behind these buildings? However, when compared to 000006, it was also too short. A cropped portion of 00006 shows the missing section in Scripps' Scaled SAN2 000006 . This pic still has some shadows but shorter.

            D. Google Earth Hidden Rules? :

Too really complicate matters, when tying to plot all these Scripps' tiles, we were in total confusion. It appeared that we were plotting the wrong pics or the Lat\Lon's were wrong. Rechecking over and over, one glance at the right spot reveals all. Isolating 000006, we examine the Google Earth Screen capture in SAN2 GE Scripps' 000006 . The underlying isn't lying as it shows the continuity of building's extending from my pic. It is noticed that only the NE corner of the pic shows, which really puzzled me for a bit. To verify the overlap, we cropped the original 000006 as seen in SAN2 Scripps' 000006 Cropped to match the small piece that Google Earth displays, which represents about the lower 1/3 of the cropped pic.

I recalled this same problem with previous Google Earth plots. Whenever there is more black than pic, it does funny things. I corrected this by converting .gif to *.png and it worked. However, not here. One last try as I converted to *.jpg but that still failed. This means that there will be no GE screen capture of Scripps' tiles.

            E. Scripps, Use with Caution :

We thought we had duplicate coverage areas with 000008, 000030 & 000032. We first show Scripps' SAN2 000030 compared to the chosen Scripps' SAN2 000032 . This comparison clearly shows the missing data. All data means nothing unless seen graphically - a test for all Lidar?

            F. SANbS.c = SAN Buildings Scripps :

The big change from SANb.c is the fact that we have x,y,z in meters again. So, for the tally() to create the SANs.bb file, we have again another coordinate system. We could have appended SAN.bb but again, it's safer to keep them separate for now - made from 2 different types of .xml files. Such a shame, as we wanted to apply the new approach in State Plane Building Bytes.dia and finish this. We tried to retain as many features as possible from State Plane Building Bytes.dia to develop our new UTM Building Bytes.dia . If you bring them both up in 2 different editors, and flip between both, you will see that we still kept the basic methods but just changed the State Plane's feet to UTM meters.

Now we could use this method for any new or old tiles using UTM's. Two important notes here:

1. The G0 tile was State Plane, these Scripps' tiles are UTM so G0 won't be in any of these Scripps' tiles.
2. Make sure the UTM of these Scripps' tiles don't change UTM zones from the G0 building's UTM's.

Based on #1., we were thinking of applying to other cities but here, G0 is in the State Plane's tile. Therefore, the only real common grounds are Lat/Lon from the NW tile corner and G0. So, we convert the Lat/Lon to UTM's for nwxm0=484178.65, nwym0=3620000.73 by Corpscon 6.0.1. See bottom of UTM Building Bytes.dia for more. Since we now have utme,utmn from cityl.dat, then dxg0=utme-nwxm0, dyg0=utmn-nwym0 in meters, we are on our way. Now we have xm[],ym[] replacing xf[],yf[]. Also note that +- y-coordinate switching between screen vs. true coordinate systems. The good news is we can still use the same 18 building bytes in tally() but with a new file name 'SANs.bb'. This time, we have to verify a few buildings away from center city or G0, so we show some Scripps' Buildings .

            G. SANSCRXB.c = SAN SCRipps' eXtract Building - 18 byte :

The only difference between this and SANxb18b.c is mostly the name. Adding an 's' for Scripps would make it 9-char for the file name, but I can only work with 8. We were careful in the ft <--> m usage especially in the cell size scale factors 'csf' & 'csm' by just using 'cs'm & fpm conversion factor when needed. The other difference is that we output to append drlwh.SAN which means, this has to be run after SANxb18b.c .

            H. Final 3 CDRCRD.c, PD%R.c & CTY3dd.c :

These 3 are the same for any city. However, CDRCRD.c is ready to read drlwh.SAN. We notice in Scripps' Buildings and Corrected SAN Type 1&2 , we may have missed some buildings in Coronado which is too close to San Diego to ignore. Well now, we had an older version of extracting buildings but now, for new and missed buildings, we want them converted and placed into drlwh.SAN for these final 3 programs. This takes us into another routine and program below.

                (1) GETB.c=Google Earth Tools for Buildings :

Such a nice name that could also stand for GET Buildings. There are 2 possible ways to save GE buildings. One method saves ground elevations, the other doesn't. Since we only have a few, we chose the 'line' method with ground elevations. To see this whole routine, it was best that we make another HTML web page file so we can make other references to it in other web pages while not being stuck in this massive web page and making it even larger. We may add the 2nd routine after all later. You may come back here afterwards?

      Google Earth Tools for Capturing Buildings - Diagrams & Explanation

What is new is that the building dimensions are derived from the local geographical corner coordinate in Building 260 deg . We finished off with CDRCRD.c, PD%R.c & CTY3dd.c which produced the final SAN input file sandst1.SAN . CTY3dd.c is reserved for later to show the finale between all cities.

    9. Baltimore,MD:

Originally, this was started before San Francisco and San Diego. Once N. Korea made boasts of missiles reaching the West Coast of the US, I shifted emphasis there and searched for available city Lidar data just explained above.

As standard routine here, we view BAL.fld . Here we only have Type I since the county of Baltimore was the only available data. It was easy to capture all that was available, as most fell within 5.25 mi anyway. The numbering of the tiles is # North, # South, # East, # West from Center City. For 1S1E classification of points, we have:

848834 Unclassified (1)
992781 Ground (2)
732 Low Point (noise) (7)
1212126 Overlap Points (12)

A view of BAL 1s1e1.xml shows it is not in the standard .xml format and must be opened up as a .txt file. It shows State Plane & altitude units in feet.

            A. BALXKML.c=BAL Xml to KML tiles :

Here we only have .xml Type 1 so no need for a '1' in the file name.

            B: BALgs.c BAL Grid Files & Scale :

We produced the grid file *.all and *.n2 along with the *.sf. Data point densities were nearly the same and we boosted the scale factor cell size to 3. Most resulted in cell size (csf) from 5'-12' with a few 4's and a few 14'+.

            C. BALgnl.c = BALgnl Grid Files :

As mentioned before, BAL was completed before SFO & SAN so the "Smoothing of the Block" was done here in BAL by double interpolation in the subroutine nedz() first, but explained in SFO above.

            D: BALgtp.c = BAL Grid Tiles to Pics :

We still show that non ground (.n2) indicates structures better than .all by doing a 'superimpose', explained above, on tiles BAL 2s4e All with BAL 2s4e N2 . See more detailed structures in BAL 2s3e N2 , BAL 2s1w N2 and BAL 1s1e N2 which we reference a lot. The banding of some '.all' still streaks with some of '.n2' as clearly seen BAL 7s3e N2 .

We have composite BAL tiles of .gnl as seen in BAL gnl 1 cycle color and BAL gnl 2 cycle color . A zoom up of this last pic with GE's BAL 3D is BAL gnl 3D 2 cycle color . Zoom up for details. The composite .n2 for BAL is seen in BAL N2 1 cycle color where the banding is really vivid. However, BAL All 1 cycle color shows better without streaks. The 1 cycle color heights are given in BAL 1 cycle color heights for both ground & buildings.

            E. BALb.c = BAL Building's:

As a test, we always want to make sure the NED's-.gnl ground level are being subtracted to correct for building heights. We see a sample run in BALb.bug . In tally(), we used the 19 building bytes routine, the same as PHLb.c and PITb.c explained back @ 2.B(2) above for PHL and PIT, with the .bb table shown in 3.F.(5). This has been all supplanted by SAN's State Plane Building Bytes.dia explained above for future cities. However, the results for 2 adjacent tile building extractions are seen in BAL 2 Tile Buildings LL with Lat/Lon Grid and BAL 2 Tile Buildings St with street names.

            F. Final 3 CDRCRD.c, PD%R.c & CTY3dd.c :

We finished off with PD%R.c & CTY3dd.c which produced the final BAL input file ctydst1.BAL . CTY3dd.c is reserved for later to show the finale between all cities.

    10. Who Let the Bugs in?: 06/18/2013

As with many programs, notice that there are always bug fixes & updates. Well, we have to admit after all this, one flew into our NED's application. Actually, it isn't a bug but a misinterpretation of the NED's header data (*.hdr). Above, we explained about the 6-pixel overlap and got our best fit picture with a 12-pixel overlap when changing tiles from west to east. We applied these statements:

ixz=(west[nn]+xgxl)/csn + 6
iyz=(north[nn]-ygy)/csn

where csn is in deg for NED's. Note, no north overlap. Since we noted that all header files gave coordinates of the SW corner (ll), this suggested the overlap. True, since if we only overlap a side and bottom or top, we have overlapping of all NED's 1 deg x 1 deg adjacent blocks. However, after viewing PIT & SFO terrain results, I was suspicious. Within the NED's download, we unzipped many files, however, the only two I really needed were the *.flt & *.hdr files. Recalling that I had viewed the other files earlier, I had decided to keep a few for further reference if needed. The one file to note is Meta Data.htm . Even though it lacks some data, the bounding coordinates should have been in the *.hdr file also. Notice the findings below:

West_Bounding_Coordinate: -90.00055555556
East_Bounding_Coordinate: -88.99944444444
North_Bounding_Coordinate: 30.00055555556
South_Bounding_Coordinate: 28.99944444444

This means that all 4 sides overlap with their adjacent NED's file. The north extends 6 pixels further north than we thought and also further east. We never considered the initial 1 Deg NED as starting 6 pixel rows earlier than it's NW corner file name identifier.

            A. How Much Error?:

There is more error north-south than west-east because of the converging longitudes. Since 1/3 Arc Sec = 9.25926e-5 deg, the north-south would be 9.25926e-5 deg x 60 nm/deg x 1.15 mi/nm = 9.25926e-5 deg x 69 mi = 6.38809e-3 mi. Now x 5280'/mi = 33.7'. For 6 pixel data points, this amounts to 202.4 ft. For 40 deg N, this would be 155' from west to east. Since we now have good results from using our new double interpolation, being over 200' in distance off in hilly terrain could result in a 20' height error with a 10% grade. Consider coming off the ocean on the northern coast in SFO, we are picking 200' elevations further north along the bottom of the hill or into the ocean than we realized which calculates to higher building heights than actual.

            B. Correcting Overlap Formulas:

The best way to examine this overlap is to use a simple diagram as seen in NED's Overlaps . The diagram and explanation is good and clear. We shift NED blocks at the end of the whole 1 deg value. This way, we don't have to shift back into the previous NED 1 deg block with our double interpolation routine. When we shift from NED block nn=0 to nn=1 for the north-south, we now use

if(iyz0 >= rows-6) { nn=1; iyz0=(north[nn]-ygy)/csn + 6 }.

            C. ReRuns?:

Yes, with SFO and even SAN, we must apply this correction. Some use of the 2-cycle colors will be used as appropriate. Also, we will now use the double interpolation msl terrain height by the inverse distance squared method for PIT,PHL,NEW & MIA as well as the variable scale factor (cell size). We didn't use any of this or even any NED's for MIA nor the tree filters and new building filters for MIA & NEW. Most all the tile displays remain the same, except for those with the new variable scale factors.

The real changes will be in the individual buildings dust calculations. We will try to be as brief as possible by just showing some new pic's and the end results finale for all the 7 cities thus far. We think it would be best to work backwards from our previous completion demonstrating more radical changes progressing our way back to the original - MIA.

    11. Rerun - 1st 7 Cities:

        A. Baltimore, MD:

Our original tile file ( BAL gnl 1-cycle color) with the associated BAL Terrain-Building Color Height Scale , is replaced by our new terrain BAL gnl 2-cycle color with no new Color Height Scale which would be half the heights repeated color. Ex: 322' color now represents 161'. The 2nd color cycle would be 480/2 + half 322' again = 240' + 161 = 301'.

        B. San Diego, CA:

Our original mosaic tile file ( SAN gnl 2-cycle color - original), is replaced by our new terrain SAN gnl 2-cycle color with the associated SAN Terrain Color Height Scale . Except for the slight difference is size, you will hardly recognize the difference with the new. If you can do a superimpose, you will see the slight shift (6 pixels) north correction when compared to the GE base map underneath. For BAL & SAN, both cities are within only one 1 deg x 1 deg NED block. A close up view of one tile is seen in SAN gnl tile cropped .

        C. San Francisco, CA:

Here we can show the effect of the double interpolation terrain on a gnl tile comparing the SFO gnl Unsmoothed with the SFO gnl Smoothed - Terrain Only . Ignoring the elevation bases of the buildings, you can see the smoothing at the color-height changing contour boundaries. Another pair to compare is SFO gnl Unsmoothed 2 with SFO gnl Smoothed 2 - Terrain Only .

        D. Pittsburgh, PA:

The original mosaic tile file is PIT gnl 2-cycle color - original , is replaced by our new terrain PIT gnl 2-cycle color with the associated PIT Terrain Color Height Scale . Again, except for the difference is size, you will hardly see the difference with the new. If you can do a superimpose, you will see the slight shift (6 pixels) north correction when compared to the GE base map underneath. This is noticeable at the green blob island in the NW corner. Here, PIT is split west-east between two 1 deg x 1 deg NED blocks. Individual building plots of 2 center city tiles are seen in PIT Individual Buildings and PIT Individual Buildings South Zoomed .

        E. Philadelphia, PA:

Here, we will just show the new terrain plot in PHL gnl 1-cycle color with the associated PHL Terrain Color Height Scale . Like PIT, the terrain is split between two 1 deg x 1 deg adjacent blocks. However, unlike PIT, they are not split west-east but north-south of the 40 deg parallel.

        F. New Orleans, LA:

Again, we will just show the new terrain plot in NEW gnl 1-cycle color with the associated NEW Terrain Color Height Scale . Like PIT, the terrain is split between two 1 deg x 1 deg adjacent west-east NED blocks and like PHL, split between the north-south NED blocks. This means we have to open 4 NED's file in NEWgnl.c where we switch when we overlap 2 NED blocks.

For the building's mosaic tile display, we have NEW .all 1-cycle color with the associated NEW Building & Terrain Color Height Scale . For the individual building plots, we have NEW Individual Buildings and NEW Individual Buildings Zoomed East . Use of the variable scale factor is seen in the center city NEW Tile Buildings cc - original and NEW Tile Buildings cc . If we zoom up on the Super Dome in the lower right corner, we see more detail and less holes in the later. We compared .all to .gnl by superimpose for seeing terrain - structures by tiles NEW 469.all and NEW 469.gnl .

        G. Miami, FL:

Miami was handled like a new city never done before. We downloaded Type 2 tiles which are displayed in MIA Tiles . I had some *.las files from a few years ago but never considered I needed the *.xml files at that time. I made a special chart relating the LAS file name to the older FIU State Plane called MIA Tile-SP .

As the flow diagram MIA.fld indicates, we have the mosaics of MIA .all with the associated MIA Terrain - Building Color Height Scale and MIA .n29 . It is noted that the later (.n29) doesn't have as much green hazy between the buildings, which is similar to all the other cities' .n29's. Therefore, we will still use these as 'structures' to subtract off the .gnl file heights from, such as shown in the mosaic MIA .gnl . One change for the tree filters, due to lower shrubs and palm trees than elsewhere, we applied the tree filters for heights < 75'. We keep zooming in on Miami's .all files with MIA 1 .all and MIA 2 .all . We isolated the MIA Harbor Port .all and MIA Beach South .all . Notice the fine details using the variable scale factor (cell size) especially if we full screen zoom in on MIA Int Airport (001145.n29) .

    12. Finale - 1st 7 Cities:

These 7 cities are ranked by the Total Tons of dust from a 1-megaton blast. Now to see the particle distribution of each, first view the scale legend in Tons of Dust % Particle Distribution . Click on the 3-letter name for the chart, the total tons of dust for pic of the chart. The 1-449u denotes the X-axis while the y axis represents the radial distance from ground zero in .1 mi increments from 0 to 5.25 miles.

  Table         Total Tons
    PHL           2.183x10^7
    PIT           1.836x10^7
    BAL           1.507x10^7
    SFO           1.351x10^7
    MIA           9.382x10^6
    SAN           4.544x10^6
    NEW           2.199x10^6

    Major Update: Lidar -> RUC -> RADFO