Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

How to extract vertexes of geometries in ESRI shapefiles using ogr library with c++

I have answered the question and also changed the title of the question from How to access vertices in a polygon layer using ogr library with c++ to How to extract vertexes of geometries in ESRI shapefiles using ogr library with c++. The code can be used to show shapes with OpenGL.
Any suggestion to make it better will be appreciated.
For example do you think that I should use projection of the layer in order to enhance the view in OpenGL? How can I use the projection and apply it to the OpenGL window?

up to now I have tried these classes in order to extract the vertices of polygons:

  • OGRPolygon : It doesn't have any method to get number of vertices, start point and end point.
  • OGRLineString : when I use getNumPoints() method, I get 1 for all of the polygon features in my layer.
  • OGRLinearRing : when I use getNumPoints() method, I get 1 for all of the polygon features in my layer.

what should I do to get the number of vertices and coordinates of vertices in a polygon shapefile in order to draw the polygon in an OpenGL window.

like image 408
Sepideh Abadpour Avatar asked Sep 11 '13 17:09

Sepideh Abadpour


People also ask

How do I get polygon coordinates from a shapefile?

If your end goal is to extract it to a spreadsheet, you can right-click your layer, select Save As... , change the Format to csv and change the GEOMETRY to AS_WKT . This will create a csv file with a field containing the coordinates of all vertices for each polygon feature :) Most welcome, glad it helped!


2 Answers

I found the answer:
You should care some information about ESRI shapefiles.

  • Point layer : Can have only features of type wkbpoint. So you can define a data structure named MyPoint2D and store the coordinates of this layer in a vector with members of type MyPoint2D.

  • Multipoint layer : Can have only features of type wkbMultipoint. So you can store each feature in a vector with members of type MyPoint2D and a collection of such vectors will hold the coordinates of all of the features in the layer.

  • Polyline layer : Can have features of type wkbLineString and wkbMultiLineString. So you can have a data structure MyLine2D, which is a vector with MyPoint2D members and a vector with members of type MyLine2D will hold the LineFeature and finally the features in the layer will be stored in a vector with members of type LineFeature.

  • Polygon Layer : About the polygons you should be aware that each polygon is consisted of an Exterior ring and several Interior rings ( this was my real problem which caused me to ask the question ). And also the layer can have features of type wkbPolygon and wkbMultiPolygon.

So at first we should define these data structures:

//data structure for points
typedef struct MyPoint2D
{
double dX;
double dY;
}MyPoint2D;

//Holds Coordinates of Point Shapefile
vector<MyPoint2D>PointLayer;

//data structure for a multipoint feature
typedef struct MultipointFeature
{
vector<MyPoint2D>PointsOfFeature;
}MultipointFeature;

//Holds Coordinates of multiPoint Shapefile
vector<MultipointFeature> MultipointLayer;

//data structure for lines
typedef struct MyLine2D
{
 vector<MyPoint2D> LineString;
}MyLine2D;

//data structure for a line feature
typedef struct LineFeature
{
vector<MyLine2D>LinesOfFeature;
}LineFeature;

//Holds Coordinates of Line Shapefile
vector<LineFeature> LineLayer;

//data structure for rings
typedef struct MyRing2D
{
vector<MyPoint2D> RingString;
bool IsClockwised;
}MyRing2D;

//data structure for polygons
typedef struct MyPolygon2D
{
vector<MyRing2D>Polygon;
}MyPolygon2D;

//data structure for a polygon feature
typedef struct PolygonFeature
{
vector<MyPolygon2D>PolygonsOfFeature;
}PolygonFeature;

//Holds Coordinates of Polygon Shapefile
vector<PolygonFeature> PolygonLayer;

//data structure to hold bounding box
typedef struct SBoundingBox
{
 float fMaxX;
 float fMaxY;
 float fMinX;
 float fMinY;
 }SBoundingBox;

//Bounding Box of Shapefile 
SBoundingBox sBoundingBox;  

And this is the whole code, I have used in order to extract the coordinates of vertexes in ESRI's shapefiles:

void OpenShapeFile(char* filename)
{
OGRErr error;
OGRDataSource *poDataSource;
poDataSource = OGRSFDriverRegistrar::Open(filename,false);
OGRLayer *poLayer;
poLayer = poDataSource ->GetLayer(0);
OGREnvelope Envelope;
error = poLayer ->GetExtent(&Envelope,true);
sBoundingBox.fMaxX = Envelope.MaxX;
sBoundingBox.fMaxY = Envelope.MaxY;
sBoundingBox.fMinX = Envelope.MinX;
sBoundingBox.fMinY = Envelope.MinY;  

OGRwkbGeometryType LayerGeometryType = poLayer ->GetGeomType();
int NumberOfFeatures = poLayer ->GetFeatureCount(true);
poLayer ->ResetReading();  

//Point Shapefile
if ( wkbFlatten ( LayerGeometryType ) == wkbPoint )
{
   OGRFeature *poFeature;
   for ( int i = 0; i < NumberOfFeatures; i++ )
   {
       poFeature = poLayer ->GetNextFeature();
       OGRGeometry *poGeometry;
       poGeometry = poFeature ->GetGeometryRef();
       if ( poGeometry != NULL && wkbFlatten ( poGeometry ->getGeometryType() ) == wkbPoint )
       {
           OGRPoint *poPoint = ( OGRPoint * )poGeometry;
           MyPoint2D pt;
           pt.dX = poPoint ->getX();
           pt.dY = poPoint ->getY();
           PointLayer.push_back(pt);
       }
       OGRFeature::DestroyFeature(poFeature);
   }
}  

//Multipoint Shapefile
if ( wkbFlatten ( LayerGeometryType ) == wkbMultiPoint )
{
   OGRFeature *poFeature;
   MultipointFeature MultiPoint;
   for ( int i = 0; i < NumberOfFeatures; i++ )
   {
       poFeature = poLayer ->GetNextFeature();
       OGRGeometry *poGeometry;
       poGeometry = poFeature ->GetGeometryRef();
       if ( poGeometry != NULL && wkbFlatten ( poGeometry ->getGeometryType() ) == wkbMultiPoint )
       {
           OGRMultiPoint *poMultipoint = ( OGRMultiPoint * )poGeometry;
           int NumberOfGeometries = poMultipoint ->getNumGeometries();
           MultiPoint.PointsOfFeature.resize( NumberOfGeometries );
           for ( int j = 0; j < NumberOfGeometries; j++ )
           {
               OGRGeometry *poPointGeometry = poMultipoint ->getGeometryRef(j);
               OGRPoint *poPoint = ( OGRPoint * )poPointGeometry;
               MyPoint2D pt;
               pt.dX = poPoint ->getX();
               pt.dY = poPoint ->getY();
               MultiPoint.PointsOfFeature.at(j) = pt;
           }
           MultipointLayer.push_back(MultiPoint);
       }
       OGRFeature::DestroyFeature(poFeature);
   }
}  

//Polyline Shapefile
if ( wkbFlatten ( LayerGeometryType ) == wkbLineString )
{
   OGRFeature *poFeature;
   LineFeature Polyline;
   OGRPoint ptTemp;
   for ( int i = 0; i < NumberOfFeatures; i++ )
   {
       poFeature = poLayer ->GetNextFeature();
       OGRGeometry *poGeometry;
       poGeometry = poFeature ->GetGeometryRef();
       if ( poGeometry != NULL && wkbFlatten ( poGeometry ->getGeometryType() ) == wkbLineString  )
       {
           OGRLineString *poLineString = ( OGRLineString * )poGeometry;
           Polyline.LinesOfFeature.resize(1);
           int NumberOfVertices = poLineString ->getNumPoints();
           Polyline.LinesOfFeature.at(0).LineString.resize(NumberOfVertices);
           for ( int k = 0; k < NumberOfVertices; k++ )
           {
               poLineString ->getPoint(k,&ptTemp);
               MyPoint2D pt;
               pt.dX = ptTemp.getX();
               pt.dY = ptTemp.getY();
               Polyline.LinesOfFeature.at(0).LineString.at(k) = pt;
           }
           LineLayer.push_back(Polyline);
       }
       else if ( poGeometry != NULL && wkbFlatten ( poGeometry ->getGeometryType() ) == wkbMultiLineString )
       {
           OGRMultiLineString *poMultiLineString = ( OGRMultiLineString * )poGeometry;
           int NumberOfGeometries = poMultiLineString ->getNumGeometries();
           Polyline.LinesOfFeature.resize(NumberOfGeometries);
           for ( int j = 0; j < NumberOfGeometries; j++ )
           {
               OGRGeometry *poLineGeometry = poMultiLineString ->getGeometryRef(j);
               OGRLineString *poLineString = ( OGRLineString * )poLineGeometry;
               int NumberOfVertices = poLineString ->getNumPoints();
               Polyline.LinesOfFeature.at(j).LineString.resize(NumberOfVertices);
               for ( int k = 0; k < NumberOfVertices; k++ )
               {
                   poLineString ->getPoint(k,&ptTemp);
                   MyPoint2D pt;
                   pt.dX = ptTemp.getX();
                   pt.dY = ptTemp.getY();
                   Polyline.LinesOfFeature.at(j).LineString.at(k) = pt;
               }
           }
           LineLayer.push_back(Polyline);
       }
       OGRFeature::DestroyFeature(poFeature);
   }
}  

//Polygon Shapefile
if ( wkbFlatten ( LayerGeometryType ) == wkbPolygon )
{
   OGRFeature *poFeature;
   PolygonFeature Polygon;
   OGRPoint ptTemp;
   for ( int i = 0; i < NumberOfFeatures; i++ )
   {
       poFeature = poLayer ->GetNextFeature();
       OGRGeometry *poGeometry;
       poGeometry = poFeature ->GetGeometryRef();
       if ( poGeometry != NULL && wkbFlatten ( poGeometry ->getGeometryType() ) == wkbPolygon )
       {
           OGRPolygon *poPolygon = ( OGRPolygon * )poGeometry;
           Polygon.PolygonsOfFeature.resize(1);
           int NumberOfInnerRings = poPolygon ->getNumInteriorRings();
           OGRLinearRing *poExteriorRing = poPolygon ->getExteriorRing();
           Polygon.PolygonsOfFeature.at(0).Polygon.resize(NumberOfInnerRings+1);
           Polygon.PolygonsOfFeature.at(0).Polygon.at(0).IsClockwised = poExteriorRing ->isClockwise();
           int NumberOfExteriorRingVertices = poExteriorRing ->getNumPoints();
           Polygon.PolygonsOfFeature.at(0).Polygon.at(0).RingString.resize(NumberOfExteriorRingVertices);
           for ( int k = 0; k < NumberOfExteriorRingVertices; k++ )
           {
               poExteriorRing ->getPoint(k,&ptTemp);
               MyPoint2D pt;
               pt.dX = ptTemp.getX();
               pt.dY = ptTemp.getY();
               Polygon.PolygonsOfFeature.at(0).Polygon.at(0).RingString.at(k) = pt;
           }
           for ( int h = 1; h <= NumberOfInnerRings; h++ )
           {
               OGRLinearRing *poInteriorRing = poPolygon ->getInteriorRing(h-1);
               Polygon.PolygonsOfFeature.at(0).Polygon.at(h).IsClockwised = poInteriorRing ->isClockwise();
               int NumberOfInteriorRingVertices = poInteriorRing ->getNumPoints();
               Polygon.PolygonsOfFeature.at(0).Polygon.at(h).RingString.resize(NumberOfInteriorRingVertices);
               for ( int k = 0; k < NumberOfInteriorRingVertices; k++ )
               {
                   poInteriorRing ->getPoint(k,&ptTemp);
                   MyPoint2D pt;
                   pt.dX = ptTemp.getX();
                   pt.dY = ptTemp.getY();
                   Polygon.PolygonsOfFeature.at(0).Polygon.at(h).RingString.at(k) = pt;
               }
           }
               PolygonLayer.push_back(Polygon);
       }
       else if ( poGeometry != NULL && wkbFlatten ( poGeometry ->getGeometryType() ) == wkbMultiPolygon )
       {
           OGRMultiPolygon *poMultiPolygon = ( OGRMultiPolygon * )poGeometry;
           int NumberOfGeometries = poMultiPolygon ->getNumGeometries();
           Polygon.PolygonsOfFeature.resize(NumberOfGeometries);
           for ( int j = 0; j < NumberOfGeometries; j++ )
           {
               OGRGeometry *poPolygonGeometry = poMultiPolygon ->getGeometryRef(j);
               OGRPolygon *poPolygon = ( OGRPolygon * )poPolygonGeometry;
               int NumberOfInnerRings = poPolygon ->getNumInteriorRings();
               OGRLinearRing *poExteriorRing = poPolygon ->getExteriorRing();
               Polygon.PolygonsOfFeature.at(j).Polygon.resize(NumberOfInnerRings+1);
               Polygon.PolygonsOfFeature.at(j).Polygon.at(0).IsClockwised = poExteriorRing ->isClockwise();
               int NumberOfExteriorRingVertices = poExteriorRing ->getNumPoints();
               Polygon.PolygonsOfFeature.at(j).Polygon.at(0).RingString.resize(NumberOfExteriorRingVertices);
               for ( int k = 0; k < NumberOfExteriorRingVertices; k++ )
               {
                   poExteriorRing ->getPoint(k,&ptTemp);
                   MyPoint2D pt;
                   pt.dX = ptTemp.getX();
                   pt.dY = ptTemp.getY();
                   Polygon.PolygonsOfFeature.at(j).Polygon.at(0).RingString.at(k) = pt;
               }
               for ( int h = 1; h <= NumberOfInnerRings; h++ )
               {
                   OGRLinearRing *poInteriorRing = poPolygon ->getInteriorRing(h-1);
                   Polygon.PolygonsOfFeature.at(j).Polygon.at(h).IsClockwised = poInteriorRing ->isClockwise();
                   int NumberOfInteriorRingVertices = poInteriorRing ->getNumPoints();
                   Polygon.PolygonsOfFeature.at(j).Polygon.at(h).RingString.resize(NumberOfInteriorRingVertices);
                   for ( int k = 0; k < NumberOfInteriorRingVertices; k++ )
                   {
                       poInteriorRing ->getPoint(k,&ptTemp);
                       MyPoint2D pt;
                       pt.dX = ptTemp.getX();
                       pt.dY = ptTemp.getY();
                       Polygon.PolygonsOfFeature.at(j).Polygon.at(h).RingString.at(k) = pt;
                   }
               }
           }
           PolygonLayer.push_back(Polygon);
       }
   }
   OGRFeature::DestroyFeature(poFeature);
}  

OGRDataSource::DestroyDataSource(poDataSource);
}
like image 77
Sepideh Abadpour Avatar answered Oct 28 '22 19:10

Sepideh Abadpour


scott_f wanted a C# implementation and probably someone in the future may also need it so I'll post it as an answer. I am lazy and don't want to rewrite the whole thing using nested lists and different data structures to the ones provided by Rhino/Grasshopper as this implementation is part of a plugin (GH_Structure, IGH_GeometricGoo, IGH_Goo are part of the Rhino/Grasshopper API). But basically IGH_Gooin this case refers to bools, IGH_GeometricGooto a geometry class that can contain both Polyline and Point3dclasses and GH_Structureis a multi-dimensional data-structure (List of Lists) that "adds" items through the .Append() method, where the first parameter is the data to add and the second is a "path", so .Append( item, new GH_Path( 0, 3 ) ) will add the item to the first list of lists at the 4th index. It's probably a good idea to do wkbFlatten, so that it works when the dimension of the wkbType is different, but I'll leave that to the reader as I didn't need it for my implementation. Note that the layer is passed as a parameter, so if you need all the data in the file, you should loop through all the layers.

class WKBToRhinoGeoUtilities
{

    private string input;
    private int layerNumber;
    private int EPSGCode;
    private int InEPSGCode;
    private OSGeo.OGR.Driver driver;
    OSGeo.OGR.DataSource dataSource;

    private GH_Path path;
    private double[] pointList = { 0, 0, 0 };

    private long numberOfFeatures;
    private int layerCount;
    private OSGeo.OGR.Layer layer;
    private OSGeo.OGR.wkbGeometryType wkbGeoType;
    private int currentProjection;
    private OSGeo.OSR.CoordinateTransformation transform;
    private OSGeo.OGR.Feature feature;
    private OSGeo.OGR.Geometry geo, ring;

    public string debug = "Working!";

    public WKBToRhinoGeoUtilities( string input, int layerNumber, int EPSGCode, 
                                   int InEPSGCode, OSGeo.OGR.DataSource dataSource )
    {

        this.input = input;
        this.layerNumber = layerNumber;
        this.EPSGCode = EPSGCode;
        this.InEPSGCode = InEPSGCode;
        this.dataSource = dataSource;

    }

    public string Setup()
    {

        GetDataSourceAttributes();

        IdentifyEPSG();

        return debug;

    }

    public GH_Structure<IGH_GeometricGoo> Run( ref GH_Structure<IGH_Goo> cullPattern )
    {

        switch( wkbGeoType )
        {

            case OSGeo.OGR.wkbGeometryType.wkbUnknown:

                return WKBUnknown( ref cullPattern );

            case OSGeo.OGR.wkbGeometryType.wkbPolygon:

                return WKBPolygon( ref cullPattern );

            case OSGeo.OGR.wkbGeometryType.wkbLineString:

                return WKBLineString( ref cullPattern );

            case OSGeo.OGR.wkbGeometryType.wkbPoint:

                return WKBPoint( ref cullPattern );

            case OSGeo.OGR.wkbGeometryType.wkbPoint25D:

                return WKBPoint( ref cullPattern );

            case OSGeo.OGR.wkbGeometryType.wkbMultiPolygon:

                return WKBMultiPolygon( ref cullPattern );

            case OSGeo.OGR.wkbGeometryType.wkbMultiLineString:

                return WKBMultiLineString( ref cullPattern );

            case OSGeo.OGR.wkbGeometryType.wkbGeometryCollection:

                return WKBGeometryCollection( ref cullPattern );

            default:

                return null;

        }

    }

    private void GetDataSourceAttributes()
    {

        if( dataSource != null )
        { 

            layerCount = dataSource.GetLayerCount();
            layer = dataSource.GetLayerByIndex( layerNumber % layerCount );
            wkbGeoType = layer.GetGeomType();
            numberOfFeatures = unchecked( ( int ) layer.GetFeatureCount( 0 ) );

        }

        else
        {

            debug = "Data Source is null!";

        }

    }

    private void IdentifyEPSG()
    {

        // Source projections, try to get it automatically otherwise resort to user input.
        OSGeo.OSR.SpatialReference source = null;
        int identifyEPSG = 0;

        try
        {

            source = layer.GetSpatialRef();
            identifyEPSG = source.AutoIdentifyEPSG();

        }

        catch
        {

            debug = "The file's EPSG code couldn't be found automatically, it will default to EPSG:3857, make sure you input the correct one!";
            source = new OSGeo.OSR.SpatialReference( "" );
            identifyEPSG = source.ImportFromEPSG( InEPSGCode );

        }

        currentProjection = Convert.ToInt32( source.GetAttrValue( "AUTHORITY", 1 ) );

        // Reprojections.
        OSGeo.OSR.SpatialReference target = new OSGeo.OSR.SpatialReference( "" );

        if( currentProjection != EPSGCode )
        { 

            target.ImportFromEPSG( EPSGCode );
            transform = new OSGeo.OSR.CoordinateTransformation( source, target );

        }

    }

    private GH_Structure<IGH_GeometricGoo> WKBUnknown( ref GH_Structure<IGH_Goo> cullPattern )
    {

        GH_Structure<IGH_GeometricGoo>  geoOutput = new GH_Structure<IGH_GeometricGoo>();
        cullPattern = new GH_Structure<IGH_Goo>();

        double[] pointList = { 0, 0, 0 };

        int counter = 0;

        do
        {

            path = new GH_Path( counter );

            feature = layer.GetNextFeature();

            bool pattern = false;

            if( feature != null )
            { 

                OSGeo.OGR.Geometry geo = feature.GetGeometryRef();

                if( geo != null )
                {

                    if( currentProjection != EPSGCode )
                    {

                        geo.Transform( transform );

                    }

                    pattern = true;

                    int pointCount = geo.GetPointCount();

                    IGH_GeometricGoo geoGoo = null;

                    Point3d[] tempPointArray = new Point3d[pointCount];

                    for( int i = 0; i < pointCount; ++i )
                    {

                        geo.GetPoint( i, pointList );

                        tempPointArray[i] = new Point3d( pointList[0], pointList[1], pointList[2] );

                    }

                    if( tempPointArray.Length > 2 )
                    {

                        Polyline tempPoly = new Polyline( tempPointArray );
                        geoGoo = GH_Convert.ToGeometricGoo( tempPoly );
                        geoOutput.Append( geoGoo, path );

                    }

                    else if( tempPointArray.Length > 0 )
                    {

                        geoGoo = GH_Convert.ToGeometricGoo( tempPointArray[0] );
                        geoOutput.Append( geoGoo, path );

                    }

                }

            }

            IGH_Goo patternGoo = GH_Convert.ToGoo( pattern );
            cullPattern.Append( patternGoo, path );

            counter++;

        } while( feature != null );

        return geoOutput;

    }

    private GH_Structure<IGH_GeometricGoo> WKBPolygon( ref GH_Structure<IGH_Goo> cullPattern )
    {

        GH_Structure<IGH_GeometricGoo> geoOutput = new GH_Structure<IGH_GeometricGoo>();
        cullPattern = new GH_Structure<IGH_Goo>();

        if ( numberOfFeatures != -1 )
        { 

            for( int i = 0; i < numberOfFeatures; ++i )
            {

                path = new GH_Path(i);

                feature = layer.GetFeature(i);

                bool pattern = false;

                if( feature != null )
                {

                    geo = feature.GetGeometryRef();

                    if( geo != null )
                    {

                        IGH_GeometricGoo geoGoo = WKBPolygonSingle( geo, ref pattern );
                        geoOutput.Append( geoGoo, path );

                    }

                }

                IGH_Goo gooPattern = GH_Convert.ToGoo( pattern );
                cullPattern.Append( gooPattern, path );

            }

        }

        else
        {

            feature = layer.GetNextFeature();

            int counter = 0;

            while( feature != null )
            {

                path = new GH_Path( counter );

                geo = feature.GetGeometryRef();

                bool pattern = false;

                if( geo != null )
                {

                    IGH_GeometricGoo geoGoo = WKBPolygonSingle( geo, ref pattern );
                    geoOutput.Append( geoGoo, path );

                }

                IGH_Goo gooPattern = GH_Convert.ToGoo( pattern );

                cullPattern.Append( gooPattern, path );

                feature = layer.GetNextFeature();

                counter++;

            }

        }

        return geoOutput;

    }

    private GH_Structure<IGH_GeometricGoo> WKBLineString( ref GH_Structure<IGH_Goo> cullPattern )
    {

        GH_Structure<IGH_GeometricGoo> geoOutput = new GH_Structure<IGH_GeometricGoo>();
        cullPattern = new GH_Structure<IGH_Goo>();

        if ( numberOfFeatures != -1 )
        { 

            for( int i = 0; i < numberOfFeatures; ++i )
            {

                path = new GH_Path(i);
                feature = layer.GetFeature(i);

                bool pattern = false;

                if( feature != null )
                {

                    geo = feature.GetGeometryRef();

                    if( geo != null )
                    { 

                        IGH_GeometricGoo geoGoo = WKBLineStringSingle( geo, ref pattern );
                        geoOutput.Append( geoGoo, path );

                    } 

                }

                IGH_Goo gooPattern = GH_Convert.ToGoo( pattern );
                cullPattern.Append( gooPattern, path );

            }

        }

        else
        {

            int counter = 0;

            feature = layer.GetNextFeature();

            while( feature != null )
            {

                path = new GH_Path( counter );

                geo = feature.GetGeometryRef();

                bool pattern = false;

                if( geo != null )
                { 

                    IGH_GeometricGoo geoGoo = WKBLineStringSingle( geo, ref pattern );
                    geoOutput.Append( geoGoo, path );

                }

                IGH_Goo gooPattern = GH_Convert.ToGoo( pattern );
                cullPattern.Append( gooPattern, path );

                feature = layer.GetNextFeature();
                counter++;

            }

        }

        return geoOutput;

    }

    private GH_Structure<IGH_GeometricGoo> WKBPoint( ref GH_Structure<IGH_Goo> cullPattern )
    {

        GH_Structure<IGH_GeometricGoo> geoOutput = new GH_Structure<IGH_GeometricGoo>();
        cullPattern = new GH_Structure<IGH_Goo>();

        if ( numberOfFeatures != -1 )
        { 

            for( int i = 0; i < numberOfFeatures; ++i )
            {

                path = new GH_Path(i);

                feature = layer.GetFeature(i);

                bool tempPattern = feature != null ? true : false;

                IGH_Goo gooPattern = GH_Convert.ToGoo( tempPattern );

                cullPattern.Append( gooPattern, path );

                if( feature != null )
                {

                    geo = feature.GetGeometryRef();

                    IGH_GeometricGoo geoGoo = WKBPointSingle( geo );

                    geoOutput.Append( geoGoo, path );

                }

            }

        }

        else
        {

            int counter = 0;

            feature = layer.GetNextFeature();

            while( feature != null )
            {

                path = new GH_Path( counter );

                bool tempPattern = feature != null ? true : false;

                IGH_Goo gooPattern = GH_Convert.ToGoo( tempPattern );

                cullPattern.Append( gooPattern, path );

                geo = feature.GetGeometryRef();

                IGH_GeometricGoo geoGoo = WKBPointSingle( geo );

                geoOutput.Append( geoGoo, path );

                feature = layer.GetNextFeature();

                counter++;

            }

        }

        return geoOutput;

    }

    private GH_Structure<IGH_GeometricGoo> WKBMultiPolygon( ref GH_Structure<IGH_Goo> cullPattern )
    {

        GH_Structure<IGH_GeometricGoo> geoOutput = new GH_Structure<IGH_GeometricGoo>();
        cullPattern = new GH_Structure<IGH_Goo>();

        if ( numberOfFeatures != -1 )
        { 

            for( int i = 0; i < numberOfFeatures; ++i )
            {

                path = new GH_Path(i);

                feature = layer.GetFeature(i);

                bool pattern = false;

                if( feature != null )
                {

                    geo = feature.GetGeometryRef();

                    for( int j = 0; j < geo.GetGeometryCount(); ++j )
                    {

                        OSGeo.OGR.Geometry inGeo = geo.GetGeometryRef(j);

                        IGH_GeometricGoo geoGoo = WKBPolygonSingle( inGeo, ref pattern );
                        geoOutput.Append( geoGoo, path );

                    }

                }

                IGH_Goo gooPattern = GH_Convert.ToGoo( pattern );
                cullPattern.Append( gooPattern, path );

            }

        }

        else
        {

            int counter = 0;
            feature = layer.GetNextFeature();

            while( feature != null )
            {

                bool pattern = false;

                path = new GH_Path( counter );

                geo = feature.GetGeometryRef();

                for( int j = 0; j < geo.GetGeometryCount(); ++j )
                {

                    OSGeo.OGR.Geometry inGeo = geo.GetGeometryRef(j);                                        

                    if( inGeo != null )
                    { 

                        IGH_GeometricGoo geoGoo = WKBPolygonSingle( inGeo, ref pattern );
                        geoOutput.Append( geoGoo, path );

                    }

                }

                IGH_Goo gooPattern = GH_Convert.ToGoo( pattern );
                cullPattern.Append( gooPattern, path );

                counter++;
                feature = layer.GetNextFeature();

            }

        }

        return geoOutput;

    }

    private GH_Structure<IGH_GeometricGoo> WKBMultiLineString( ref GH_Structure<IGH_Goo> cullPattern )
    {

        GH_Structure<IGH_GeometricGoo> geoOutput = new GH_Structure<IGH_GeometricGoo>();
        cullPattern = new GH_Structure<IGH_Goo>();

        if ( numberOfFeatures != -1 )
        { 

            for( int i = 0; i < numberOfFeatures; ++i )
            {

                path = new GH_Path(i);
                feature = layer.GetFeature(i);

                bool pattern = false;

                if( feature != null )
                {

                    geo = feature.GetGeometryRef();

                    if( geo != null )
                    { 

                        for( int j = 0; j < geo.GetGeometryCount(); ++j )
                        {

                            ring = geo.GetGeometryRef(j);

                            IGH_GeometricGoo geoGoo = WKBLineStringSingle( ring, ref pattern );
                            geoOutput.Append( geoGoo, path );

                        }

                    }

                }

                IGH_Goo gooPattern = GH_Convert.ToGoo( pattern );
                cullPattern.Append( gooPattern, path );

            }

        }

        else
        {

            feature = layer.GetNextFeature();
            int counter = 0;

            while( feature != null )
            {

                path = new GH_Path( counter );

                geo = feature.GetGeometryRef();

                bool pattern = false;

                if( geo != null )
                { 

                    for( int j = 0; j < geo.GetGeometryCount(); ++j )
                    {

                        ring = geo.GetGeometryRef(j);

                        IGH_GeometricGoo geoGoo = WKBLineStringSingle( ring, ref pattern );
                        geoOutput.Append( geoGoo, path );

                    }

                }

                IGH_Goo gooPattern = GH_Convert.ToGoo( pattern );
                cullPattern.Append( gooPattern, path );

                feature = layer.GetNextFeature();

                counter++;

            }

        }

        return geoOutput;

    }

    private GH_Structure<IGH_GeometricGoo> WKBGeometryCollection( ref GH_Structure<IGH_Goo> cullPattern )
    {

        GH_Structure<IGH_GeometricGoo> geoOutput = new GH_Structure<IGH_GeometricGoo>();
        cullPattern = new GH_Structure<IGH_Goo>();
        GH_Structure<IGH_Goo> cullPatternTemp = new GH_Structure<IGH_Goo>();

        if ( numberOfFeatures != -1 )
        {

            for( int i = 0; i < numberOfFeatures; ++i )
            {

                path = new GH_Path(i);
                feature = layer.GetFeature(i);

                if( feature != null )
                {

                    geo = feature.GetGeometryRef();

                    if( geo != null )
                    {

                        bool tempPattern = feature != null ? true : false;

                        IGH_GeometricGoo geoGoo = null;

                        if( geo != null )
                        {

                            for( int j = 0; j < geo.GetGeometryCount(); ++j )
                            {

                                path = new GH_Path( j, i );

                                ring = geo.GetGeometryRef(i);

                                OSGeo.OGR.wkbGeometryType geoType = ring.GetGeometryType();

                                switch( geoType )
                                {

                                    case OSGeo.OGR.wkbGeometryType.wkbPoint:

                                        geoGoo = WKBPointSingle( ring );

                                        break;

                                    case OSGeo.OGR.wkbGeometryType.wkbPoint25D:

                                        geoGoo = WKBPointSingle( ring );

                                        break;

                                    case OSGeo.OGR.wkbGeometryType.wkbPolygon:

                                        geoGoo = WKBPolygonSingle( ring, ref tempPattern );

                                        break;

                                    case OSGeo.OGR.wkbGeometryType.wkbLineString:

                                        geoGoo = WKBLineStringSingle( ring, ref tempPattern );

                                        break;

                                }

                                geoOutput.Append( geoGoo, path );

                                IGH_Goo gooPattern = GH_Convert.ToGoo( tempPattern );
                                cullPattern.Append( gooPattern, path );

                            }

                        }

                    }

                }

            }

        }

        else
        {

            int counter = 0;

            feature = layer.GetNextFeature();

            while( feature != null )
            {

                geo = feature.GetGeometryRef();

                bool tempPattern = feature != null ? true : false;

                IGH_GeometricGoo geoGoo = null;

                if( geo != null )
                {

                    for( int i = 0; i < geo.GetGeometryCount(); ++i )
                    {

                        path = new GH_Path( counter, i );

                        ring = geo.GetGeometryRef(i);

                        OSGeo.OGR.wkbGeometryType geoType = ring.GetGeometryType();

                        switch( geoType )
                        {

                            case OSGeo.OGR.wkbGeometryType.wkbPoint:

                                geoGoo = WKBPointSingle( ring );

                                break;

                            case OSGeo.OGR.wkbGeometryType.wkbPoint25D:

                                geoGoo = WKBPointSingle( ring );

                                break;

                            case OSGeo.OGR.wkbGeometryType.wkbPolygon:

                                geoGoo = WKBPolygonSingle( ring, ref tempPattern );

                                break;

                            case OSGeo.OGR.wkbGeometryType.wkbLineString:

                                geoGoo = WKBLineStringSingle( ring, ref tempPattern );

                                break;

                        }

                        geoOutput.Append( geoGoo, path );

                        IGH_Goo gooPattern = GH_Convert.ToGoo( tempPattern );
                        cullPattern.Append( gooPattern, path );

                    }

                }

                feature = layer.GetNextFeature();

                counter++;

            }

        }

        return geoOutput;

    }

    private IGH_GeometricGoo WKBPointSingle( OSGeo.OGR.Geometry geo )
    {

        IGH_GeometricGoo geoGoo = null;

        if( currentProjection != EPSGCode )
        {

            geo.Transform( transform );

        }

        geo.GetPoint( 0, pointList );

        Point3d tempPoint = new Point3d( pointList[0], pointList[1], pointList[2] );
        geoGoo = GH_Convert.ToGeometricGoo( tempPoint );

        return geoGoo;

    }

    private IGH_GeometricGoo WKBPolygonSingle( OSGeo.OGR.Geometry geo, ref bool pattern )
    {

        IGH_GeometricGoo geoGoo = null;

        ring = geo.GetGeometryRef( 0 );

        if( currentProjection != EPSGCode )
        {

            ring.Transform( transform );

        }

        int pointCount = ring.GetPointCount();

        Point3d[] tempPointArray = new Point3d[pointCount];

        for (int j = 0; j < pointCount; ++j)
        {

            ring.GetPoint( j, pointList );

            tempPointArray[j] = new Point3d( pointList[0], pointList[1], pointList[2] );

        }

        if( tempPointArray.Length > 1 )
        { 

            Polyline polyOut = new Polyline( tempPointArray );

            geoGoo = GH_Convert.ToGeometricGoo( polyOut );

            pattern = true;

        }

        return geoGoo;

    }

    private IGH_GeometricGoo WKBLineStringSingle( OSGeo.OGR.Geometry geo, ref bool pattern )
    {

        IGH_GeometricGoo geoGoo = null;

        if( currentProjection != EPSGCode )
        {

            geo.Transform( transform );

        }

        int pointCount = geo.GetPointCount();

        Point3d[] tempPointArray = new Point3d[pointCount];

        for( int j = 0; j < pointCount; ++j )
        {

            geo.GetPoint( j, pointList );

            tempPointArray[j] = new Point3d( pointList[0], pointList[1], pointList[2] );

        }

        if( tempPointArray.Length > 1 )
        { 

            Polyline polyOut = new Polyline( tempPointArray );
            geoGoo = GH_Convert.ToGeometricGoo( polyOut );

            pattern = true;

        }

        IGH_Goo gooPattern = GH_Convert.ToGoo( pattern );

        return geoGoo;

    }

}

You'd call it by something like (make sure you call

GdalConfiguration.ConfigureOgr();
GdalConfiguration.ConfigureGdal();

):

OSGeo.OGR.Driver driver = Utilities.GetExtractorDriver( input, ref debug );

OSGeo.OGR.DataSource dataSource = driver.Open( input, 0 );

WKBToRhinoGeoUtilities wkbToRhino = new WKBToRhinoGeoUtilities( input, layerNumber, EPSGCode, InEPSGCode, dataSource );
debug = wkbToRhino.Setup();
AddRuntimeMessage( GH_RuntimeMessageLevel.Blank, debug );

if( debug != "Driver is null!" || debug != "Data Source is null!" || debug != "File does not exist!" )
{

    GH_Structure<IGH_Goo> cullPattern = new GH_Structure<IGH_Goo>();
    GH_Structure<IGH_GeometricGoo> geoOutput = wkbToRhino.Run( ref cullPattern );

    DA.SetDataTree( 0, geoOutput );
    DA.SetDataTree( 1, cullPattern );

}

Where I have a Utilities class with this static method:

public static OSGeo.OGR.Driver GetExtractorDriver( string inputString, ref string debug )
{

    if( File.Exists( inputString ) )
    {

        debug = "File exists!";

        string isRunning = RunGDAL.output;

        OSGeo.OGR.Driver driver = null;

        string[] fileTypeArray = inputString.Split( '.' );
        string fileType = fileTypeArray[fileTypeArray.Length - 1];

        if( fileType == "shp" )
        {

            driver = OSGeo.OGR.Ogr.GetDriverByName( "ESRI Shapefile" );

        }

        else if( fileType == "json" || fileType == "geojson" )
        {

            driver = OSGeo.OGR.Ogr.GetDriverByName( "GeoJSON" );

        }

        else if( fileType == "topojson" )
        {

            driver = OSGeo.OGR.Ogr.GetDriverByName( "TopoJSON" );

        }

        else if( fileType == "gpkg" )
        {

            driver = OSGeo.OGR.Ogr.GetDriverByName( "GPKG" );

        }

        else if( fileType == "dxf" )
        {

            driver = OSGeo.OGR.Ogr.GetDriverByName( "DXF" );

        }

        else if( fileType == "kml" )
        {

            driver = OSGeo.OGR.Ogr.GetDriverByName( "KML" );

        }

        else if( fileType == "osm" )
        {

            driver = OSGeo.OGR.Ogr.GetDriverByName( "OSM" );

        }

        return driver;

    }

    else
    {

        debug = "File does not exist!";

        return null;

    }

}
like image 34
Felipe Gutierrez Avatar answered Oct 28 '22 19:10

Felipe Gutierrez