[Liblas-devel] PROJ4 / Georeference info?

[Howard Butler]

This works almost exactly for me...

> >>> from liblas import file
> >>> input = file.File('TO_core_last_clip.las')
> >>> h = input.header
> >>> h.proj4
> ''
> >>> h.proj4 =  '+proj=utm +zone=17 +ellps=WGS84 +units=m'
> >>> h.proj4
> '+proj=utm +zone=17 +ellps=WGS84 +units=m'
> >>> output = file.File('junk.las',mode='w', header = h)
> >>> output.header.proj4
> '+proj=utm +zone=17 +ellps=WGS84 +units=m'
> >>> output.write(input.read(0))
> >>> output.close()
> >>> del output
> >>> output = file.File('junk.las')
> >>> 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 '
> >>>
>

The only thing I'm doing different is making sure to manually delete  
the reference to output when working in the interpreter to ensure that  
the file is actually flushed and not using 'header' as the name for my  
header that I'm writing to the file to ensure that I don't have a  
namespace collision or something.

Other than that, is it possible our versions are out of sync?


On Oct 30, 2008, at 2:25 PM, Jason Woolard wrote:

> hi again,
>
> Your example still returns an empty string for me:
>
> >>> input = file.File('C:\\liblas_test\\debug.las')
> >>> header = input.header
> >>> header.proj4
> ''
> >>> header.proj4 = '+proj=utm +zone=17 +ellps=WGS84 +units=m'
> >>> header.proj4
> '+proj=utm +zone=17 +ellps=WGS84 +units=m'
> >>> output = file.File('C:\\liblas_test\ 
> \junk_srs.las',mode='w',header=header)
> >>> output.header.proj4
> '+proj=utm +zone=17 +ellps=WGS84 +units=m'
> >>> output.write(input.read(0))
> >>> output.close()
> >>> output = file.File('C:\\liblas_test\\junk_srs.las',mode='r')
> >>> output.header.proj4
> ''
> >>>
>
> Am I maybe using Proj4 incorrectly? I'm just trying to add the  
> Projection/geoference info into the header where it doesn't exist.  
> Also, If I understood the posts from you and Mateusz it is necessary  
> to write something (even a zero record) if you are making any  
> changes to the input LAS file?
>
> Thanks again,
> Jason
>
>
>
>
> Howard Butler wrote:
>>
>> 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 '
>>>
>>
>>