[Liblas-devel] Extract Color from NAIP with python, libLAS, and GDAL

[Howard Butler]

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 <hobu.inc at gmail.com> 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
> > Liblas-devel at lists.osgeo.org
> > http://lists.osgeo.org/mailman/listinfo/liblas-devel
> 
>