[pdal] Odd Inplace Reprojection Offsets

[Howard Butler]

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