LAS-NED.dia for newb.c Also see SP_LL.dia *.grd for NEW (UTM), or *.gnl PHL,PIT (State Plane). Clipping LAS overlaps: NEW las file is supposed to be 7.5' or 7.5'/60'/deg = 1/8 =.125 deg. So 8 x 8 = 64 las files could fit inside a 1 x 1 deg ned file. See topo.dia and pic samples link in the main html file. However, we have overlapping las files as pictured below. If a building is within the overlap [*], it would be counted twice. So, we want no overlap, except at 1 point on the border of .125 deg @ | in [|]. We also want to end at the 1 deg boundary of the NED's files. To be sure we assign the las to the correct 1 deg neds, if we average xg[0] & xg[1] of the LAS file (*.grd). Let las 1=xgmid=(xg[0]+xg[1])/2, then las 1 can be within a specific NED since the overlap isn't too large. Here we have west = (int) (xgmid+1) as the ned's W corner (using +). That is, all 89.?????? will fall between 90W -> 89W as seen in the 3rd diagram below. For a real case, we have xg[0],yg[0] = 90.126796, 29.936126 xg[1],yg[1] = 90.060663, 30.001376 xgmid=90.093729 ygmid=29.968751 So the integer values of xgmid,ygmid west= 91 north= 30. What topo quad out of the 8 might xgmid be? For each quad, we have .125 deg. ix=(west-xgmid)/.125 = .093729/.125 = 7.498 = 7 iy=(ygmid-north)/.125= .031249/.125 = .003 = 0 We use a variable quad[n] for n=0->8 multiple values of .125 or quad[0]=0 -> quad[8]=1. So quad[ix=7] beginning boundary = .875 and for quad[iy=0] = .000. The ending boundary for quad[7] = .875 + .125 = 1.000 (quad[8]). For quad[0] = .000 + .125 = .125 (quad[1]). Since we want to start and end exactly at the .125 boundaries to avoid overlaps, we skip into the *.las file from xg[0] to quad[7] as xg[0]-(west-quad[ix]). The ratio (xrat) of this to the full scale x of *.las is expressed as xrat=(xg[0]-(west-quad[ix]))/(xg[0]-xg[1]) = 90.126796-(91-.875)/(90.126796-90.060663) = (90.126796-90.125 )/.066133 = .001796/.066133 = .0271574. If the full x range is lx=1594 data points, then we start at ix1=lx*xrat = 1398 x .0271574 = 37.96 rounded up to 38. Using the same analogy for the ending point before lx @ xg[1], we have xrat=(xg[1]-(west-quad[ix+1]))/(xg[0]-xg[1]) = .917 x 1398=1282 (ix2) instead of using lx=1398. Recall west is exactly 91.000, the initial .125 quad. This might be depicted by the las2 in the diagram below. For beginning iy1, we would have yrat=(yg[1]-(north+quad[iy]))/(yg[1]-yg[0]) = .021 x 1594 = 34 and ending iy2, yrat=(yg[0]-(north-quad[iy+1]))/(yg[1]-yg[0]) = .937 x 1398 = 1493 rounded up. If the cell size 'cs' is 14.574 m, then for ix, we are skipping the 1st 34 data points and ending 1398-1282= 106 data points from lx, then 14.574 m x 34 points = 495m and ending 14.574 m x 106 data points= 1545m. xg[0] xg[1] [ [|] |] [ [|] |] [ las [*] las |] [ 1 [|] 2 |] [ [|] |] [ [|] |] [ [|] |] [ [|] |] xg[0] xg[1] | | | .875 1.000 ix1=38 ix2=1282 Linking to NED's: Also See NED.dia The UTM, state plane and geographical coordinates start at the equator and increase northward by latitude y. However, for the PC screen's pixel display with north at the top is pixel row yp = 0 so must be flipped. Examining the las file below, we have random paths of x,y's so we use iy,ix to calculate the file positions (filp) to place the building height data. The iy,ix's are from xmin,ymin, which is the SW corner. This means that the highest filp in the converted las (*.grd) file increases to the bottom of the file but represents further north. This means we have to flip images vertically for Google Earth displays, as seen on the right of ==>, since iy's are transferred to pixel rows yp. However, an alternate method would be allow filp to increase south like in NED's by y=(maxy-iy)*{} instead of y=iy*{} where {} = the cell size units. For NEW, the units are m as we then use UTM's to Lat/Lon shown below. This would allow no image flipping and make a link to NED's filp easier. N = ymax N N = ymax 1-------------2 iy=0 1-------------2 yp=0 1-------------2 iy=0 | | | | | | | | y=iy*{} | | | | | | | | | | y=(maxy-iy)*{} | las | | | | las | | *.grd | ==> | *.gif | | *.grd | E| *.pcx |W E| |W or E| *.pcx |W |^ | | flipped | |f | |f | | vertically | |i | |i | | | |l | | |l | | | |p | |p | | | |! | 0-------------3 iy=ly 0-------------3 yp=py 0-------------3 iy=ly xmin,ymin S S xmin,ymin S Now we zoom to the LAS building in *.grd into the NED's file: * = ixmin,ixmax,iyt,iyb building corners within the *.las file. Here we now use LAS' lx,ly as the # of data points in x,y. With the LAS's end UTM data points (lx,ly) and the building's extremities, we find the ratios xrat,yrat position to the lx,ly end points. We then apply these ratios to the corresponding xg[0],xg[1] ends points to find xgmid,ygmid (now used as the building's center). Even though we only use ix1 -> ix2, they are still counted in the ix=0->lx since lx is the x-range of the UTM data points. We could have ratioed ix1 & ix2 from the quad[ix] multiples of .125 deg also. 0 *.grd lx ym[1] ................. With the UTM zone = 15 for both sides :ixmin *iyt : of 90 deg, we can use ix,iy's min & max : * : to utm's then llutmf.exe in.utm -> ll.out. : xgmid,ygmid : ym[] = *.grd's vertical extent, y[] = : o : building's vertical extent iyb -> iyt. : iyb * : y[0]=ym[1] ->:y[0]* ixmax : -iyb*cs/fpm :...............: ym[0] west - NEDS - 91W N30W091.grd 90W :+----------------+ 30N If we let 'o' represent the geographic center of a an ideal :| las | rectangular building in a las file, @ coordinate xgmid,ygmid :| ix1__ix2 | derived above, then the # of deg into the NED's file=west-xgmid. :| | | | We keep the + since the NW NED's corner 'west' names the NED's :| | o | | file. Now since each NED's cell size 'cs'=9.25926e-5 deg or :| |____| | 1/3 arc sec ix=(west+xgmid)/9.25926e-5+6 is # of data points :| xg[0] xg[1] | along x to locate the building in the NED's file N30W091.grd. :| | The : represents the 6 pixels or data points overlap whose Lon :| | is west+6/cs to the left or about .000555... deg west & south. :| | This is the real left boundary in the NED's *.hdr files. :| | However, they just tried to calculate the Lon which we don't +------^----------+ 29N need as we now know it's 6 ix data points more. ................... ix=0 ix ix=10812 (xgmid-90)/cs lx data points