C++ PointViewPtr::spatialReference方法代码示例

TIndexKernel::FileInfo TIndexKernel::getFileInfo(KernelFactory& factory,
    const std::string& filename)
{
    FileInfo fileInfo;

    PipelineManager manager;
    manager.commonOptions() = m_manager.commonOptions();
    manager.stageOptions() = m_manager.stageOptions();

    // Need to make sure options get set.
    Stage& reader = manager.makeReader(filename, "");

    if (m_fastBoundary)
    {
        QuickInfo qi = reader.preview();

        std::stringstream polygon;
        polygon << "POLYGON ((";

        polygon <<         qi.m_bounds.minx << " " << qi.m_bounds.miny;
        polygon << ", " << qi.m_bounds.maxx << " " << qi.m_bounds.miny;
        polygon << ", " << qi.m_bounds.maxx << " " << qi.m_bounds.maxy;
        polygon << ", " << qi.m_bounds.minx << " " << qi.m_bounds.maxy;
        polygon << ", " << qi.m_bounds.minx << " " << qi.m_bounds.miny;
        polygon << "))";
        fileInfo.m_boundary = polygon.str();
        if (!qi.m_srs.empty())
            fileInfo.m_srs = qi.m_srs.getWKT();
    }
    else
    {
        Stage& hexer = manager.makeFilter("filters.hexbin", reader);

        PointTable table;
        hexer.prepare(table);
        PointViewSet set = hexer.execute(table);

        MetadataNode m = table.metadata();
        m = m.findChild("filters.hexbin:boundary");
        fileInfo.m_boundary = m.value();

        PointViewPtr v = *set.begin();
        if (!v->spatialReference().empty())
            fileInfo.m_srs = v->spatialReference().getWKT();
    }

    FileUtils::fileTimes(filename, &fileInfo.m_ctime, &fileInfo.m_mtime);
    fileInfo.m_filename = filename;

    return fileInfo;
}

PointViewSet MergeFilter::run(PointViewPtr in)
{
    PointViewSet viewSet;

    // If the SRS of all the point views aren't the same, print a warning
    // unless we're explicitly overriding the SRS.
    if (getSpatialReference().empty() &&
      (in->spatialReference() != m_view->spatialReference()))
        log()->get(LogLevel::Warning) << getName() << ": merging points "
            "with inconsistent spatial references." << std::endl;
    m_view->append(*in.get());
    viewSet.insert(m_view);
    return viewSet;
}

// Make sure that spatialreference works for random readers
TEST(json, issue_2159)
{
    class XReader : public Reader
    {
        std::string getName() const
            { return "readers.x"; }

        virtual void addDimensions(PointLayoutPtr layout)
        {
            using namespace Dimension;

            layout->registerDims( { Id::X, Id::Y, Id::Z } );
        }

        virtual point_count_t read(PointViewPtr v, point_count_t count)
        {
            using namespace Dimension;

            for (PointId idx = 0; idx < count; ++idx)
            {
                v->setField(Id::X, idx, idx);
                v->setField(Id::Y, idx, 10 * idx);
                v->setField(Id::Z, idx, 1.152);
            }
            return count;
        }
    };

    XReader xr;
    Options rOpts;
    rOpts.add("count", "1000");
    rOpts.add("spatialreference", "EPSG:4326");
    xr.setOptions(rOpts);

    StatsFilter f;
    f.setInput(xr);

    PointTable t;
    f.prepare(t);
    PointViewSet s = f.execute(t);
    PointViewPtr v = *(s.begin());
    SpatialReference srs = v->spatialReference();
    EXPECT_EQ(srs, SpatialReference("EPSG:4326"));
}

TEST(LasWriterTest, fix1063_1064_1065)
{
    std::string outfile = Support::temppath("out.las");
    std::string infile = Support::datapath("las/test1_4.las");

    FileUtils::deleteFile(outfile);


    std::string cmd = "pdal translate --writers.las.forward=all "
        "--writers.las.a_srs=\"EPSG:4326\" " + infile + " " + outfile;
    std::string output;
    Utils::run_shell_command(Support::binpath(cmd), output);

    Options o;
    o.add("filename", outfile);

    LasReader r;
    r.setOptions(o);

    PointTable t;
    r.prepare(t);
    PointViewSet s = r.execute(t);
    EXPECT_EQ(s.size(), 1u);
    PointViewPtr v = *s.begin();
    EXPECT_EQ(v->size(), 1000u);

    // https://github.com/PDAL/PDAL/issues/1063
    for (PointId idx = 0; idx < v->size(); ++idx)
        EXPECT_EQ(8, v->getFieldAs<int>(Dimension::Id::ClassFlags, idx));

    // https://github.com/PDAL/PDAL/issues/1064
    MetadataNode m = r.getMetadata();
    m = m.findChild("global_encoding");
    EXPECT_EQ(17, m.value<int>());

    // https://github.com/PDAL/PDAL/issues/1065
    SpatialReference ref = v->spatialReference();
    std::string wkt = "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4326\"]]";
    EXPECT_EQ(ref.getWKT(), wkt);
}

std::vector<PointId> SMRFilter::processGround(PointViewPtr view)
{
    log()->get(LogLevel::Info) << "processGround: Running SMRF...\n";

    // The algorithm consists of four conceptually distinct stages. The first is
    // the creation of the minimum surface (ZImin). The second is the processing
    // of the minimum surface, in which grid cells from the raster are
    // identified as either containing bare earth (BE) or objects (OBJ). This
    // second stage represents the heart of the algorithm. The third step is the
    // creation of a DEM from these gridded points. The fourth step is the
    // identification of the original LIDAR points as either BE or OBJ based on
    // their relationship to the interpolated

    std::vector<PointId> groundIdx;

    BOX2D bounds;
    view->calculateBounds(bounds);
    SpatialReference srs(view->spatialReference());

    // Determine the number of rows and columns at the given cell size.
    m_numCols = ((bounds.maxx - bounds.minx) / m_cellSize) + 1;
    m_numRows = ((bounds.maxy - bounds.miny) / m_cellSize) + 1;

    MatrixXd cx(m_numRows, m_numCols);
    MatrixXd cy(m_numRows, m_numCols);
    for (auto c = 0; c < m_numCols; ++c)
    {
        for (auto r = 0; r < m_numRows; ++r)
        {
            cx(r, c) = bounds.minx + (c + 0.5) * m_cellSize;
            cy(r, c) = bounds.miny + (r + 0.5) * m_cellSize;
        }
    }

    // STEP 1:

    // As with many other ground filtering algorithms, the first step is
    // generation of ZImin from the cell size parameter and the extent of the
    // data. The two vectors corresponding to [min:cellSize:max] for each
    // coordinate – xi and yi – may be supplied by the user or may be easily and
    // automatically calculated from the data. Without supplied ranges, the SMRF
    // algorithm creates a raster from the ceiling of the minimum to the floor
    // of the maximum values for each of the (x,y) dimensions. If the supplied
    // cell size parameter is not an integer, the same general rule applies to
    // values evenly divisible by the cell size. For example, if cell size is
    // equal to 0.5 m, and the x values range from 52345.6 to 52545.4, the range
    // would be [52346 52545].

    // The minimum surface grid ZImin defined by vectors (xi,yi) is filled with
    // the nearest, lowest elevation from the original point cloud (x,y,z)
    // values, provided that the distance to the nearest point does not exceed
    // the supplied cell size parameter. This provision means that some grid
    // points of ZImin will go unfilled. To fill these values, we rely on
    // computationally inexpensive image inpainting techniques. Image inpainting
    // involves the replacement of the empty cells in an image (or matrix) with
    // values calculated from other nearby values. It is a type of interpolation
    // technique derived from artistic replacement of damaged portions of
    // photographs and paintings, where preservation of texture is an important
    // concern (Bertalmio et al., 2000). When empty values are spread through
    // the image, and the ratio of filled to empty pixels is quite high, most
    // methods of inpainting will produce satisfactory results. In an evaluation
    // of inpainting methods on ground identification from the final terrain
    // model, we found that Laplacian techniques produced error rates nearly
    // three times higher than either an average of the eight nearest neighbors
    // or D’Errico’s spring-metaphor inpainting technique (D’Errico, 2004). The
    // spring-metaphor technique imagines springs connecting each cell with its
    // eight adjacent neighbors, where the inpainted value corresponds to the
    // lowest energy state of the set, and where the entire (sparse) set of
    // linear equations is solved using partial differential equations. Both of
    // these latter techniques were nearly the same with regards to total error,
    // with the spring technique performing slightly better than the k-nearest
    // neighbor (KNN) approach.

    MatrixXd ZImin = eigen::createMinMatrix(*view.get(), m_numRows, m_numCols,
                                            m_cellSize, bounds);

    // MatrixXd ZImin_painted = inpaintKnn(cx, cy, ZImin);
    // MatrixXd ZImin_painted = TPS(cx, cy, ZImin);
    MatrixXd ZImin_painted = expandingTPS(cx, cy, ZImin);

    if (!m_outDir.empty())
    {
        std::string filename = FileUtils::toAbsolutePath("zimin.tif", m_outDir);
        eigen::writeMatrix(ZImin, filename, "GTiff", m_cellSize, bounds, srs);

        filename = FileUtils::toAbsolutePath("zimin_painted.tif", m_outDir);
        eigen::writeMatrix(ZImin_painted, filename, "GTiff", m_cellSize, bounds, srs);
    }

    ZImin = ZImin_painted;

    // STEP 2:

    // The second stage of the ground identification algorithm involves the
    // application of a progressive morphological filter to the minimum surface
    // grid (ZImin). At the first iteration, the filter applies an image opening
    // operation to the minimum surface. An opening operation consists of an
    // application of an erosion filter followed by a dilation filter. The
    // erosion acts to snap relative high values to relative lows, where a
    // supplied window radius and shape (or structuring element) defines the
//.........这里部分代码省略.........

TIndexKernel::FileInfo TIndexKernel::getFileInfo(KernelFactory& factory,
    const std::string& filename)
{
    FileInfo fileInfo;

    StageFactory f;

    std::string driverName = f.inferReaderDriver(filename);
    Stage *s = f.createStage(driverName);
    Options ops;
    ops.add("filename", filename);
    setCommonOptions(ops);
    s->setOptions(ops);
    applyExtraStageOptionsRecursive(s);
    if (m_fastBoundary)
    {
        QuickInfo qi = s->preview();

        std::stringstream polygon;
        polygon << "POLYGON ((";

        polygon <<         qi.m_bounds.minx << " " << qi.m_bounds.miny;
        polygon << ", " << qi.m_bounds.maxx << " " << qi.m_bounds.miny;
        polygon << ", " << qi.m_bounds.maxx << " " << qi.m_bounds.maxy;
        polygon << ", " << qi.m_bounds.minx << " " << qi.m_bounds.maxy;
        polygon << ", " << qi.m_bounds.minx << " " << qi.m_bounds.miny;
        polygon << "))";
        fileInfo.m_boundary = polygon.str();
        if (!qi.m_srs.empty())
            fileInfo.m_srs = qi.m_srs.getWKT();
    }
    else
    {
        PointTable table;

        Stage *hexer = f.createStage("filters.hexbin");
        if (! hexer)
        {

            std::ostringstream oss;

            oss << "Unable to create hexer stage to create boundaries. "
                << "Is PDAL_DRIVER_PATH environment variable set?";
            throw pdal_error(oss.str());
        }
        hexer->setInput(*s);

        hexer->prepare(table);
        PointViewSet set = hexer->execute(table);

        MetadataNode m = table.metadata();
        m = m.findChild("filters.hexbin:boundary");
        fileInfo.m_boundary = m.value();

        PointViewPtr v = *set.begin();
        if (!v->spatialReference().empty())
            fileInfo.m_srs = v->spatialReference().getWKT();
    }

    FileUtils::fileTimes(filename, &fileInfo.m_ctime, &fileInfo.m_mtime);
    fileInfo.m_filename = filename;

    return fileInfo;
}