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;
> }
>
_______________________________________________
Liblas-devel mailing list
[email protected]
http://lists.osgeo.org/mailman/listinfo/liblas-devel
