5

I would like to know if there is a Python library that can be used to convert GeoJSON Multi-polygon to a Polygon.

Note: i want to convert .geojson file Multi-Polygon feature type to Polygon feature type.

  • 1
    What have you already tried / looked at? How did those attempts go? If you tell us more background, we can avoid telling you things you already know. – BradHards Oct 18 '15 at 19:38
  • shapely could help – Zoltan Oct 18 '15 at 20:19
  • @BradHards, i have fair knowledge of GIS, am a software developer, i am working on a software that does extraction of .geojson file content, parse, and load into Postgresql to be used by a web app. i want to run Douglas Pecker Algorithm en.wikipedia.org/wiki/… on the polygon but the Python module, rdp that am using does not support Multi-polygon – Jideobi Benedine Ofomah Oct 19 '15 at 15:56
5
from osgeo import ogr

str = '''
{ "type": "MultiPolygon",
    "coordinates": [
        [
            [[40, 40], [20, 45], [45, 30], [40, 40]]
        ],
        [
            [[20, 35], [10, 30], [10, 10], [30, 5], [45, 20], [20, 35]],
            [[30, 20], [20, 15], [20, 25], [30, 20]]
        ]
    ]
}
'''
g = ogr.CreateGeometryFromJson(str) 

print "Hi! I'm a %s with an Area  %s" % (g.GetGeometryName(), g.Area())
print "I have inside me %s feature(s)!\n" % g.GetGeometryCount()
for idx, f in enumerate(g):
    print "I'm feature n.%s and I am a %s.\t I have an Area of %s - You can get my json repr with f.ExportToJson()" % (idx, f.GetGeometryName(),f.Area())

>>> Hi! I'm a MULTIPOLYGON with an Area 712.5
>>> I have inside me 2 feature(s)!

>>> I'm feature n.0 and I am a POLYGON. I have an Area of 87.5 - You can get my json repr with f.ExportToJson()
>>> I'm feature n.1 and I am a POLYGON. I have an Area of 625.0 - You can get my json repr with f.ExportToJson()

3

This is a quick way to get a list of single part geometries from a multipart feature using shapely. The buffer is in there to correct any invalid geometry in both the original feature and in the exploded parts.

import json
from shapely.geometry import shape, mapping
# Sample data
str = '''
{ "type": "MultiPolygon",
 "coordinates": [
     [[[40, 40], [20, 45], [45, 30], [40, 40]]],
     [[[20, 35], [10, 30], [10, 10], [30, 5], [45, 20], [20, 35]],[[30, 20], [20, 15], [20, 25], [30, 20]]]
 ]
}
'''
g = [x.buffer(0) for x in shape(json.loads(str)).buffer(0).geoms]

What is going on here:

First, starting from the inside, json.loads() loads our geojson string into a python object. shape() then turns this into a shapely feature, a Multipolygon in this case. buffer(0) on this object corrects any invalid geometry (for example, turning a single-part bowtie into a multipolygon). geoms is a property of the resulting geometry from the buffer operation that contains a collection of all the geometries in the multipart geometry. Iterate over this ([x.buffer(0) for x in ...]), again calling buffer on each geometry, to construct a list of valid single part geometries.

What is missing here is attaching any attributes to each of the new geometry parts. If needed, I can add that in a little bit (once I get a better sample put together).

2

In PyQGIS, you only need to convert from Multi Part to Single Part. I tried out my approach with the Multi Part Polygon of next image:

enter image description here

The following code did the work:

layer = iface.mapCanvas().currentLayer()
features = layer.getFeatures()

for feature in features:
    geom = feature.geometry()
    print geom.exportToGeoJSON()

    # check if feature geometry is multipart
    if geom.isMultipart():
        new_features = []
        temp_feature = QgsFeature(feature)

        # set the geometry of each part
        for part in geom.asGeometryCollection ():
            temp_feature.setGeometry(part)
            new_features.append(QgsFeature(temp_feature)) 

for feature in new_features:
    print feature.geometry().exportToWkt()
    print feature.geometry().exportToGeoJSON()

After running the code at the Python Console I got:

enter image description here

The first print in the red rectangles is for the GeoJSON Multi-polygon and the other are for single polygons in GeoJSON format.

To test the effective conversion to single parts, I used the WKT notation for the geometries with the QuickWKT plugin (see next image):

enter image description here

The Multi Part geometry was effectively converted to Single Part. I can use these features to get a memory layer or a shapefile.

1

There are many Python modules to do that: all those which use the geointerface protocol (GeoJSON-like protocol for geo-spatial (GIS) vector data, proposed by Sean Gillies)

Shapely is one of them but we begin with PyGeoIf ("shapely ultralight", easiest to install)

from pygeoif import geometry
import geojson
str = '''
{ "type": "MultiPolygon",
 "coordinates": [
     [[[40, 40], [20, 45], [45, 30], [40, 40]]],
     [[[20, 35], [10, 30], [10, 10], [30, 5], [45, 20], [20, 35]],[[30, 20], [20, 15], [20, 25], [30, 20]]]
 ]
}
'''
mp = geojson.loads(str)
Multipoly = geometry.as_shape(mp)
for poly in Multipoly.geoms:
     print poly
POLYGON((40.0 40.0, 20.0 45.0, 45.0 30.0, 40.0 40.0))
POLYGON((20.0 35.0, 10.0 30.0, 10.0 10.0, 30.0 5.0, 45.0 20.0, 20.0 35.0),(30.0 20.0, 20.0 15.0, 20.0 25.0, 30.0 20.0))
# __geointerface__
for poly in Multipoly.geoms:
    print poly.__geo_interface__
{'type': 'Polygon', 'coordinates': (((40.0, 40.0), (20.0, 45.0), (45.0, 30.0), (40.0, 40.0)),)}
{'type': 'Polygon', 'coordinates': (((20.0, 35.0), (10.0, 30.0), (10.0, 10.0), (30.0, 5.0), (45.0, 20.0), (20.0, 35.0)), ((30.0, 20.0), (20.0, 15.0), (20.0, 25.0), (30.0, 20.0)))}

Likewise, you could also use Shapely

 from shapely.geometry import shape, mapping
 Multipoly = shape(mp)
 for poly in Multipoly:
     print poly.wkt
 POLYGON ((40 40, 20 45, 45 30, 40 40))
 POLYGON ((20 35, 10 30, 10 10, 30 5, 45 20, 20 35), (30 20, 20 15, 20 25, 30 20))
for poly in Multipoly:
    print mapping(poly)
{'type': 'Polygon', 'coordinates': (((40.0, 40.0), (20.0, 45.0), (45.0, 30.0), (40.0, 40.0)),)}
{'type': 'Polygon', 'coordinates': (((20.0, 35.0), (10.0, 30.0), (10.0, 10.0), (30.0, 5.0), (45.0, 20.0), (20.0, 35.0)), ((30.0, 20.0), (20.0, 15.0), (20.0, 25.0), (30.0, 20.0)))}

And you can use the QGIS QuickWKT plugin to show the results (as Xulnik)

enter image description here

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.