[Liblas-devel] PROJ4 / Georeference info?

[Howard Butler]

On Oct 29, 2008, at 10:27 PM, Jason Woolard wrote:

>
> Thanks for the great work here. The library is working quite well. I  
> have another quick question though. Sorry if I'm missing something  
> obvious.
>
> I'm trying to add projection and georeference info to LAS header  
> files that currently don't have any projection info. I was thinking  
> something like this would work:
>
> >>> from liblas import file
> >>> infile = liblas.file.File('C:\\liblas_test\\debug.las',mode='r')
> >>> h = infile.header
> >>> h.proj4 ='+proj=utm +zone=17 +ellps=WGS84 +units=m'
> >>> h.proj4
> '+proj=utm +zone=17 +ellps=WGS84 +units=m'
> >>> outfile = liblas.file.File('C:\\liblas_test\ 
> \debug_junk.las',mode='w', header=h)
> >>> outfile.header.proj4
> '+proj=utm +zone=17 +ellps=WGS84 +units=m'
> >>> outfile.close()
> >>> outfile = liblas.file.File('C:\\liblas_test\ 
> \debug_junk.las',mode='r')
> >>> outfile.header.proj4
> ''
> >>>



You have to write a point to the file or it is not valid.  I guess we  
should update the wiki to reflect this...  I also wonder if maybe we  
should have the writer try to throw an exception when nPoints is 0...

> >>> input = file.File('srs.las')
> >>> header = input.header
> >>> header.proj4
> '+proj=utm +zone=17 +ellps=WGS84 +units=m '

> >>> output = file.File('junk_srs.las',mode='w',header=header)
> >>> output.header.proj4
> '+proj=utm +zone=17 +ellps=WGS84 +units=m '
> >>> output.write(input.read(0)) # read the first point and write it  
> to the file
> >>> output.close()
> >>> output = file.File('junk_srs.las',mode='r')
> >>> output.header.proj4
> '+proj=tmerc +lat_0=0.000000000 +lon_0=-81.000000000 +k=0.999600  
> +x_0=500000.000 +y_0=0.000 +ellps=WGS84 +units=m '
>