# Using GDAL to get vertices of a multipolygon

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?

• You can get the number of polygons with `getNumGeometries()` and get each polygon with `getGeometryRef(i)`. Commented Feb 23, 2015 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. Commented Feb 23, 2015 at 23:12

## 3 Answers

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
``````

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
Commented Feb 24, 2015 at 16:14

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.

}

}

}

}

}

}
``````