[pdal] Odd Inplace Reprojection Offsets

[Roger Bedell]

Hello,
I am working on importing some NRCS Lidar data into pgPointCloud. One of
our requirements is to use 4326 to store the data, since the data covers
several UTM zones, and I want a single internal projection. So, using the
very convenient Inplace Reprojection Filter, I imported some data,
converting from EPSG:26914 to EPSG:4326, very similarly to the example in
the docs for the filter. Unfortunately, the Latitude (Y) data doesn't work
properly.

"{"pcid":1,"pt":[90,1,1,0,0,2,6,0,3,-2.29314e+007,0,0,0,-96.9688,-21.4745,264.218,204,0]}"

Note, the Latitude is -21.4745 instead of around 36.09.

I think that I have located the difficulty. If you look at the Schema, you
see this:

  <pc:dimension>
    <pc:position>15</pc:position>
    <pc:size>4</pc:size>
    <pc:description>Y coordinate as a long integer. You must use the scale
and offset information of the header to determine the double
value.</pc:description>
    <pc:name>Y</pc:name>
    <pc:interpretation>int32_t</pc:interpretation>
    <pc:scale>0.00000001000000</pc:scale>
    <pc:offset>0.000297639019</pc:offset>
    <pc:active>true</pc:active>
    <pc:uuid>0d1d8966-3e1c-4303-86de-69421f04003b</pc:uuid>
    <pc:parent_uuid>00000000-0000-0000-0000-000000000000</pc:parent_uuid>
  </pc:dimension>

Why is pc:offset = 0.0002976, when I explicitly asked for an offset of 33
in the pipeline XML?:

<?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.

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>

Then, running pcinfo on the result, you see this:

    <Metadata name="scale_x" type="double">1e-007</Metadata>
    <Metadata name="scale_y" type="double">1e-007</Metadata>
    <Metadata name="scale_z" type="double">0.01</Metadata>
    <Metadata name="offset_x" type="double">-85.48949644011782</Metadata>
    <Metadata name="offset_y" type="double">0.0003427359249286009</Metadata>
    <Metadata name="offset_z" type="double">0</Metadata>

One additional item, I requested offsets of -84.0 and 38.0 for the
UTM17.las sample, but instead, got -85.4894... and 0.0003427...

Shouldn't I get what I asked for?

Thanks,

Roger Bedell
Coordinate Solutions Inc.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/pdal/attachments/20130714/af6ab4f9/attachment.html>