Using Geometry collections with Fiona to write output of Multipoints and Points?

Question

I am getting the following error on this function:

ValueError: Record's geometry type does not match collection schema's geometry type: 'Point' != 'MultiPoint'

The problem is that I have multipoints mixed with regular points. When I try to write the outputs at I get a schema confusion problem

The record causing problems at 'output.write' is this one (note: that for a MultiPoint it only has one x,y):

{'geometry': {'type': 'MultiPoint', 'coordinates': [(13531245.475704141, 2886003.2689278126)]}, 'type': 'Feature', 'id': '0', 'properties': OrderedDict([(u'H20Status', u'Fair')])}

When I go:

                        if geom['type'] == 'MultiPoint':
                        # break Multipoints into points
                        # i.e. <class 'shapely.geometry.multipoint.MultiPoint'>

                            for pt in my_shape:
                                pt = self.move_point(pt, compare, dist)

'pt' becomes a 'Point':

POINT (13531245.47570414 2886003.268927813)

Then when I try to write it back to what it was (a Multipoint):

                    # write to outfile - WHERE EVERYTHING FAILS
                    output.write({'geometry':mapping(pt), 'properties': item['properties']})

I get this:

ValueError: Record's geometry type does not match collection schema's geometry type: 'Point' != 'MultiPoint'

If I try to re-write 'pt' to a 'MultiPoint' BEFORE the write operation (and again, this is a multipoint with ONE coordinate set) like so,

 mpoint = MultiPoint(pt)

I get this:

  File "C:\Python27\lib\site-packages\shapely\geometry\multipoint.py", line 143, in geos_multipoint_from_py
    assert len(array['shape']) == 2
AssertionError

*Full working code after help from @gene and @Vince *

def process_file(self, in_file, out_file, compare_file, dist):
    # open input file and compare file
    # loop over each
    with fiona.open(in_file, 'r') as input:
        meta = input.meta
        # The outFile has the same crs, schema as inFile
        with fiona.open(out_file, 'w', **meta) as output:
            with fiona.open(compare_file, 'r') as compare:
            # Read shapely geometries from file
            # Loop through all shapely objects
            # type(item) = 'dict'
                for item in input:
                    geom = item['geometry']
                    my_shape = shape(geom)

                    # check if it is a multipoint or point
                    if geom['type'] == 'MultiPoint':
                        # break Multipoints into points
                        # i.e. <class 'shapely.geometry.multipoint.MultiPoint'>

                        for pt in my_shape:
                            single_point = self.move_point(pt, compare, dist)
                            mpoint = MultiPoint([(single_point.x, single_point.y)])
                            mpoint_for_merge = shape(mapping(mpoint))
                            # write to outfile - WHERE EVERYTHING NOW WORKS
                            output.write({'geometry':mapping(mpoint_for_merge), 'properties': item['properties']})

                    elif geom['type'] == 'Point':
                        # return of move_point is a shapely geom
                        # i.e. <class 'shapely.geometry.point.Point'>
                        my_shape = self.move_point(my_shape, compare, dist)
                        # write to outfile
                        output.write({'geometry':mapping(my_shape), 'properties': item['properties']})

                    else:
                        raise ValueError('unhandled geometry type: ' + repr(geom.type))

Point and multipoint objects are distinct record types; shapefiles will only permit one or the other. The workaround is to make "multipoint" objects with a count of 1. You wouldn't want to do this if ALL features were single-point, since the envelope search algorithm isn't as effective as the point search protocol, but if you have a mixture you don't have a choice. — Vince, Mar 24 '15 at 17:39
Fiona enables to save MultiPoint geometries (there are transformed by the module) — gene, Mar 24 '15 at 17:58
@Vince I think you are right. I try to write a 1 pt multipoint though...and I get a error. No probs reading it. Error above./ — user14696, Mar 24 '15 at 20:14
Singleton multipoints are somewhat esoteric, so the framework may not have been tested for that case (I had 4 year old shape-writing source that had this wrong the whole time, so I know it could happen). I'm not familiar with Fiona, just the shapefile spec. — Vince, Mar 24 '15 at 21:12

gene · Accepted Answer · 2015-03-24T20:48:25.270

3

Keep things simple, a MultiPoint is a simple list of Points and the type of geometry can be obtained by Fiona (key ['geometry']['type']) of by Shapely (geometry.type)

with fiona.open(in_file, 'r') as input:
   with fiona.open(compare_file, 'r') as compare:
       meta = input.meta
       with fiona.open(out_file, 'w', **meta) as output:
          for item in input:
              geom = item['geometry']
              my_shape = shape(geom)
              # check if it is a Multipoint 
              if geom['type'] == 'MultiPoint':
                  for pt in my_shape:
                       my_shape = self.move_point(pt, compare, dist)
              else: # the Points
                   my_shape = self.move_point(my_shape, compare, dist)

You can also use

if my_shape.type == 'MultiPoint':

New

Your problem is to transform a Point into a MultiPoint. With Shapely you need to use the coordinates of the Point and not the Point itself (look at Collections of Points)

Multi = MultiPoint([(13531245.475704141, 2886003.2689278126)])
print Multi
MULTIPOINT (13531245.47570414 2886003.268927813)
point = Multi[0]
print point
POINT (13531245.47570414 2886003.268927813)

Now, you can't use mpoint = MultiPoint(point) but

mpoint = MultiPoint([(point.x, point.y)])
print mapping(mpoint)
{'type': 'MultiPoint', 'coordinates': ((13531245.475704141, 2886003.2689278126),)}
# and
print shape(mapping(mpoint))
MULTIPOINT (13531245.47570414 2886003.268927813)

edited Mar 24 '15 at 20:48

answered Mar 24 '15 at 17:52

gene

54,868
3
110
187

it is weird...but it is stopping at 'else' instead of 'if item['type'] == 'MultiPoint':' with this 'MULTIPOINT (13531245.47570414 2886003.268927813)'...this error...NotImplementedError: Multi-part geometries do not provide a coordinate sequence – user14696 Mar 24 '15 at 18:13
Sorry, look my correction and MultiLineString to separate individual lines using Python with GDAL/OGR, Fiona, Shapely – gene Mar 24 '15 at 18:19
I looked at this more (and thank you for the feedback on the line example). It seems like Fiona doesn't like me mashing a one-coordinate multipoint back into a multipoint....is this something you've seen before? I updated my question to show it. – user14696 Mar 24 '15 at 19:54
No, but why would you want that ? What is the schema of your shapefile ? (Points or MultiPoints ?). – gene Mar 24 '15 at 20:21
Shapefile type: MultiPoint ...is the in_file (which out_file is following the schema of). I ran cascaded_union ...got object of type 'Point' has no len().... :/ – user14696 Mar 24 '15 at 20:30
look above in New – gene Mar 24 '15 at 20:41

Using Geometry collections with Fiona to write output of Multipoints and Points?

1 Answers1

New