1

I am successfully using GDAL to get the vertices of a polygon from a ShapeFile using .NET code. This link showed me how to do that:

http://geoinformaticstutorial.blogspot.com/2012/10/accessing-vertices-from-polygon-with.html

However if the shape type is a multipolygon, then this process does not get any vertices. How can I get the vertices of a multipolygon?

  • 1
    You can get the number of polygons with getNumGeometries() and get each polygon with getGeometryRef(i). – Mike T Feb 23 '15 at 22:16
  • and then follow the example for each part in the polygon. Polylines can also do the same if they have multiple parts. – Michael Stimson Feb 23 '15 at 23:12
0

I suppose that you are happy with your script but you may find use for the Spatialite SQL dialect in the future.

Let's take a shapefile with one multipolygon to start with.

ogrinfo multipolygon.shp -al
INFO: Open of `multipolygon.shp'
      using driver `ESRI Shapefile' successful.

Layer name: multipolygon
Geometry: Polygon
Feature Count: 1
Extent: (403784.354400, 6924130.250945) - (409732.186830, 6926861.991572)
Layer SRS WKT:
(unknown)
OGRFeature(multipolygon):0
  MULTIPOLYGON (((403784.35439963144 6926261.3961154055,405838.0033819516 692686
1.9915724993,407000.44620213285 6925680.1747053144,405140.53768984287 6924692.09
83081609,403784.35439963144 6926261.3961154055)),((407484.79737720836 6926745.74
72904809,409732.18682955875 6926745.7472904809,409344.70588949829 6924323.991415
103,406922.95001412078 6924130.2509450726,407484.79737720836 6926745.7472904809)
))

Next utilize the ST_DissolvePoints function of Spatialite https://www.gaia-gis.it/gaia-sins/xmlBlob/spatialite-sql-latest.html through SQLite SQL dialect of GDAL http://www.gdal.org/ogr_sql_sqlite.html. You must have GDAL that is compiled with new enough Spatialite because the ST_DissolvePoints function is pretty new.

ogrinfo multipolygon.shp -dialect sqlite -sql "select dissolvepoints(geometry) from multipolygon"
INFO: Open of `multipolygon.shp'
      using driver `ESRI Shapefile' successful.

Layer name: SELECT
Geometry: Unknown (any)
Feature Count: 1
Extent: (403784.354400, 6924130.250945) - (409732.186830, 6926861.991572)
Layer SRS WKT:
(unknown)
Geometry Column = dissolvepoints(geometry)
OGRFeature(SELECT):0
  MULTIPOINT (403784.35439963144 6926261.3961154055,405838.0033819516 6926861.99
15724993,407000.44620213285 6925680.1747053144,405140.53768984287 6924692.098308
1609,403784.35439963144 6926261.3961154055,407484.79737720836 6926745.7472904809
,409732.18682955875 6926745.7472904809,409344.70588949829 6924323.991415103,4069
22.95001412078 6924130.2509450726,407484.79737720836 6926745.7472904809)

I know that the result is not exactly what you aim at but perhaps you will find some use for the SQLite dialect in the future.

  • Thanks, but I'm going with what I know especially since it's working for me. Hopefully someone finds your post useful. – Rono Feb 24 '15 at 16:14
0

Here is the code I used to get this to work. This is VB.NET code:

Public Function SchoolDistrictBoundary() As List(Of List(Of Coordinate))
Dim district As Feature = schoolDistrict()
Dim result = New List(Of List(Of Coordinate))

If district IsNot Nothing Then
  Dim geometry = district.GetGeometryRef()
  Dim type = geometry.GetGeometryType()

  If type = wkbGeometryType.wkbPolygon Then
    'only one polygon
    Dim coordSet = New List(Of Coordinate)
    Dim ring = geometry.GetGeometryRef(0)
    For x = 0 To ring.GetPointCount() - 1
      Dim point() As Double = {0, 0}
      ring.GetPoint(x, point)

      Dim coord = New Coordinate()
      coord.Longitude = point(0)
      coord.Latitude = point(1)
      coordSet.Add(coord)
    Next
    result.Add(coordSet)

  ElseIf type = wkbGeometryType.wkbMultiPolygon Then
    'multiple polygons
    Dim boundaries = geometry.GetBoundary()
    For shape = 0 To boundaries.GetGeometryCount() - 1
      Dim coordSet = New List(Of Coordinate)
      Dim polygon = boundaries.GetGeometryRef(shape)
      For x = 0 To polygon.GetPointCount() - 1
        Dim point() As Double = {0, 0}
        polygon.GetPoint(x, point)

        Dim coord = New Coordinate()
        coord.Longitude = point(0)
        coord.Latitude = point(1)

        coordSet.Add(coord)
      Next

      result.Add(coordSet)
    Next

  End If
End If
Return result
End Function
0

Here's a C# version, ds is a OSGeo.OGR.Driver driver and layerIdX is the layer you want:

OSGeo.OGR.Feature feature = null;
OSGeo.OGR.Geometry geo = null, ring = null;

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

int layerIdX = 2;

OSGeo.OGR.Layer layer = ds.GetLayerByIndex(layerIdX);

if( geoType == "wkbMultiPolygon") 
{

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

        feature = layer.GetFeature(i);

        if( feature != null )
        {

            geo = feature.GetGeometryRef();

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

                var inGeo = geo.GetGeometryRef(j);

                for( int k = 0; k < inGeo.GetGeometryCount(); ++k )
                {

                    ring = inGeo.GetGeometryRef(k);

                    int pointCount = ring.GetPointCount();

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

                        ring.GetPoint( l, pointList );
                        // pointList[0] is the Longitude.
                        // pointList[1] is the Latitude.
                        // pointList[2] is the Altitude. 

                    }

                }

            }

        }

    }

}

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.