[Liblas-devel] Re: GetProj4() string and NAD83

[Howard Butler]

On Jan 10, 2012, at 8:50 AM, Adam Stylinski wrote:

> I am not a geologist so I'm sorry for my limited experience and knowledge when it comes to world coordinate systems and these APIs.  I'm probably doing something completely wrong, though obviously for only use cases that aren't WGS84.  

I tested with the IowaDNR-CloudPeakSoft-1.0-UTM15N.laz file, and I am getting the following when trying to reproject it to EPSG:4326:

> ERROR 1: point not within available datum shift grids
> error: Could not project point for ReprojectionTransform::point not within available datum shift grids0

It appears it is only this file that is causing the trouble.  I'm not quite sure what is going on though. I do have the grid files in place to support the datum transformation.  When I overrode the coordinate system from the one that the GeoTIFF keys define (with the NAVD88 datum) to a simple EPSG:26915, it reprojected fine for me.  Something must be up in either a) the file itself, or b) how proj is reprojecting the data or c) libLAS reprojection machinery.

Here's the code I used to test reprojecting. Instead of using proj directly, I am using the transforms mechanism of liblas::Reader. This allows me to control scaling and reproject the data as they are read.

I'll try to dig into this some more...

Howard



> #include <liblas/liblas.hpp>
> #include <string>
> #include <fstream>
> #include <iostream>
> 
> 
> int main()
> 
> {
>     liblas::ReaderFactory f;
> 
>     std::ifstream ifs("./test/data/IowaDNR-CloudPeakSoft-1.0-UTM15N.laz", std::ios::in|std::ios::binary);
>     liblas::Reader reader = f.CreateWithStream(ifs);
> 
>     liblas::Header const& header = reader.GetHeader();
>     liblas::SpatialReference const& in_ref = header.GetSRS();
>     
>     // liblas::SpatialReference in_ref;
>     // in_ref.SetFromUserInput("EPSG:26915");
>     std::cout << "Input SRS: " << in_ref.GetProj4() << std::endl;
>     
>     liblas::SpatialReference out_ref;
>     out_ref.SetFromUserInput("EPSG:4326");
>     
>     
>     // Copy the header's properties for the reprojection transform
>     // We will instead change the scale factors for the X, Y fields due 
>     // to the fact that the input file is in UTM and has a scale factor of 
>     // 0.001
>     liblas::Header reprojection_header(header);
>     reprojection_header.SetScale(0.0000001, 0.0000001, header.GetScaleZ());
>     reprojection_header.SetOffset(0, 0, header.GetOffsetZ());
> 
>     std::vector<liblas::TransformPtr> transforms;
> 
>     liblas::TransformPtr srs_transform = liblas::TransformPtr(new liblas::ReprojectionTransform(in_ref, out_ref, &reprojection_header));
>     transforms.push_back(srs_transform);
>     reader.SetTransforms(transforms);
>     
>     std::cout.precision(7); std::cout.setf(std::ios_base::fixed);
>     while (reader.ReadNextPoint())
>     {
>         liblas::Point const& p = reader.GetPoint();
>         std::cout << "x: " << p.GetX() << " y: " << p.GetY() << " z: " << p.GetZ() << std::endl;
>         break;
>     }
>     
>     return 0;
> }
>