Don't think I can actually put up these examples, but I do have other examples that are in the public sector. Additionally, I believe I can toss the corresponding image up on my server, no problem. I'll hunt it down and throw it up there tonight!<div>
<br></div><div>aaron<br><br><div class="gmail_quote">On Fri, Jul 23, 2010 at 7:13 AM, Howard Butler <span dir="ltr"><<a href="mailto:hobu.inc@gmail.com">hobu.inc@gmail.com</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex;">
Aaron,<br>
<br>
Is there any chance I could download both of the files from somewhere and do a test with the latest development version? Additionally, is there any opportunity to add both of these files to the libLAS sample library <<a href="http://liblas.org/samples" target="_blank">http://liblas.org/samples</a>>? We don't have any examples with coincident lidar and imagery data, and this would be a great addition for things like documentation, etc.<br>
<br>
Howard<br>
<div><div></div><div class="h5"><br>
On Jul 22, 2010, at 8:01 PM, Aaron Reyna wrote:<br>
<br>
> Hi all,<br>
><br>
> 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.<br>
><br>
> 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.<br>
><br>
> Warmly,<br>
><br>
> aaron<br>
><br>
><br>
> #############################################<br>
> import os, sys<br>
> from liblas import file<br>
> from liblas import header<br>
> from liblas import color<br>
> from liblas import *<br>
><br>
> try:<br>
> from osgeo import ogr, gdal<br>
> from osgeo.gdalconst import *<br>
> os.chdir('C:\\crap\\county\\Crook')<br>
> except ImportError:<br>
> import ogr, gdal<br>
> from gdalconst import *<br>
> os.chdir(r'C:\\crap\\county\\Crook')<br>
> print "loaded"<br>
><br>
> lassy = 'D:\\aaron_working\\DESCHUTTES\\Decluttered_no_birds\\DES_04593.las'<br>
><br>
> # register all of the GDAL drivers<br>
> gdal.AllRegister()<br>
><br>
> # open the image<br>
> img = gdal.Open('naip_1_1_1n_s_or013_2005_1_C4593.img', GA_ReadOnly)<br>
> if img is None:<br>
> print 'Could not open aster.img'<br>
> sys.exit(1)<br>
> print "loaded img"<br>
> # get image size<br>
> rows = img.RasterYSize<br>
> cols = img.RasterXSize<br>
> bands = img.RasterCount<br>
><br>
> # get georeference info<br>
> transform = img.GetGeoTransform()<br>
> xOrigin = transform[0]<br>
> yOrigin = transform[3]<br>
> pixelWidth = transform[1]<br>
> pixelHeight = transform[5]<br>
><br>
> data=file.File(lassy, mode='r')<br>
><br>
> print "creating .LAS file"<br>
> h = header.Header()<br>
><br>
> h.dataformat_id = 1<br>
><br>
> h.minor_version = 2<br>
> newdata=file.File('D:\\aaron_working\\DESCHUTTES\\Decluttered_no_birds\\DES_04593aaaab.las', mode='w', header=h)<br>
><br>
> for p in data:<br>
> pt = point.Point()<br>
> xL = p.x<br>
> yL = p.y<br>
> pt.x = p.x<br>
> pt.y = p.y<br>
> pt.z = p.z<br>
> pt.intensity = p.intensity<br>
> pt.flightline_edge = p.flightline_edge<br>
> pt.scan_flags = p.scan_flags<br>
> pt.number_of_returns = p.number_of_returns<br>
> pt.classification = p.classification<br>
> pt.scan_angle = p.scan_angle<br>
> pt.user_data = p.user_data<br>
><br>
> # compute pixel offset<br>
> xOffset = int((xL - xOrigin) / pixelWidth)<br>
> yOffset = int((yL - yOrigin) / pixelHeight)<br>
> # loop through the bands<br>
> for j in range(bands):<br>
> band1 = img.GetRasterBand(1) # 1-based index<br>
> # read data and add the value to the string<br>
> RED = band1.ReadAsArray(xOffset, yOffset, 1, 1)<br>
><br>
> band2 = img.GetRasterBand(2)<br>
> GREEN = band1.ReadAsArray(xOffset, yOffset, 1, 1)<br>
><br>
> band3 = img.GetRasterBand(3)<br>
> BLUE = band1.ReadAsArray(xOffset, yOffset, 1, 1)<br>
> r16RED = int(RED)* 256<br>
> r16GREEN = int(GREEN) * 256<br>
> r16BLUE = int(BLUE) * 256<br>
><br>
> #print r16RED, r16GREEN, r16BLUE<br>
><br>
> pt.color = color.Color(r16RED, r16GREEN, r16BLUE)<br>
> newdata.write(pt)<br>
><br>
> newdata.close()<br>
> data.close()<br>
</div></div>> _______________________________________________<br>
> Liblas-devel mailing list<br>
> <a href="mailto:Liblas-devel@lists.osgeo.org">Liblas-devel@lists.osgeo.org</a><br>
> <a href="http://lists.osgeo.org/mailman/listinfo/liblas-devel" target="_blank">http://lists.osgeo.org/mailman/listinfo/liblas-devel</a><br>
<br>
</blockquote></div><br></div>