Excellent. Send me the download link via private email and I'll put it up on http://liblas.org/samples
I have some ideas on how we could make this kind of a process much faster via GDAL's C API in las2las as well... Howard On Jul 23, 2010, at 11:56 AM, Aaron Reyna wrote: > 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! > > aaron > > On Fri, Jul 23, 2010 at 7:13 AM, Howard Butler <[email protected]> wrote: > Aaron, > > 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 > <http://liblas.org/samples>? We don't have any examples with coincident > lidar and imagery data, and this would be a great addition for things like > documentation, etc. > > Howard > > On Jul 22, 2010, at 8:01 PM, Aaron Reyna wrote: > > > Hi all, > > > > 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. > > > > 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. > > > > Warmly, > > > > aaron > > > > > > ############################################# > > import os, sys > > from liblas import file > > from liblas import header > > from liblas import color > > from liblas import * > > > > try: > > from osgeo import ogr, gdal > > from osgeo.gdalconst import * > > os.chdir('C:\\crap\\county\\Crook') > > except ImportError: > > import ogr, gdal > > from gdalconst import * > > os.chdir(r'C:\\crap\\county\\Crook') > > print "loaded" > > > > lassy = 'D:\\aaron_working\\DESCHUTTES\\Decluttered_no_birds\\DES_04593.las' > > > > # register all of the GDAL drivers > > gdal.AllRegister() > > > > # open the image > > img = gdal.Open('naip_1_1_1n_s_or013_2005_1_C4593.img', GA_ReadOnly) > > if img is None: > > print 'Could not open aster.img' > > sys.exit(1) > > print "loaded img" > > # get image size > > rows = img.RasterYSize > > cols = img.RasterXSize > > bands = img.RasterCount > > > > # get georeference info > > transform = img.GetGeoTransform() > > xOrigin = transform[0] > > yOrigin = transform[3] > > pixelWidth = transform[1] > > pixelHeight = transform[5] > > > > data=file.File(lassy, mode='r') > > > > print "creating .LAS file" > > h = header.Header() > > > > h.dataformat_id = 1 > > > > h.minor_version = 2 > > newdata=file.File('D:\\aaron_working\\DESCHUTTES\\Decluttered_no_birds\\DES_04593aaaab.las', > > mode='w', header=h) > > > > for p in data: > > pt = point.Point() > > xL = p.x > > yL = p.y > > pt.x = p.x > > pt.y = p.y > > pt.z = p.z > > pt.intensity = p.intensity > > pt.flightline_edge = p.flightline_edge > > pt.scan_flags = p.scan_flags > > pt.number_of_returns = p.number_of_returns > > pt.classification = p.classification > > pt.scan_angle = p.scan_angle > > pt.user_data = p.user_data > > > > # compute pixel offset > > xOffset = int((xL - xOrigin) / pixelWidth) > > yOffset = int((yL - yOrigin) / pixelHeight) > > # loop through the bands > > for j in range(bands): > > band1 = img.GetRasterBand(1) # 1-based index > > # read data and add the value to the string > > RED = band1.ReadAsArray(xOffset, yOffset, 1, 1) > > > > band2 = img.GetRasterBand(2) > > GREEN = band1.ReadAsArray(xOffset, yOffset, 1, 1) > > > > band3 = img.GetRasterBand(3) > > BLUE = band1.ReadAsArray(xOffset, yOffset, 1, 1) > > r16RED = int(RED)* 256 > > r16GREEN = int(GREEN) * 256 > > r16BLUE = int(BLUE) * 256 > > > > #print r16RED, r16GREEN, r16BLUE > > > > pt.color = color.Color(r16RED, r16GREEN, r16BLUE) > > newdata.write(pt) > > > > newdata.close() > > data.close() > > _______________________________________________ > > Liblas-devel mailing list > > [email protected] > > http://lists.osgeo.org/mailman/listinfo/liblas-devel > > _______________________________________________ Liblas-devel mailing list [email protected] http://lists.osgeo.org/mailman/listinfo/liblas-devel
