[pdal] Odd Inplace Reprojection Offsets

[Howard Butler]

Roger,

I think Paul submitted a patch to do this in https://github.com/PDAL/PDAL/issues/143 

Because of the way we use the Inplacereprojection filter from oracle, we load cloud objects per-file (with each having its own scale/offset), so we haven't run into this before.

Howard


On Jul 15, 2013, at 11:12 AM, Roger Bedell <sylvanascent at gmail.com> wrote:

> Hello Howard,
> Thanks for your quick reply. I would like to make a suggestion.
> 
> - if the offset is specified in the options, use it exactly as specified.
> -if not specified, derive a suitable offset from the data.
> 
> The reason is that with pgPointCloud, it is possible to add many LAS files to the same table. Only the user knows what might be a suitable offset beforehand.
> 
> I can try my hand at this if you dont have the time, if you can point me in the right direction.
> 
> Regards,
> Roger
> 
> Sent from my iPad
> 
> On Jul 15, 2013, at 16:00, Howard Butler <hobu.inc at gmail.com> wrote:
> 
>> 
>> On Jul 14, 2013, at 12:08 PM, Roger Bedell <rbedell at coordinatesolutions.com> wrote:
>>> Note, the Latitude is -21.4745 instead of around 36.09. 
>> 
>> The combination of scale * raw in32_t  + offset  values are overflowing past 32 bits. You probably need to dial the scale back to 0.0000001 or 0.000001.
>> 
>> 
>>>   <pc:scale>0.00000001000000</pc:scale>
>>>   <pc:offset>0.000297639019</pc:offset>
>>> Why is pc:offset = 0.0002976, when I explicitly asked for an offset of 33 in the pipeline XML?:
>> 
>> It's trying to calculate an offset for you based on the offset of the incoming data. Try turning up the debug noise with -d and -v 7 to see it spits out any illuminating noise.
>> 
>>> 
>>> <?xml version="1.0" encoding="utf-8"?>
>>> <Pipeline version="1.0">
>>> <Writer type="drivers.pgpointcloud.writer">
>>>   <Option name="connection">host='localhost' dbname='LIDAR2' user='postgres'</Option>
>>>   <Option name="table">qtest4326</Option>
>>>   <Option name="srid">4326</Option>
>>>   <Option name="overwrite">true</Option>
>>>   <Filter type="filters.chipper">
>>>     <Option name="capacity">400</Option>
>>>     <Filter type="filters.cache">
>>>       <Option name="max_cache_blocks">3</Option>
>>>       <Option name="cache_block_size">32184</Option>
>>>       <Filter type="filters.inplacereprojection">
>>>         <Option name="in_srs">EPSG:26914</Option>
>>>         <Option name="out_srs">EPSG:4326</Option>
>>>         <Option name="offset_x">-104</Option>
>>>         <Option name="offset_y">33</Option>
>>>         <Option name="offset_z">0.0</Option>
>>>         <Option name="scale_x">0.00000001</Option>
>>>         <Option name="scale_y">0.00000001</Option>
>>>         <Option name="scale_z">0.001</Option>
>>>         <Reader type="drivers.las.reader">
>>>           <Option name="filename">NRCSTEST10k.las</Option>
>>>           <Option name="spatialreference">EPSG:26914</Option>
>>>         </Reader>
>>>       </Filter>
>>>     </Filter>
>>>   </Filter>
>>> </Writer>
>>> </Pipeline>
>>> 
>>> When I increase the scale_x and scale_y to 0.0000001, it works better. To me this seems like the offset isn't being written correctly, and with the smaller scale, the 32 bit int overflows due to the incorrect offset_y value.
>> 
>> InplaceReprojection filter tries to be smart and allow you to merely set the scale value and back-calculate an offset based on the offset of the X dimension that's coming from the LAS reader. If that calculation overflows, bad things happen. I used to trap for overflows, but the cost of that trap happening per-dimension-per-point was prohibitive. Maybe I should try to put it back in debug mode or something.
>> 
>>> I was able to reproduce this using your UTM17.LAS sample data. Using this pipeline XML:
>>> 
>>> 
>>> <Pipeline version="1.0">
>>> <Writer type="drivers.las.writer">
>>>   <Option name="filename">UTM17_out.las</Option>
>>>   <Option name="spatialreference">EPSG:4326</Option>
>>>   <Filter type="filters.inplacereprojection">
>>>         <Option name="in_srs">EPSG:32617</Option>
>>>         <Option name="out_srs">EPSG:4326</Option>
>>>         <Option name="offset_x">-84.0</Option>
>>>         <Option name="offset_y">38.0</Option>
>>>         <Option name="offset_z">0.0</Option>
>>>         <Option name="scale_x">0.0000001</Option>
>>>         <Option name="scale_y">0.0000001</Option>
>>>         <Option name="scale_z">0.01</Option>
>>>         <Reader type="drivers.las.reader">
>>>           <Option name="filename">UTM17.las</Option>
>>>           <Option name="spatialreference">EPSG:32617</Option>
>>>         </Reader>
>>>       </Filter>
>>> </Writer>
>>> </Pipeline>
>> 
>> With current PDAL/PDAL:master, this seems to produce what is expected for me. Any chance you can try with the latest code? If you're still stuck, please catch me on IRC.
>> 
>> Thanks for being willing to dive into PDAL. It is still very much software that is growing, and there's still a number of pointy bits like this throughout.
>> 
>> Howard
>> 
>>