C++ SpatialReference::getWKT方法代码示例

inline ptree toPTree(const SpatialReference& ref)
{
    ptree srs;

    srs.put("proj4", ref.getProj4());
    srs.put("prettywkt", ref.getWKT(SpatialReference::eHorizontalOnly, true));
    srs.put("wkt", ref.getWKT(SpatialReference::eHorizontalOnly, false));
    srs.put("compoundwkt", ref.getWKT(SpatialReference::eCompoundOK, false));
    srs.put("prettycompoundwkt", ref.getWKT(SpatialReference::eCompoundOK,
       true));

    return srs;
}

std::string transformWkt(std::string wkt, const SpatialReference& from,
    const SpatialReference& to)
{
    //ABELL - Should this throw?  Return empty string?
    if (from.empty() || to.empty())
        return wkt;

    gdal::SpatialRef fromRef(from.getWKT());
    gdal::SpatialRef toRef(to.getWKT());
    gdal::Geometry geom(wkt, fromRef);
    geom.transform(toRef);
    return geom.wkt();
}

void Stage::setSpatialReference(MetadataNode& m,
    const SpatialReference& spatialRef)
{
    m_spatialReference = spatialRef;

    auto pred = [](MetadataNode m){ return m.name() == "spatialreference"; };

    MetadataNode spatialNode = m.findChild(pred);
    if (spatialNode.empty())
    {
        m.add(spatialRef.toMetadata());
        m.add("spatialreference", spatialRef.getWKT(), "SRS of this stage");
        m.add("comp_spatialreference", spatialRef.getWKT(),
            "SRS of this stage");
    }
}

void Stage::setSpatialReference(MetadataNode& m,
    const SpatialReference& spatialRef)
{
    m_spatialReference = spatialRef;

    auto pred = [](MetadataNode m){ return m.name() == "spatialreference"; };

    MetadataNode spatialNode = m.findChild(pred);
    if (spatialNode.empty())
    {
        m.add(Utils::toMetadata(spatialRef));
        m.add("spatialreference",
           spatialRef.getWKT(SpatialReference::eHorizontalOnly, false),
           "SRS of this stage");
        m.add("comp_spatialreference",
            spatialRef.getWKT(SpatialReference::eCompoundOK, false),
            "SRS of this stage");
    }
}

bool SpatialReference::equals(const SpatialReference& input) const
{
    OGRSpatialReferenceH current =
        OSRNewSpatialReference(getWKT(eCompoundOK, false).c_str());
    OGRSpatialReferenceH other =
        OSRNewSpatialReference(input.getWKT(eCompoundOK, false).c_str());

    int output = OSRIsSame(current, other);
    OSRDestroySpatialReference(current);
    OSRDestroySpatialReference(other);

    return (output == 1);
}

Polygon Polygon::transform(const SpatialReference& ref) const
{
    if (m_srs.empty())
        throw pdal_error("Polygon::transform failed due to m_srs being empty");
    if (ref.empty())
        throw pdal_error("Polygon::transform failed due to ref being empty");

    gdal::SpatialRef fromRef(m_srs.getWKT());
    gdal::SpatialRef toRef(ref.getWKT());
    gdal::Geometry geom(wkt(12, true), fromRef);
    geom.transform(toRef);
    return Polygon(geom.wkt(), ref, m_ctx);
}

/// Add a Well-known Text VLR associated with the spatial reference.
/// \param  srs - Associated spatial reference.
/// \return  Whether the VLR was added.
bool LasWriter::addWktVlr(const SpatialReference& srs)
{
    std::string wkt = srs.getWKT(SpatialReference::eCompoundOK);
    if (wkt.empty())
        return false;

    std::vector<uint8_t> wktBytes(wkt.begin(), wkt.end());
    // This tacks a NULL to the end of the data, which is required by the spec.
    wktBytes.resize(wktBytes.size() + 1, 0);
    addVlr(TRANSFORM_USER_ID, WKT_RECORD_ID, "OGC Tranformation Record",
        wktBytes);

    // The data in the vector gets moved to the VLR, so we have to recreate it.
    std::vector<uint8_t> wktBytes2(wkt.begin(), wkt.end());
    wktBytes2.resize(wktBytes2.size() + 1, 0);
    addVlr(LIBLAS_USER_ID, WKT_RECORD_ID,
        "OGR variant of OpenGIS WKT SRS", wktBytes2);
    return true;
}

/// Set VLRs from the active spatial reference.
/// \param  srs - Active spatial reference.
void LasWriter::setVlrsFromSpatialRef(const SpatialReference& srs)
{
    VlrList vlrs;

#ifdef PDAL_HAVE_LIBGEOTIFF
    GeotiffSupport geotiff;
    geotiff.resetTags();

    std::string wkt = srs.getWKT(SpatialReference::eCompoundOK, false);
    geotiff.setWkt(wkt);

    addGeotiffVlr(geotiff, GEOTIFF_DIRECTORY_RECORD_ID,
        "GeoTiff GeoKeyDirectoryTag");
    addGeotiffVlr(geotiff, GEOTIFF_DOUBLES_RECORD_ID,
        "GeoTiff GeoDoubleParamsTag");
    addGeotiffVlr(geotiff, GEOTIFF_ASCII_RECORD_ID,
        "GeoTiff GeoAsciiParamsTag");
    addWktVlr(srs);
#endif // PDAL_HAVE_LIBGEOTIFF
}

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);
}

TEST(RialtoReaderTest, test)
{
    const std::string filename(Support::temppath("rialto4.gpkg"));
    
    FileUtils::deleteFile(filename);

    {
        // set up test database
        PointTable table;
        PointViewPtr inputView(new PointView(table));
        RialtoTest::Data* actualData = RialtoTest::sampleDataInit(table, inputView);
        RialtoTest::createDatabase(table, inputView, filename, 2);
    }

  RialtoReader reader;
  Options options;
  options.add("filename", filename);
  //options.add("verbose", LogLevel::Debug);
  reader.setOptions(options);

  {
      PointTable table;
      reader.prepare(table);
      
      const SpatialReference& srs = reader.getSpatialReference();
      const std::string& wkt = srs.getWKT();
      const SpatialReference srs4326("EPSG:4326");
      const std::string wkt4326 = srs4326.getWKT();
      EXPECT_EQ(wkt, wkt4326);
      
      const GpkgMatrixSet& info = reader.getMatrixSet();
      EXPECT_EQ(3u, info.getNumDimensions());
      const std::vector<GpkgDimension>& dims = info.getDimensions();
      EXPECT_EQ(3u, dims.size());
      EXPECT_EQ(89.0, dims[1].getMaximum());
  }
  
  {
      PointTable table;
      reader.prepare(table);

      PointViewSet viewSet = reader.execute(table);
      EXPECT_EQ(viewSet.size(), 1u);
      PointViewPtr view = *viewSet.begin();
      EXPECT_EQ(8u, view->size());
  }

  {
      BOX3D bounds(0.0, 0.0, 0.0, 10.0, 10.0, 10.0);
      options.remove("bounds");
      options.add("bounds", bounds);
      reader.setOptions(options);
  
      PointTable table;
      reader.prepare(table);

      PointViewSet viewSet = reader.execute(table);
      EXPECT_EQ(viewSet.size(), 1u);
      PointViewPtr view = *viewSet.begin();
      EXPECT_EQ(view->size(), 0u);
  }
  
  {
      BOX3D bounds(1.0, 1.0, -10000, 89.0, 89.0, 10000.0);
      options.remove("bounds");
      options.add("bounds", bounds);
      reader.setOptions(options);
  
      PointTable table;
      reader.prepare(table);

      PointViewSet viewSet = reader.execute(table);
      EXPECT_EQ(viewSet.size(), 1u);
      PointViewPtr view = *viewSet.begin();
      EXPECT_EQ(view->size(), 1u);
  }
}