[GIS] Using Geometry collections with Fiona to write output of Multipoints and Points

fionapyqgisqgis-python-consoleshapely

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

Best Answer

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)
Related Question