C++ OGRGeometry::clone方法代码示例

OGRErr      OGRWarpedLayer::CreateFeature( OGRFeature *poFeature )
{
    OGRGeometry* poOrigGeom = poFeature->GetGeometryRef();
    OGRErr eErr;

    if( poOrigGeom != NULL )
    {
        if( m_poReversedCT == NULL )
            return OGRERR_FAILURE;

        OGRGeometry* poTransformedGeom = poOrigGeom->clone();
        if( poTransformedGeom->transform(m_poReversedCT) != OGRERR_NONE )
        {
            delete poTransformedGeom;
            return OGRERR_FAILURE;
        }
        poFeature->StealGeometry();
        poFeature->SetGeometryDirectly(poTransformedGeom);

        eErr = m_poDecoratedLayer->CreateFeature(poFeature);

        poFeature->StealGeometry();
        poFeature->SetGeometryDirectly(poOrigGeom);
        delete poTransformedGeom;
    }
    else
        eErr = m_poDecoratedLayer->CreateFeature(poFeature);

    return eErr;
}

OGRErr GTMWaypointLayer::CreateFeature (OGRFeature *poFeature)
{
    FILE* fp = poDS->getOutputFP();
    if (fp == NULL)
        return CE_Failure;

    OGRGeometry *poGeom = poFeature->GetGeometryRef();
    if ( poGeom == NULL )
    {
        CPLError( CE_Failure, CPLE_AppDefined,
                  "Features without geometry not supported by GTM writer in waypoints layer." );
        return OGRERR_FAILURE;
    }

    if (NULL != poCT)
    {
        poGeom = poGeom->clone();
        poGeom->transform( poCT );
    }


    switch( poGeom->getGeometryType() )
    {
    case wkbPoint:
    case wkbPoint25D:
    {
        OGRPoint* point = (OGRPoint*)poGeom;
        double lat = point->getY();
        double lon = point->getX();
        CheckAndFixCoordinatesValidity(lat, lon);
        poDS->checkBounds((float)lat, (float)lon);
        writeDouble(fp, lat);
        writeDouble(fp, lon);
        float altitude = 0.0;
        if (poGeom->getGeometryType() == wkbPoint25D)
            altitude = (float) point->getZ();

        WriteFeatureAttributes(poFeature, altitude);
        break;
    }

    default:
    {
        CPLError( CE_Failure, CPLE_NotSupported,
                  "Geometry type of `%s' not supported for 'waypoint' element.\n",
                  OGRGeometryTypeToName(poGeom->getGeometryType()) );
        return OGRERR_FAILURE;
    }
    }

    if (NULL != poCT)
        delete poGeom;

    return OGRERR_NONE;

}

/*!
  \brief Get feature (private)
  
  \return pointer to OGRFeature
  \return NULL not found
*/
OGRFeature *OGRVFKLayer::GetFeature(VFKFeature *poVFKFeature)
{
    OGRGeometry *poGeom;
    
    /* skip feature with unknown geometry type */
    if (poVFKFeature->GetGeometryType() == wkbUnknown)
	return NULL;
    
    /* get features geometry */
    poGeom = CreateGeometry(poVFKFeature);
    if (poGeom != NULL)
	poGeom->assignSpatialReference(poSRS);
    
    /* does it satisfy the spatial query, if there is one? */
    if (m_poFilterGeom != NULL && poGeom && !FilterGeometry(poGeom)) {
	return NULL;
    }
    
    /* convert the whole feature into an OGRFeature */
    OGRFeature *poOGRFeature = new OGRFeature(GetLayerDefn());
    poOGRFeature->SetFID(poVFKFeature->GetFID());
    // poOGRFeature->SetFID(++m_iNextFeature);
    
    for (int iField = 0; iField < poDataBlock->GetPropertyCount(); iField++) {
	if (poVFKFeature->GetProperty(iField)->IsNull())
	    continue;
	OGRFieldType fType = poOGRFeature->GetDefnRef()->GetFieldDefn(iField)->GetType();
	if (fType == OFTInteger) 
	    poOGRFeature->SetField(iField,
				   poVFKFeature->GetProperty(iField)->GetValueI());
	else if (fType == OFTReal)
	    poOGRFeature->SetField(iField,
				   poVFKFeature->GetProperty(iField)->GetValueD());
	else
	    poOGRFeature->SetField(iField,
				   poVFKFeature->GetProperty(iField)->GetValueS());
    }
    
    /* test against the attribute query */
    if (m_poAttrQuery != NULL &&
	!m_poAttrQuery->Evaluate(poOGRFeature)) {
	delete poOGRFeature;
	return NULL;
    }
    
    if (poGeom)
	poOGRFeature->SetGeometryDirectly(poGeom->clone());
    
    return poOGRFeature;
}

/*!
  \brief Get feature (private)

  \return pointer to OGRFeature or NULL not found
*/
OGRFeature *OGRVFKLayer::GetFeature(IVFKFeature *poVFKFeature)
{
    OGRGeometry *poGeom;

    /* skip feature with unknown geometry type */
    if (poVFKFeature->GetGeometryType() == wkbUnknown)
        return NULL;

    /* get features geometry */
    poGeom = CreateGeometry(poVFKFeature);
    if (poGeom != NULL)
        poGeom->assignSpatialReference(poSRS);

    /* does it satisfy the spatial query, if there is one? */
    if (m_poFilterGeom != NULL && poGeom && !FilterGeometry(poGeom)) {
        return NULL;
    }

    /* convert the whole feature into an OGRFeature */
    OGRFeature *poOGRFeature = new OGRFeature(GetLayerDefn());
    poOGRFeature->SetFID(poVFKFeature->GetFID());
    // poOGRFeature->SetFID(++m_iNextFeature);

    poVFKFeature->LoadProperties(poOGRFeature);

    /* test against the attribute query */
    if (m_poAttrQuery != NULL &&
            !m_poAttrQuery->Evaluate(poOGRFeature)) {
        delete poOGRFeature;
        return NULL;
    }

    if (poGeom)
        poOGRFeature->SetGeometryDirectly(poGeom->clone());

    return poOGRFeature;
}

OGRErr GTMTrackLayer::ICreateFeature (OGRFeature *poFeature)
{
    VSILFILE* fpTmpTrackpoints = poDS->getTmpTrackpointsFP();
    if (fpTmpTrackpoints == nullptr)
        return OGRERR_FAILURE;

    VSILFILE* fpTmpTracks = poDS->getTmpTracksFP();
    if (fpTmpTracks == nullptr)
        return OGRERR_FAILURE;

    OGRGeometry *poGeom = poFeature->GetGeometryRef();
    if ( poGeom == nullptr )
    {
        CPLError( CE_Failure, CPLE_AppDefined,
                  "Features without geometry not supported by GTM writer in "
                  "track layer." );
        return OGRERR_FAILURE;
    }

    if (nullptr != poCT)
    {
        poGeom = poGeom->clone();
        poGeom->transform( poCT );
    }

    switch( poGeom->getGeometryType() )
    {
    case wkbLineString:
    case wkbLineString25D:
    {
        WriteFeatureAttributes(poFeature);
        OGRLineString* line = poGeom->toLineString();
        for(int i = 0; i < line->getNumPoints(); ++i)
        {
            double lat = line->getY(i);
            double lon = line->getX(i);
            float altitude = 0;
            CheckAndFixCoordinatesValidity(lat, lon);
            poDS->checkBounds((float)lat, (float)lon);
            if (line->getGeometryType() == wkbLineString25D)
                altitude = static_cast<float>(line->getZ(i));
            WriteTrackpoint( lat, lon, altitude, i==0 );
        }
        break;
    }

    case wkbMultiLineString:
    case wkbMultiLineString25D:
    {
        for( auto&& line: poGeom->toMultiLineString() )
        {
            WriteFeatureAttributes(poFeature);
            int n = line->getNumPoints();
            for(int i = 0; i < n; ++i)
            {
                double lat = line->getY(i);
                double lon = line->getX(i);
                float altitude = 0;
                CheckAndFixCoordinatesValidity(lat, lon);
                if (line->getGeometryType() == wkbLineString25D)
                    altitude = static_cast<float>(line->getZ(i));
                WriteTrackpoint( lat, lon, altitude, i==0 );
            }
        }
        break;
    }

    default:
    {
        CPLError( CE_Failure, CPLE_NotSupported,
                  "Geometry type of `%s' not supported for 'track' element.\n",
                  OGRGeometryTypeToName(poGeom->getGeometryType()) );
        if (nullptr != poCT)
            delete poGeom;
        return OGRERR_FAILURE;
    }
    }

    if (nullptr != poCT)
        delete poGeom;

    return OGRERR_NONE;
}

//.........这里部分代码省略.........

    /*-----------------------------------------------------------------
     * MITAB won't accept new features unless they are in a type derived
     * from TABFeature... so we have to do our best to map to the right
     * feature type based on the geometry type.
     *----------------------------------------------------------------*/
    poGeom = poFeature->GetGeometryRef();
    if( poGeom != NULL )
        eGType = poGeom->getGeometryType();
    else
        eGType = wkbNone;

    switch( wkbFlatten(eGType) )
    {
    /*-------------------------------------------------------------
     * POINT
     *------------------------------------------------------------*/
    case wkbPoint:
        poTABFeature = new TABPoint(poFeature->GetDefnRef());
        if(poFeature->GetStyleString())
        {
            poTABPointFeature = (TABPoint*)poTABFeature;
            poTABPointFeature->SetSymbolFromStyleString(
                poFeature->GetStyleString());
        }
        break;
    /*-------------------------------------------------------------
     * REGION
     *------------------------------------------------------------*/
    case wkbPolygon:
    case wkbMultiPolygon:
        poTABFeature = new TABRegion(poFeature->GetDefnRef());
        if(poFeature->GetStyleString())
        {
            poTABRegionFeature = (TABRegion*)poTABFeature;
            poTABRegionFeature->SetPenFromStyleString(
                poFeature->GetStyleString());

            poTABRegionFeature->SetBrushFromStyleString(
                poFeature->GetStyleString());
        }
        break;
    /*-------------------------------------------------------------
     * LINE/PLINE/MULTIPLINE
     *------------------------------------------------------------*/
    case wkbLineString:
    case wkbMultiLineString:
        poTABFeature = new TABPolyline(poFeature->GetDefnRef());
        if(poFeature->GetStyleString())
        {
            poTABPolylineFeature = (TABPolyline*)poTABFeature;
            poTABPolylineFeature->SetPenFromStyleString(
                poFeature->GetStyleString());
        }
        break;
    /*-------------------------------------------------------------
     * Collection types that are not directly supported... convert
     * to multiple features in output file through recursive calls.
     *------------------------------------------------------------*/
    case wkbGeometryCollection:
    case wkbMultiPoint:
    {
        OGRErr eStatus = OGRERR_NONE;
        int i;
        OGRGeometryCollection *poColl = (OGRGeometryCollection*)poGeom;
        OGRFeature *poTmpFeature = poFeature->Clone();

        for (i=0; eStatus==OGRERR_NONE && i<poColl->getNumGeometries(); i++)
        {
            poTmpFeature->SetGeometry(poColl->getGeometryRef(i));
            eStatus = CreateFeature(poTmpFeature);
        }
        delete poTmpFeature;
        return eStatus;
    }
    break;
    /*-------------------------------------------------------------
     * Unsupported type.... convert to MapInfo geometry NONE
     *------------------------------------------------------------*/
    case wkbUnknown:
    default:
        poTABFeature = new TABFeature(poFeature->GetDefnRef());
        break;
    }

    if( poGeom != NULL )
        poTABFeature->SetGeometryDirectly(poGeom->clone());

    for (int i=0; i< poFeature->GetDefnRef()->GetFieldCount(); i++)
    {
        poTABFeature->SetField(i,poFeature->GetRawFieldRef( i ));
    }


    eErr = CreateFeature(poTABFeature);

    delete poTABFeature;

    return eErr;
}

OGRErr      OGRGFTTableLayer::ISetFeature( OGRFeature *poFeature )
{
    GetLayerDefn();

    if (!poDS->IsReadWrite())
    {
        CPLError(CE_Failure, CPLE_AppDefined,
                 "Operation not available in read-only mode");
        return OGRERR_FAILURE;
    }

    if (osTableId.size() == 0)
    {
        CPLError(CE_Failure, CPLE_NotSupported,
                "Cannot set feature to non-created table");
        return OGRERR_FAILURE;
    }

    if (poDS->GetAccessToken().size() == 0)
    {
        CPLError(CE_Failure, CPLE_AppDefined,
                 "Operation not available in unauthenticated mode");
        return OGRERR_FAILURE;
    }

    if (poFeature->GetFID() == OGRNullFID)
    {
        CPLError( CE_Failure, CPLE_AppDefined,
                  "FID required on features given to SetFeature()." );
        return OGRERR_FAILURE;
    }

    CPLString      osCommand;

    osCommand += "UPDATE ";
    osCommand += osTableId;
    osCommand += " SET ";

    int iField;
    int nFieldCount = poFeatureDefn->GetFieldCount();
    for(iField = 0; iField < nFieldCount + bHiddenGeometryField; iField++)
    {
        if (iField > 0)
            osCommand += ", ";

        if (iField == nFieldCount)
        {
            osCommand += EscapeAndQuote(GetGeometryColumn());
        }
        else
        {
            const char* pszFieldName =
                poFeatureDefn->GetFieldDefn(iField)->GetNameRef();
            osCommand += EscapeAndQuote(pszFieldName);
        }

        osCommand += " = ";

        OGRGeometry* poGeom = poFeature->GetGeometryRef();
        if (iGeometryField != iLatitudeField && iField == iGeometryField &&
            (iField == nFieldCount || poGeom != NULL || !poFeature->IsFieldSet( iField )))
        {
            if (poGeom == NULL)
                osCommand += "''";
            else
            {
                char* pszKML;
                if (poGeom->getSpatialReference() != NULL &&
                    !poGeom->getSpatialReference()->IsSame(poSRS))
                {
                    OGRGeometry* poGeom4326 = poGeom->clone();
                    poGeom4326->transformTo(poSRS);
                    pszKML = poGeom4326->exportToKML();
                    delete poGeom4326;
                }
                else
                {
                    pszKML = poGeom->exportToKML();
                }
                osCommand += "'";
                osCommand += pszKML;
                osCommand += "'";
                CPLFree(pszKML);
            }
            continue;
        }

        if( !poFeature->IsFieldSet( iField ) )
        {
            osCommand += "''";
        }
        else
        {
            OGRFieldType eType = poFeatureDefn->GetFieldDefn(iField)->GetType();
            if (eType != OFTInteger && eType != OFTReal)
            {
                CPLString osTmp;
                const char* pszVal = poFeature->GetFieldAsString(iField);

                if (!CPLIsUTF8(pszVal, -1))
//.........这里部分代码省略.........

OGRErr OGRGFTTableLayer::ICreateFeature( OGRFeature *poFeature )

{
    if (!poDS->IsReadWrite())
    {
        CPLError(CE_Failure, CPLE_AppDefined,
                 "Operation not available in read-only mode");
        return OGRERR_FAILURE;
    }

    if (osTableId.size() == 0)
    {
        CreateTableIfNecessary();
        if (osTableId.size() == 0)
        {
            CPLError(CE_Failure, CPLE_NotSupported,
                    "Cannot add feature to non-created table");
            return OGRERR_FAILURE;
        }
    }

    if (poDS->GetAccessToken().size() == 0)
    {
        CPLError(CE_Failure, CPLE_AppDefined,
                 "Operation not available in unauthenticated mode");
        return OGRERR_FAILURE;
    }

    CPLString      osCommand;

    osCommand += "INSERT INTO ";
    osCommand += osTableId;
    osCommand += " (";

    int iField;
    int nFieldCount = poFeatureDefn->GetFieldCount();
    for(iField = 0; iField < nFieldCount; iField++)
    {
        if (iField > 0)
            osCommand += ", ";

        const char* pszFieldName =
            poFeatureDefn->GetFieldDefn(iField)->GetNameRef();
        osCommand += EscapeAndQuote(pszFieldName);
    }
    if (bHiddenGeometryField)
    {
        if (iField > 0)
            osCommand += ", ";
        osCommand += EscapeAndQuote(GetGeometryColumn());
    }
    osCommand += ") VALUES (";
    for(iField = 0; iField < nFieldCount + bHiddenGeometryField; iField++)
    {
        if (iField > 0)
            osCommand += ", ";

        OGRGeometry* poGeom = poFeature->GetGeometryRef();
        /* If there's a geometry, let's use it in priority over the textual */
        /* content of the field. */
        if (iGeometryField != iLatitudeField && iField == iGeometryField &&
            (iField == nFieldCount || poGeom != NULL || !poFeature->IsFieldSet( iField )))
        {
            if (poGeom == NULL)
                osCommand += "''";
            else
            {
                char* pszKML;
                if (poGeom->getSpatialReference() != NULL &&
                    !poGeom->getSpatialReference()->IsSame(poSRS))
                {
                    OGRGeometry* poGeom4326 = poGeom->clone();
                    poGeom4326->transformTo(poSRS);
                    pszKML = poGeom4326->exportToKML();
                    delete poGeom4326;
                }
                else
                {
                    pszKML = poGeom->exportToKML();
                }
                osCommand += "'";
                osCommand += pszKML;
                osCommand += "'";
                CPLFree(pszKML);
            }
            continue;
        }

        if( !poFeature->IsFieldSet( iField ) )
        {
            osCommand += "''";
        }
        else
        {
            OGRFieldType eType = poFeatureDefn->GetFieldDefn(iField)->GetType();
            if (eType != OFTInteger && eType != OFTReal)
            {
                CPLString osTmp;
                const char* pszVal = poFeature->GetFieldAsString(iField);

//.........这里部分代码省略.........

OGRErr GTMTrackLayer::CreateFeature (OGRFeature *poFeature)
{
    FILE* fpTmpTrackpoints = poDS->getTmpTrackpointsFP();
    if (fpTmpTrackpoints == NULL)
        return CE_Failure;

    FILE* fpTmpTracks = poDS->getTmpTracksFP();
    if (fpTmpTracks == NULL)
        return CE_Failure;

    OGRGeometry *poGeom = poFeature->GetGeometryRef();
    if ( poGeom == NULL )
    {
        CPLError( CE_Failure, CPLE_AppDefined, 
                  "Features without geometry not supported by GTM writer in track layer." );
        return OGRERR_FAILURE;
    }
   
    if (NULL != poCT)
    {
        poGeom = poGeom->clone();
        poGeom->transform( poCT );
    }
    
    switch( poGeom->getGeometryType() )
    {
    case wkbLineString:
    case wkbLineString25D:
    {
        WriteFeatureAttributes(poFeature);
        OGRLineString* line = (OGRLineString*)poGeom;
        for(int i = 0; i < line->getNumPoints(); ++i)
        {
            double lat = line->getY(i);
            double lon = line->getX(i);
            float altitude = 0;
            CheckAndFixCoordinatesValidity(lat, lon);
            poDS->checkBounds((float)lat, (float)lon);
            if (line->getGeometryType() == wkbLineString25D)
	      altitude = (float)line->getZ(i);
            WriteTrackpoint( lat, lon, altitude, i==0 );
        }
        break;
    }

    case wkbMultiLineString:
    case wkbMultiLineString25D:
    {
        int nGeometries = ((OGRGeometryCollection*)poGeom)->getNumGeometries ();
        for(int j = 0; j < nGeometries; ++j)
        {
            WriteFeatureAttributes(poFeature);
            OGRLineString* line = (OGRLineString*) ( ((OGRGeometryCollection*)poGeom)->getGeometryRef(j) );
            int n = (line) ? line->getNumPoints() : 0;
            for(int i = 0; i < n; ++i)
            {
                double lat = line->getY(i);
                double lon = line->getX(i);
                float altitude = 0;
                CheckAndFixCoordinatesValidity(lat, lon);
                if (line->getGeometryType() == wkbLineString25D)
		  altitude = (float) line->getZ(i);
                WriteTrackpoint( lat, lon, altitude, i==0 );
            }
        }
        break;
    }
    
    default:
    {
        CPLError( CE_Failure, CPLE_NotSupported,
                  "Geometry type of `%s' not supported for 'track' element.\n",
                  OGRGeometryTypeToName(poGeom->getGeometryType()) );
        if (NULL != poCT)
            delete poGeom;
        return OGRERR_FAILURE;
    }
    }
    
    if (NULL != poCT)
        delete poGeom;

    return OGRERR_NONE;
}

int main (int argc, const char * argv[]) {
  
  if (argc < 2 || argc > 3 || strcmp(argv[1], "-h") == 0 || strcmp(argv[1], "--help") == 0) {
    std::cout << "=== prepair Help ===\n" << std::endl;
    std::cout << "Usage:   triface 'POLYGON(...)'" << std::endl;
    std::cout << "OR" << std::endl;
    std::cout << "Usage:   triface -f infile.txt (infile.txt must contain one WKT on the 1st line)" << std::endl;
    return 0;
  }
  
  // Read input
  unsigned int bufferSize = 10000000;
  char *inputWKT = (char *)malloc(bufferSize*sizeof(char *));
  
  for (int argNum = 1; argNum < argc; ++argNum) {
    if (strcmp(argv[argNum], "-f") == 0) {
      
      if (argNum + 1 <= argc - 1 && argv[argNum+1][0] != '-') {
        std::ifstream infile(argv[argNum+1], std::ifstream::in);
        infile.getline(inputWKT, bufferSize);
        ++argNum;
      } else {
        std::cerr << "Error: Missing input file name." << std::endl;
        return 1;
      }
    }
    else 
      strcpy(inputWKT, argv[argNum]);
  }
  
//  std::cout << "Processing: " << inputWKT << std::endl;
  
  OGRGeometry *geometry;
  OGRGeometryFactory::createFromWkt(&inputWKT, NULL, &geometry);
  if (geometry == NULL) {
    std::cout << "Error: WKT is not valid" << std::endl;
    return 1;
  }

  if (geometry->getGeometryType() != wkbPolygon25D) {
    std::cout << "Error: input geometry is not a 3D polygon" << std::endl;
    return 1;
  }
  
  //-- project to proper plane + get flattened geometry
  int proj = get_projection_plane(geometry);
  OGRGeometry *flatgeom = geometry->clone();
  if (proj == 1) {
    OGRPolygon *polygon = (OGRPolygon *)flatgeom;
    for (int curp = 0; curp < polygon->getExteriorRing()->getNumPoints(); ++curp) 
      polygon->getExteriorRing()->setPoint(curp, polygon->getExteriorRing()->getX(curp), polygon->getExteriorRing()->getZ(curp), 0);
    for (int currentRing = 0; currentRing < polygon->getNumInteriorRings(); ++currentRing) {
      for (int curp = 0; curp < polygon->getInteriorRing(currentRing)->getNumPoints(); ++curp)
        polygon->getInteriorRing(currentRing)->setPoint(curp, polygon->getInteriorRing(currentRing)->getX(curp), polygon->getInteriorRing(currentRing)->getZ(curp), 0);
    }
  }
  else if (proj == 0) {
    OGRPolygon *polygon = (OGRPolygon *)geometry;
    for (int curp = 0; curp < polygon->getExteriorRing()->getNumPoints(); ++curp)
      polygon->getExteriorRing()->setPoint(curp, polygon->getExteriorRing()->getY(curp), polygon->getExteriorRing()->getZ(curp), 0);
    for (int currentRing = 0; currentRing < polygon->getNumInteriorRings(); ++currentRing) {
      for (int curp = 0; curp < polygon->getInteriorRing(currentRing)->getNumPoints(); ++curp)
        polygon->getInteriorRing(currentRing)->setPoint(curp, polygon->getInteriorRing(currentRing)->getY(curp), polygon->getInteriorRing(currentRing)->getZ(curp), 0);
    }
  }
  flatgeom->flattenTo2D();
//  std::cout << "geom: " << geometry->getCoordinateDimension() << std::endl;
//  std::cout << "flatgeom: " << flatgeom->getCoordinateDimension() << std::endl;
  
  //-- check if flattened geometry is valid
  if (flatgeom->IsValid() == FALSE) {
    std::cout << "Error: input polygon is not valid." << std::endl;
    return 1;
  }
  
  if (proj == 2) {
    Triangulationxy triangulation;
    triangulateandtag_xy(geometry, triangulation);
    for (Triangulationxy::Finite_faces_iterator currentFace = triangulation.finite_faces_begin(); currentFace != triangulation.finite_faces_end(); ++currentFace) {
      std::cout << "--triangle--" << std::endl;
      Point p = currentFace->vertex(0)->point();
      std::cout << p.x() << ", " << p.y() << ", " << p.z() << std::endl;
      p = currentFace->vertex(1)->point();
      std::cout << p.x() << ", " << p.y() << ", " << p.z() << std::endl;
      p = currentFace->vertex(2)->point();
      std::cout << p.x() << ", " << p.y() << ", " << p.z() << std::endl;
    }
  }
  else if (proj == 1) {
    Triangulationxz triangulation;
    triangulateandtag_xz(geometry, triangulation);
    for (Triangulationxz::Finite_faces_iterator currentFace = triangulation.finite_faces_begin(); currentFace != triangulation.finite_faces_end(); ++currentFace) {
      std::cout << "--triangle--" << std::endl;
      Point p = currentFace->vertex(0)->point();
      std::cout << p.x() << ", " << p.y() << ", " << p.z() << std::endl;
      p = currentFace->vertex(1)->point();
      std::cout << p.x() << ", " << p.y() << ", " << p.z() << std::endl;
      p = currentFace->vertex(2)->point();
      std::cout << p.x() << ", " << p.y() << ", " << p.z() << std::endl;
    }
//.........这里部分代码省略.........

int Raster::VectortoRaster(const char * sVectorSourcePath,
                   const char * sRasterOutputPath,
                   const char * psFieldName,
                   RasterMeta * p_rastermeta ){

    OGRRegisterAll();

    OGRDataSource * pDSVectorInput;
    pDSVectorInput = OGRSFDriverRegistrar::Open( sVectorSourcePath, FALSE );
    if (pDSVectorInput == NULL)
        return INPUT_FILE_ERROR;

    OGRLayer * poLayer = pDSVectorInput->GetLayer(0);

    // The type of the field.
    OGRFeature * feat1 = poLayer->GetFeature(0);
    int fieldindex = feat1->GetFieldIndex(psFieldName);
    OGRFieldType fieldType = feat1->GetFieldDefnRef(fieldindex)->GetType();
    OGRFeature::DestroyFeature( feat1 );

    // The data type we're going to use for the file
    GDALDataType OutputDataType = GDT_Byte;

    // Handle field types according to their type:
    switch (fieldType) {
    case OFTString:
        CSVWriteVectorValues(poLayer, psFieldName, sRasterOutputPath);
        break;
    case OFTInteger:
        break;
    case OFTReal:
        OutputDataType = GDT_Float64;
    default:
        throw RasterManagerException(VECTOR_FIELD_NOT_VALID, "Type of field not recognized.");
        break;
    }

    // Get our projection and set the rastermeta accordingly.
    // -------------------------------------------------------
    char *pszWKT = NULL;

    OGRSpatialReference* poSRS = poLayer->GetSpatialRef();

    poSRS->exportToWkt(&pszWKT);
    p_rastermeta->SetProjectionRef(pszWKT);
    CPLFree(pszWKT);
    OSRDestroySpatialReference(poSRS);

    // Create the output dataset for writing
    GDALDataset * pDSOutput = CreateOutputDS(sRasterOutputPath, p_rastermeta);

    // Create a list of burn-in values
    // -------------------------------------------------------


    std::vector<OGRGeometryH> ogrBurnGeometries;
    std::vector<double> dBurnValues;

    poLayer->ResetReading();
    OGRFeature * ogrFeat = poLayer->GetNextFeature();   // <------- DRMEM!!! UNADDRESSABLE ACCESS of freed memory

    while( ogrFeat != NULL ){

        OGRGeometry * ogrGeom = ogrFeat->GetGeometryRef();

        // No geometry found. Move along.
        if( ogrGeom == NULL )
        {
            OGRFeature::DestroyFeature( ogrFeat );
            continue;
        }

        OGRGeometry * geoClone = ogrGeom->clone();  // <------- DRMEM!!! UNADDRESSABLE ACCESS of freed memory AND LEAK 20 direct bytes

        // Push a clone of this geometry onto the list of shapes to burn
        ogrBurnGeometries.push_back( (OGRGeometryH) geoClone );

        if (fieldType == OFTString){
            // If it's a string type we burn the FID. The value is then placed in a CSV file
            dBurnValues.push_back( ogrFeat->GetFID() );
        }
        else {
            // If it's a float type or a byte type we write it directly.
            dBurnValues.push_back( ogrFeat->GetFieldAsDouble(psFieldName) );
        }
        // GetNextFeature() creates a clone so we must delete it.
        OGRFeature::DestroyFeature( ogrFeat );      // <------- DRMEM!!! UNADDRESSABLE ACCESS beyond heap bounds:
        ogrFeat = poLayer->GetNextFeature();        // <------- DRMEM!!! UNADDRESSABLE ACCESS of freed memory
    }

    // Do the Actual Burning of Geometries.
    // -------------------------------------------------------

    int band = 1;
    char **papszRasterizeOptions = NULL;
    papszRasterizeOptions =
        CSLSetNameValue( papszRasterizeOptions, "ALL_TOUCHED", "TRUE" );

    CPLErr err = GDALRasterizeGeometries( pDSOutput, 1, &band,
                                          ogrBurnGeometries.size(),
//.........这里部分代码省略.........

//.........这里部分代码省略.........

                    /* add last point */
                    ((OGRLineString *) poGeom)->getPoint(0, &pt);
                    poGeomString.addPoint(&pt);

                }
            }
            else if (strlen(ftype) > 2 && EQUALN(ftype, "15", 2)) { /* -> circle with radius */
                float r;
                char s[3]; /* 15 */

                r = 0;
                if (2 != sscanf(ftype, "%s %f", s, &r) || r < 0) {
                    CPLDebug("OGR-VFK", "%s: invalid circle (unknown or negative radius) "
                             "fid = " CPL_FRMT_GIB, m_poDataBlock->GetName(), m_nFID);
                    m_bValid = FALSE;
                }
                else {
                    double c_x, c_y;
                    
                    ((OGRLineString *) poGeom)->getPoint(0, &pt);
                    c_x = pt.getX();
                    c_y = pt.getY();

                    /* define first point on a circle */
                    pt.setX(c_x + r);
                    pt.setY(c_y);
                    poGeomString.addPoint(&pt);

                    /* define second point on a circle */
                    pt.setX(c_x);
                    pt.setY(c_y + r);
                    poGeomString.addPoint(&pt);

                    /* define third point on a circle */
                    pt.setX(c_x - r);
                    pt.setY(c_y);
                    poGeomString.addPoint(&pt);

                    /* define fourth point on a circle */
                    pt.setX(c_x);
                    pt.setY(c_y - r);
                    poGeomString.addPoint(&pt);

                    /* define last point (=first) on a circle */
                    pt.setX(c_x + r);
                    pt.setY(c_y);
                    poGeomString.addPoint(&pt);
                }

            }
            else if (EQUAL(ftype, "11")) {                          /* curve */
                int npoints;

                npoints = ((OGRLineString *) poGeom)->getNumPoints();
                if (npoints > 2) { /* circular otherwise line string */
                    for (int i = 0; i < npoints; i++) {
                        ((OGRLineString *) poGeom)->getPoint(i, &pt);
                        poGeomString.addPoint(&pt);
                    }
                }
            }

            if (!poGeomString.IsEmpty())
                poGeomCurved = poGeomString.CurveToLine();
            
            if (poGeomCurved) {
                int npoints;
                
                npoints = ((OGRLineString *) poGeomCurved)->getNumPoints();
                CPLDebug("OGR-VFK", "%s: curve (type=%s) to linestring (npoints=%d) fid = " CPL_FRMT_GIB,
                         m_poDataBlock->GetName(), ftype,
                         npoints, m_nFID);
                if (npoints > 1)
                    m_paGeom = (OGRGeometry *) poGeomCurved->clone();
                delete poGeomCurved;
            }
        }

        if (!m_paGeom) {
            /* check degenerated linestrings */
            if (m_nGeometryType == wkbLineString) {
                int npoints;
                
                npoints = ((OGRLineString *) poGeom)->getNumPoints();
                if (npoints < 2) {
                    CPLError(CE_Warning, CPLE_AppDefined,
                             "%s: invalid linestring (%d vertices) fid = " CPL_FRMT_GIB,
                             m_poDataBlock->GetName(), npoints, m_nFID);
                    m_bValid = FALSE;
                }
            }

            if (m_bValid)
                m_paGeom = (OGRGeometry *) poGeom->clone(); /* make copy */
        }
    }

    return m_bValid;
}

OGRDataSource *AoIIntersection::intersectAoIWithLayers ( OGRDataSource *ogrSourceData, OGRPolygon *aoiPoly, IntersectionSummary *summary, const char *outFmt )
{
  OGRDataSource *ogrIntersection = NULL;

  // Spatial reference setup
  // make a spatial reference for the area of interest polygon
  OGRSpatialReference aoiRef;
  aoiRef.SetWellKnownGeogCS( "WGS84" );

  // make a spatial reference for the coord sys we will use to calculate area in acres - Albers USA equal area conic
  bool acreageCalcAvailable = true;
  char *aecWkt = "PROJCS[\"USA_Contiguous_Lambert_Conformal_Conic\",GEOGCS[\"GCS_North_American_1983\",DATUM[\"North_American_Datum_1983\",SPHEROID[\"GRS_1980\",6378137,298.257222101]],PRIMEM[\"Greenwich\",0],UNIT[\"Degree\",0.017453292519943295]],PROJECTION[\"Lambert_Conformal_Conic_2SP\"],PARAMETER[\"False_Easting\",0],PARAMETER[\"False_Northing\",0],PARAMETER[\"Central_Meridian\",-96],PARAMETER[\"Standard_Parallel_1\",33],PARAMETER[\"Standard_Parallel_2\",45],PARAMETER[\"Latitude_Of_Origin\",39],UNIT[\"Meter\",1],AUTHORITY[\"EPSG\",\"102004\"]]";
  OGRSpatialReference aecRef;
  OGRErr ogrErr = aecRef.importFromWkt( &aecWkt );
  if ( ogrErr != OGRERR_NONE )
  {
    setError ( NO_SPATIAL_REFERENCE );
    acreageCalcAvailable = false;
  }

  // begin creating the output data structure
  // OGRDataSource is the root
  ogrIntersection = buildIntersectionDataSource( outFmt );
  if (! ogrIntersection )
  {
    setError( NO_OUTPUT_DATASOURCE );
    return 0;
  }

  int acreIndex = 0, areaIndex = 0;
  OGRFieldDefn *acreFldDefn = NULL, *areaPctFldDefn = NULL;
  OGRFeatureDefn *featureDefn = buildFeatureDefinition( acreIndex, areaIndex, acreFldDefn, areaPctFldDefn );

  // walk the layers in the input data
  //
  OGRLayer *inputLayer;
  summary->numLayers = ogrSourceData->GetLayerCount();

  for (int layerCt = 0; layerCt < summary->numLayers; ++layerCt)
  {
    inputLayer = ogrSourceData->GetLayer( layerCt );
    if ( inputLayer == NULL )
    {
      setError( NO_INPUT_LAYER );
      // clean up
      delete ogrIntersection;
      return 0;
    }
      
    // make a clone of aoi polygon to be manipulated
    OGRPolygon *aoiClone = (OGRPolygon *)aoiPoly->clone();
    if ( ! aoiClone )
    {
      setError( NO_AOI_CLONE );
      // clean up
      delete ogrIntersection;
      return 0;
    }

    // ensure that the area of interest polygon is in the same spatial reference as the data layer
    // find the spatial reference for the layer
    OGRSpatialReference *dataRef = inputLayer->GetSpatialRef();
    if ( dataRef )
    {
      OGRCoordinateTransformation *aoiTransform = OGRCreateCoordinateTransformation( &aoiRef, dataRef );
      if( aoiTransform == NULL )
      {
        setError( NO_AOI_TRANSFORM );
        // clean up
        delete ogrIntersection;
        delete aoiClone;
        return 0;
      }
      aoiClone->transform( aoiTransform );
      delete aoiTransform;
    }

    // find the transform from data layer's CS to Albers USA  
    // for acreage calculation
    OGRCoordinateTransformation *aecTransform = NULL;
    acreageCalcAvailable = false;
    if ( dataRef )
    {
      aecTransform = OGRCreateCoordinateTransformation( dataRef, &aecRef );
      if( aecTransform == NULL )
      {
        setError( NO_ACRE_TRANSFORM );
      }
      else
        acreageCalcAvailable = true;
    }

    // the area enclosed by the AoI
    // used for computing the percentage of the AoI intersected by polygons
    summary->aoiArea = aoiClone->getExteriorRing()->get_Area();

    // create a layer for outputting the intersecting polygons
    OGRLayer *intersectionLayer = ogrIntersection->CreateLayer( inputLayer->GetLayerDefn()->GetName(), dataRef, wkbPolygon, 0 );
    if ( ! intersectionLayer )
    {
//.........这里部分代码省略.........

//.........这里部分代码省略.........
                {
                    OGRMakeWktCoordinate(szLowerCorner, sGeomBounds.MinX, sGeomBounds.MinY, sGeomBounds.MinZ, nCoordDimension);
                    OGRMakeWktCoordinate(szUpperCorner, sGeomBounds.MaxX, sGeomBounds.MaxY, sGeomBounds.MaxZ, nCoordDimension);
                }
                if (bWriteSpaceIndentation)
                    VSIFPrintfL(fp, "      ");
                poDS->PrintLine( fp, "<gml:boundedBy><gml:Envelope%s%s><gml:lowerCorner>%s</gml:lowerCorner><gml:upperCorner>%s</gml:upperCorner></gml:Envelope></gml:boundedBy>",
                                (nCoordDimension == 3) ? " srsDimension=\"3\"" : "",pszSRSName, szLowerCorner, szUpperCorner);
                CPLFree(pszSRSName);
            }

            char** papszOptions = (bIsGML3Output) ? CSLAddString(NULL, "FORMAT=GML3") : NULL;
            if (bIsGML3Output && !poDS->IsLongSRSRequired())
                papszOptions = CSLAddString(papszOptions, "GML3_LONGSRS=NO");
            const char* pszSRSDimensionLoc = poDS->GetSRSDimensionLoc();
            if( pszSRSDimensionLoc != NULL )
                papszOptions = CSLSetNameValue(papszOptions, "SRSDIMENSION_LOC", pszSRSDimensionLoc);
            if (poDS->IsGML32Output())
            {
                if( poFeatureDefn->GetGeomFieldCount() > 1 )
                    papszOptions = CSLAddString(papszOptions,
                        CPLSPrintf("GMLID=%s.%s." CPL_FRMT_GIB,
                                   poFeatureDefn->GetName(),
                                   poFieldDefn->GetNameRef(),
                                   poFeature->GetFID()));
                else
                    papszOptions = CSLAddString(papszOptions,
                        CPLSPrintf("GMLID=%s.geom." CPL_FRMT_GIB,
                                   poFeatureDefn->GetName(), poFeature->GetFID()));
            }
            if( !bIsGML3Output && OGR_GT_IsNonLinear(poGeom->getGeometryType()) )
            {
                OGRGeometry* poGeomTmp = OGRGeometryFactory::forceTo(
                    poGeom->clone(),OGR_GT_GetLinear(poGeom->getGeometryType()));
                pszGeometry = poGeomTmp->exportToGML(papszOptions);
                delete poGeomTmp;
            }
            else
                pszGeometry = poGeom->exportToGML(papszOptions);
            CSLDestroy(papszOptions);
            if (bWriteSpaceIndentation)
                VSIFPrintfL(fp, "      ");
            if( bRemoveAppPrefix )
                poDS->PrintLine( fp, "<%s>%s</%s>",
                                 poFieldDefn->GetNameRef(),
                                 pszGeometry,
                                 poFieldDefn->GetNameRef() );
            else
                poDS->PrintLine( fp, "<%s:%s>%s</%s:%s>",
                                 pszPrefix, poFieldDefn->GetNameRef(),
                                 pszGeometry,
                                 pszPrefix, poFieldDefn->GetNameRef() );
            CPLFree( pszGeometry );
        }
    }

    // Write all "set" fields.
    for( int iField = 0; iField < poFeatureDefn->GetFieldCount(); iField++ )
    {

        OGRFieldDefn *poFieldDefn = poFeatureDefn->GetFieldDefn( iField );

        if( poFeature->IsFieldSet( iField ) && iField != nGMLIdIndex )
        {
            OGRFieldType eType = poFieldDefn->GetType();
            if (eType == OFTStringList )

//.........这里部分代码省略.........
            {
		        if(pTrackCancel)
                    pTrackCancel->PutMessage(wxString::Format(_("Band %d requested, but only bands 1 to %d available."), panBandList[i], pSrcRasterDataSet->GetBandCount()), -1, enumGISMessageErr);
                return false;
            }
        }
    }

    bool bDefaultfilter(false);
    if(pFilter == NULL)
    {
        pFilter = new wxGxRasterFilter(enumRasterTiff);
        bDefaultfilter = true;
    }

    CPLString szDriver(pFilter->GetDriver().mb_str());
    CPLString szExt(pFilter->GetExt().mb_str());
    CPLString szBaseName = CPLGetBasename(pSrcRasterDataSet->GetPath());

	GDALDriver* pDriver = (GDALDriver*)GDALGetDriverByName( szDriver );
	if( pDriver == NULL )
    {
		if(pTrackCancel)
            pTrackCancel->PutMessage(wxString::Format(_("Output driver '%s' not recognised."), szDriver.c_str()), -1, enumGISMessageErr);
        if(bDefaultfilter)
            wxDELETE(pFilter);
		return false;
    }

    pSrcFeatureDataSet->Reset();
    OGRFeatureSPtr pFeature;
    size_t nNameCounter(nCounterBegin);
    while((pFeature = pSrcFeatureDataSet->Next()) != NULL)
    {
        if(pTrackCancel && !pTrackCancel->Continue())
        {
            wxString sErr(_("Interrupted by user"));
            CPLString sFullErr(sErr.mb_str());
            CPLError( CE_Warning, CPLE_AppDefined, sFullErr );

            if(pTrackCancel)
                pTrackCancel->PutMessage(wxString(sFullErr, wxConvLocal), -1, enumGISMessageErr);

            if(bDefaultfilter)
                wxDELETE(pFilter);
            return false;
        }

        OGRGeometry *pGeom = pFeature->GetGeometryRef();
		if(wkbFlatten(pSrcFeatureDataSet->GetGeometryType()) != wkbUnknown && !pGeom)
            continue;

		OGRGeometry *pNewGeom(NULL);
        if( !bSame && poCT )
        {
            if(pGeom)
            {
				pNewGeom = pGeom->clone();
                OGRErr eErr = pNewGeom->transform(poCT);
                if(eErr != OGRERR_NONE)
	                wxDELETE(pNewGeom);
            }
        }
		else
			pNewGeom = pGeom->clone();

		OGREnvelope GeomEnv;
		pNewGeom->getEnvelope(&GeomEnv);
		OGREnvelope RasterEnv = pSrcRasterDataSet->GetEnvelope();
		GeomEnv.Intersect(RasterEnv);
		if(GeomEnv.IsInit())
		{
            CPLString szPath;
			if(bUseCounter)
            {
                szPath.Printf("%s_%d", szBaseName.c_str(), nNameCounter++);
                CPLString sNewName(CheckUniqName(szDstFolderPath, szPath, szExt).mb_str(wxConvUTF8));
                szPath = CPLFormFilename(szDstFolderPath, sNewName, szExt);
            }
            else
            {
                CPLString szName = pFeature->GetFieldAsString(nFieldNo);
                CPLString sNewName(CheckUniqName(szDstFolderPath, szName, szExt).mb_str(wxConvUTF8));
                szPath = CPLFormFilename(szDstFolderPath, sNewName, szExt);
            }
			CreateSubRaster(pSrcRasterDataSet, GeomEnv, pNewGeom, pDriver, szPath, eOutputType, nBandCount, panBandList, dfOutResX, dfOutResY, bCopyNodata, bSkipSourceMetadata, papszOptions, pTrackCancel);
		}
        nCounter++;
        if(pProgressor)
            pProgressor->SetValue(nCounter);
    }

    if(poCT)
        OCTDestroyCoordinateTransformation(poCT);

    if(bDefaultfilter)
        wxDELETE(pFilter);

	return true;
}