Hi all,<div><br></div><div>I put together this python script to extract color from a NAIP and apply it to the las file. I can pull out color values from the NAIP, and I scaled them up to 16 bit, but I seem to be having a problem writing that color value to file. Actually, everything will come through except the color, which seems weird. I thought I'd toss it up to the list to see if anyone had any idea's. Thanks in advance for any suggestions.</div>
<div><br></div><div>Also, keep in mind I have only been writing code about a year. I know, it's sloppy and slow... but that never stopped me before.</div><div><br></div><div>Warmly,</div><div><br></div><div>aaron <br>
<div><br></div><div><br></div><div>#############################################</div><div><div>import os, sys</div><div>from liblas import file</div><div>from liblas import header</div><div>from liblas import color</div>
<div>from liblas import *</div><div><br></div><div>try:</div><div> from osgeo import ogr, gdal</div><div> from osgeo.gdalconst import *</div><div> os.chdir('C:\\crap\\county\\Crook')</div><div>except ImportError:</div>
<div> import ogr, gdal</div><div> from gdalconst import *</div><div> os.chdir(r'C:\\crap\\county\\Crook')</div><div>print "loaded"</div><div><br></div><div>lassy = 'D:\\aaron_working\\DESCHUTTES\\Decluttered_no_birds\\DES_04593.las'</div>
<div><br></div><div># register all of the GDAL drivers</div><div>gdal.AllRegister()</div><div><br></div><div># open the image</div><div>img = gdal.Open('naip_1_1_1n_s_or013_2005_1_C4593.img', GA_ReadOnly)</div><div>
if img is None:</div><div> print 'Could not open aster.img'</div><div> sys.exit(1)</div><div>print "loaded img"</div><div># get image size</div><div>rows = img.RasterYSize</div><div>cols = img.RasterXSize</div>
<div>bands = img.RasterCount</div><div><br></div><div># get georeference info</div><div>transform = img.GetGeoTransform()</div><div>xOrigin = transform[0]</div><div>yOrigin = transform[3]</div><div>pixelWidth = transform[1]</div>
<div>pixelHeight = transform[5]</div><div><br></div><div>data=file.File(lassy, mode='r')</div><div><br></div><div>print "creating .LAS file"</div><div>h = header.Header()</div><div><br></div><div>h.dataformat_id = 1</div>
<div><br></div><div>h.minor_version = 2</div><div>newdata=file.File('D:\\aaron_working\\DESCHUTTES\\Decluttered_no_birds\\DES_04593aaaab.las', mode='w', header=h)</div><div><br></div><div>for p in data:</div>
<div> pt = point.Point() </div><div> xL = p.x</div><div> yL = p.y</div><div> pt.x = p.x</div><div> pt.y = p.y</div><div> pt.z = p.z</div><div> pt.intensity = p.intensity</div><div> pt.flightline_edge = p.flightline_edge</div>
<div> pt.scan_flags = p.scan_flags</div><div> pt.number_of_returns = p.number_of_returns</div><div> pt.classification = p.classification</div><div> pt.scan_angle = p.scan_angle</div><div> pt.user_data = p.user_data</div>
<div> </div><div> # compute pixel offset</div><div> xOffset = int((xL - xOrigin) / pixelWidth)</div><div> yOffset = int((yL - yOrigin) / pixelHeight)</div><div> # loop through the bands</div><div> for j in range(bands):</div>
<div> band1 = img.GetRasterBand(1) # 1-based index</div><div> # read data and add the value to the string</div><div> RED = band1.ReadAsArray(xOffset, yOffset, 1, 1)</div><div><br></div><div> band2 = img.GetRasterBand(2)</div>
<div> GREEN = band1.ReadAsArray(xOffset, yOffset, 1, 1)</div><div><br></div><div> band3 = img.GetRasterBand(3)</div><div> BLUE = band1.ReadAsArray(xOffset, yOffset, 1, 1)</div><div> r16RED = int(RED)* 256</div>
<div> r16GREEN = int(GREEN) * 256</div><div> r16BLUE = int(BLUE) * 256</div><div> </div><div> #print r16RED, r16GREEN, r16BLUE</div><div> </div><div> pt.color = color.Color(r16RED, r16GREEN, r16BLUE)</div>
<div> newdata.write(pt)</div><div> </div><div>newdata.close()</div><div>data.close()</div></div></div>