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

[Howard Butler]

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