def merge(buildingIn, addressIn, mergedOut):
    addresses = []

    with collection(addressIn, "r") as input:
        for address in input:
            shape = asShape(address['geometry'])
            shape.original = address
            addresses.append(shape)

    # Load and index all buildings.
    buildings = []
    buildingShapes = []
    buildingIdx = index.Index()
    with collection(buildingIn, "r") as input:
        for building in input:
            shape = asShape(building['geometry'])
            building['properties']['addresses'] = []
            buildings.append(building)
            buildingShapes.append(shape)
            buildingIdx.add(len(buildings) - 1, shape.bounds)

    # Map addresses to buildings.
    for address in addresses:
        for i in buildingIdx.intersection(address.bounds):
            if buildingShapes[i].contains(address):
                buildings[i]['properties']['addresses'].append(
                    address.original)

    with open(mergedOut, 'w') as outFile:
	    outFile.writelines(json.dumps(buildings, indent=4))
	    print 'Exported ' + mergedOut

def generator_function(f, verbose):
    counter = 0
    for line in f:
        try:
            obj = geojson.loads(line)
        except:
            print "Unexpected error:", sys.exc_info()
            continue
        properties = property_map(obj.get('properties'))
        geometry = obj.get('geometry')
        geom_type = geometry.get('type')
        if geom_type == 'Polygon':
            poly = asShape(geometry)
            bounds = poly.bounds
            feature = geojson.Feature(id=counter, geometry=poly, properties=properties)

            print counter, bounds, properties.get('name')
            counter += 1
            yield (counter, bounds, json.loads(geojson.dumps(feature)))
        elif geom_type == 'MultiPolygon':
            mpoly = asShape(geometry)
            for poly in mpoly:
                bounds = poly.bounds
                feature = geojson.Feature(id=counter, geometry=poly, properties=properties)

                print counter, bounds, properties.get('name')
                counter += 1
                yield (counter, bounds, json.loads(geojson.dumps(feature)))
        else:
            print "unsupported type", geom_type
            continue

    def test_linestring_geojson(self):
        '''Create a line that goes from west to east (clip on)'''
        
        self.defineGeometry('LINESTRING')

        geom = LineString( [(-180, 32), (180, 32)] )

        self.insertTestRow(geom.wkt)

        # we should have a line that clips at 0...

        # for western hemisphere....
        tile_mimetype, tile_content = utils.request(self.config_file_content, "vector_test", "geojson", 0, 0, 0)
        self.assertEqual(tile_mimetype, "text/json")
        geojson_result = json.loads(tile_content)
        west_hemisphere_geometry = asShape(geojson_result['features'][0]['geometry'])
        expected_geometry = LineString([(-180, 32), (0, 32)])
        self.assertTrue(expected_geometry.almost_equals(west_hemisphere_geometry))

        # for eastern hemisphere....
        tile_mimetype, tile_content = utils.request(self.config_file_content, "vector_test", "geojson", 0, 1, 0)
        self.assertEqual(tile_mimetype, "text/json")
        geojson_result = json.loads(tile_content)
        east_hemisphere_geometry = asShape(geojson_result['features'][0]['geometry'])
        expected_geometry = LineString([(0, 32), (180, 32)])
        self.assertTrue(expected_geometry.almost_equals(east_hemisphere_geometry))

def intersection(shape, other):
    if not hasattr(shape,'__geo_interface__'): raise TypeError("%r does not appear to be a shape"%shape)
    if not hasattr(other,'__geo_interface__'): raise TypeError("%r does not appear to be a shape"%shape)
    o = geom.asShape(shape)
    o2 = geom.asShape(other)
    res = o.intersection(o2)
    return asShape(res)

def chunk(featureFileName, sectionFileName, pattern, key = None):

    # Load and index
    with collection(featureFileName, "r") as featureFile:
        featureIdx = index.Index()
        features = []
        for feature in featureFile:
            features.append(feature)
            featureIdx.add(len(features) - 1, asShape(feature['geometry']).bounds)

        # Break up by sections and export
        with collection(sectionFileName, "r") as sectionFile:
            i = 0
            for section in sectionFile:
                fileName = pattern % i
                if key:
                    fileName = pattern % section['properties'][key]
                    properties = {}
                    try:
                        with collection(fileName, 'w', 'ESRI Shapefile',
                                schema = featureFile.schema,
                                crs = featureFile.crs) as output:
                            sectionShape = asShape(section['geometry'])
                            for j in featureIdx.intersection(sectionShape.bounds):
                                if asShape(features[j]['geometry']).intersects(sectionShape):
                                    properties = features[j]['properties']
                                    output.write(features[j])
                            print "Exported %s" % fileName
                            i = i + 1
                    except ValueError:
                        print "Error exporting " + fileName
                        pprint(properties)
                        pprint(featureFile.schema)

def symmetric_difference(shape, other):
    if not hasattr(shape,'__geo_interface__'): raise TypeError("%r does not appear to be a shape"%shape)
    if not hasattr(other,'__geo_interface__'): raise TypeError("%r does not appear to be a shape"%shape)
    o = geom.asShape(shape)
    o2 = geom.asShape(other)
    res = o.symmetric_difference(o2)
    return asShape(res)

def shp_to_shply_multiply(shapefile_path):
    """
    Convert Polygon Shapefile to Shapely MultiPolygon
    :param shapefile_path: path to a shapefile on disk
    :return: shapely MultiPolygon
    """
    in_ply = shapefile.Reader(shapefile_path)

    # using pyshp reading geometry
    ply_shp = in_ply.shapes()

    if len(ply_shp) > 1:
        # using python list comprehension syntax
        # shapely asShape to convert to shapely geom
        ply_list = [asShape(feature) for feature in ply_shp]

        # create new shapely multipolygon
        out_multi_ply = MultiPolygon(ply_list)

        # # equivalent to the 2 lines above without using list comprehension
        # new_feature_list = []
        # for feature in features:
        #     temp = asShape(feature)
        #     new_feature_list.append(temp)
        # out_multi_ply = MultiPolygon(new_feature_list)

        print "converting to MultiPolygon: "
    else:
        print "one or no features found"
        shply_ply = asShape(ply_shp)
        out_multi_ply = MultiPolygon(shply_ply)

    return out_multi_ply

def mix_tma_g(air_tma, air_g):
    eps = 1e-4

    geom_tma = asShape(air_tma['geometry'])
    geom_g = asShape(air_g['geometry'])

    if geom_g.intersects(geom_tma):
        geom_tma_minus_g = geom_tma.difference(geom_g)
        geom_intersection = geom_g.intersection(geom_tma)
        # geom_g_minus_tma = geom_g.difference(geom_tma)

        if geom_intersection.area < eps:
            return

        if air_tma['properties']['lower'] >= air_g['properties']['upper']:
            return

        air_tma['geometry'] = geom_tma_minus_g

        air_intersection = Feature(
            geometry=geom_intersection,
            properties=deepcopy(air_tma['properties']))

        air_intersection['properties']['name'] = '{} - {}'.format(
            air_tma['properties']['name'], air_g['properties']['name'])

        air_intersection['properties']['lower'] = air_g['properties']['upper']
        del(air_intersection['properties']['lower_raw'])
        air_intersection['properties']['class'] = 'TMA_G_PG'
        air_intersection['properties']['fill'] = '#457aa3'

        return air_intersection

def reprojectPoint(obj, inProjection = 'epsg:3857', outProjection = 'epsg:4326'):
    
    #if our object is a raw string
    if type(obj) == str:
        obj = geojson.loads(obj)
    #if our object is a geojson dict
    if type(obj) == dict:
        #is our geojson a feature class?
        try:
            point = asShape(obj['geometry'])
        #or a regular geometry?
        except KeyError:
            point = asShape(obj)
    else:
        point = obj
        
    #pyproj transformation
    if type(inProjection) != str:
        inP = inProjection
    else:
        inP = pyproj.Proj(init=inProjection)
    if type(outProjection) != str:
        outP = outProjection
    else:
        outP = pyproj.Proj(init=outProjection)
        
    newCoords =  [pyproj.transform(inP,outP,i[0],i[1]) for i in list(point.coords)]
    return Point(newCoords) 

 def getNetworkArea(self, bufferDist):
     from shapely.geometry import asShape
     segs = self.getEdges()
     segUnion = asShape(segs[0]).buffer(bufferDist)
     for n in xrange(len(segs) - 1):
         print n
         segUnion = segUnion.union(asShape(segs[n + 1]).buffer(bufferDist))
     return segUnion.area

 def distance(self, b):
     try:
         g = IGeoreferenced(self.context)
         geom = asShape({'type': g.type, 'coordinates': g.coordinates})
         y0 = geom.centroid.y
         other = asShape(b.zgeo_geometry)
         d = geom.distance(other)
         return int(math.cos(math.pi*y0/180.0)*d/F/1000)
     except:
         log.warn("Failed to find distance between %s and %s" % (
             self.context, b.getPath()))
         raise

    def _add_some_objects(self, some_objects, dst):
        for an_object in some_objects:
            if asShape(an_object.geometry).within(self.buffer):
                the_distance = self.shape.project(asShape(an_object.geometry))

                # explanation: the buffer extends beyond the endpoints of the cross-section
                # points beyond the endpoints but within the buffer are
                # projected at 0. and length distance with a sharp angle
                # these points are not added to the cross-section
                # points exactly at 0. or length distance are also not added
                if (the_distance > 0.) and (the_distance < self.length):
                    dst.append((the_distance, an_object))

    def test_points_geojson(self):
        """
        Create 3 points (2 on west, 1 on east hemisphere) and retrieve as geojson.
        2 points should be returned in western hemisphere and 1 on eastern at zoom level 1
        (clip on)
        """

        self.defineGeometry("POINT")

        point_sf = Point(-122.42, 37.78)
        point_berlin = Point(13.41, 52.52)
        point_lima = Point(-77.03, 12.04)

        self.insertTestRow(point_sf.wkt, "San Francisco")
        self.insertTestRow(point_berlin.wkt, "Berlin")
        self.insertTestRow(point_lima.wkt, "Lima")

        ########
        # northwest quadrant should return San Francisco and Lima

        tile_mimetype, tile_content = utils.request(self.config_file_content, "vectile_test", "json", 0, 0, 1)
        geojson_result = json.loads(tile_content)

        self.assertTrue(tile_mimetype.endswith("/json"))
        self.assertEqual(geojson_result["type"], "FeatureCollection")
        self.assertEqual(len(geojson_result["features"]), 2)

        cities = []

        # Make sure that the right cities have been returned and that the geometries match

        for feature in geojson_result["features"]:
            if feature["properties"]["name"] == "San Francisco":
                cities.append(feature["properties"]["name"])
                self.assertTrue(point_sf.almost_equals(asShape(feature["geometry"])))

            elif feature["properties"]["name"] == "Lima":
                cities.append(feature["properties"]["name"])
                self.assertTrue(point_lima.almost_equals(asShape(feature["geometry"])))

        self.assertTrue("San Francisco" in cities)
        self.assertTrue("Lima" in cities)

        ##########
        # northeast quadrant should return Berlin

        tile_mimetype, tile_content = utils.request(self.config_file_content, "vectile_test", "json", 0, 1, 1)
        geojson_result = json.loads(tile_content)

        self.assertTrue(tile_mimetype.endswith("/json"))
        self.assertEqual(geojson_result["type"], "FeatureCollection")
        self.assertEqual(len(geojson_result["features"]), 1)
        self.assertTrue("Berlin" in geojson_result["features"][0]["properties"]["name"])

    def test_points_geojson(self):
        '''
        Create 3 points (2 on west, 1 on east hemisphere) and retrieve as geojson.
        2 points should be returned in western hemisphere and 1 on eastern at zoom level 0
        (clip on)
        '''
        
        self.defineGeometry('POINT')

        point_sf = Point(-122.4183, 37.7750)
        point_berlin = Point(13.4127, 52.5233)
        point_lima = Point(-77.0283, 12.0433)

        self.insertTestRow(point_sf.wkt, 'San Francisco')
        self.insertTestRow(point_berlin.wkt, 'Berlin')
        self.insertTestRow(point_lima.wkt, 'Lima')

        ########
        # western hemisphere should return San Francisco and Lima

        tile_mimetype, tile_content = utils.request(self.config_file_content, "vector_test", "geojson", 0, 0, 0)
        geojson_result = json.loads(tile_content)

        self.assertEqual(tile_mimetype, "text/json")
        self.assertEqual(geojson_result['type'], 'FeatureCollection')
        self.assertEqual(len(geojson_result['features']), 2)

        cities = []

        # Make sure that the right cities have been returned and that the geometries match

        for feature in geojson_result['features']:
            if feature['properties']['name'] == 'San Francisco':
                cities.append(feature['properties']['name'])
                self.assertTrue(point_sf.almost_equals(asShape(feature['geometry'])))

            elif feature['properties']['name'] == 'Lima':
                cities.append(feature['properties']['name'])
                self.assertTrue(point_lima.almost_equals(asShape(feature['geometry'])))

        self.assertTrue('San Francisco' in cities)
        self.assertTrue('Lima' in cities)

        ##########
        # eastern hemisphere should return Berlin

        tile_mimetype, tile_content = utils.request(self.config_file_content, "vector_test", "geojson", 0, 1, 0)
        geojson_result = json.loads(tile_content)

        self.assertEqual(tile_mimetype, "text/json")
        self.assertEqual(geojson_result['type'], 'FeatureCollection')
        self.assertEqual(len(geojson_result['features']), 1)
        self.assertTrue('Berlin' in geojson_result['features'][0]['properties']['name'])

def import_aquifers(self):
    aquifers = geojson.load(open(
        'src/iwlearn.project/iwlearn/project/dataimport/aquifers.json',
        'r'))
    parent = self.portal_url.getPortalObject()['iw-projects']['basins']['aquifers']
    for aquifer in aquifers['features']:
        ext = idn.normalize(aquifer['properties']['FIRST_ISAR'])
        new_obj_id = idn.normalize(aquifer['properties']['NAME']) + ext
        if new_obj_id in parent:
            mpoly = []
            new_obj=parent[new_obj_id]
            geo = IGeoManager(new_obj)
            add_geom = asShape(aquifer['geometry']).simplify(0.2)
            my_geom = wkt.loads(geo.wkt)
            if my_geom.geom_type == 'MultiPolygon':
                for poly in my_geom.geoms:
                    if poly.contains(add_geom):
                        continue
            elif my_geom.contains(add_geom):
                continue
            elif add_geom.contains(my_geom):
                q = add_geom.__geo_interface__
            else:
                if add_geom.geom_type == 'Polygon':
                    mpoly.append(add_geom)
                elif add_geom.geom_type == 'MultiPolygon':
                    mpoly += list(add_geom.geoms)
                if my_geom.geom_type == 'Polygon':
                    mpoly.append(my_geom)
                elif my_geom.geom_type == 'MultiPolygon':
                    mpoly += list(my_geom.geoms)
                q = MultiPolygon(mpoly).__geo_interface__
            geo.setCoordinates(q['type'], q['coordinates'])
            print new_obj_id, '*'
        else:
            self.portal_types.constructContent('Basin', parent, new_obj_id)
            new_obj=parent[new_obj_id]
            print new_obj_id
            new_obj.setTitle(aquifer['properties']['NAME'])
            new_obj.setDescription("Area: %s; Length: %s" % (
                            aquifer['properties']['Shape_Area'],
                            aquifer['properties']['Shape_Leng']))
            new_obj.setBasin_type('Aquifer')
            color='c1742c'
            style = IGeoCustomFeatureStyle(new_obj)
            style.geostyles.data['use_custom_styles']=True
            style.geostyles.data['polygoncolor']=color
            style.geostyles.update(style.geostyles)
            geo = IGeoManager(new_obj)

            q = asShape(aquifer['geometry']).simplify(0.2).__geo_interface__
            geo.setCoordinates(q['type'], q['coordinates'])

def reprojectLine(obj, inProjection = 'epsg:3857', outProjection = 'epsg:4326'):
    
    #if our object is a raw string
    if type(obj) == str:
        obj = geojson.loads(obj)
    #if our object is a geojson dict
    if type(obj) == dict:
        #is our geojson a feature class?
        try:
            line = asShape(obj['geometry'])
        #or a regular geometry?
        except KeyError:
            line = asShape(obj)

def read_input_layer(l):
    c = fiona.open(INFILE, 'r', driver='OpenFileGDB', layer=l)
    #-- store the geometries in a list (shapely objects)
    dPolys = {}
    for gid, each in enumerate(c):
        if MERGEFEATURES is True:
            k = each['properties']['TOP10_ID']
        else:
            k = gid
        if k not in dPolys:
            dPolys[k] = [asShape(each['geometry'])]
        else:
            dPolys[k].append(asShape(each['geometry']))
    return dPolys

def import_oceans(self):
    oceans = geojson.load(open(
        'src/iwlearn.project/iwlearn/project/dataimport/oceans.json',
        'r'))
    parent = self.portal_url.getPortalObject()['iw-projects']['basins']['oceans']
    for ocean in oceans['features']:
        mpoly = []
        geom = asShape(ocean['geometry'])
        for g in geom.geoms:
            if g.area > 5:
                mpoly.append(g)
        mp = MultiPolygon(mpoly).simplify(0.2)
        q = mp.__geo_interface__
        new_obj_id = idn.normalize(ocean['properties']['NAME'])
        print new_obj_id
        self.portal_types.constructContent('Basin', parent, new_obj_id)
        new_obj=parent[new_obj_id]
        new_obj.setTitle(ocean['properties']['NAME'])
        new_obj.setBasin_type('Ocean')
        color='aa22ff'
        style = IGeoCustomFeatureStyle(new_obj)
        style.geostyles.data['use_custom_styles']=True
        style.geostyles.data['polygoncolor']=color
        style.geostyles.update(style.geostyles)
        geo = IGeoManager(new_obj)
        geo.setCoordinates(q['type'], q['coordinates'])

        def check_geometry(r, feature, o):
            # we need both the "original" and "new" geometry to be
            # within the restriction area
            geom_attr, srid = self._get_geom_col_info(layer)
            geom_attr = getattr(o, geom_attr)
            geom = feature.geometry
            allowed = DBSession.query(func.count(RestrictionArea.id))
            allowed = allowed.join(RestrictionArea.roles)
            allowed = allowed.join(RestrictionArea.layers)
            allowed = allowed.filter(RestrictionArea.readwrite.is_(True))
            allowed = allowed.filter(Role.id == self.request.user.role.id)
            allowed = allowed.filter(Layer.id == layer.id)
            allowed = allowed.filter(or_(
                RestrictionArea.area.is_(None),
                RestrictionArea.area.ST_Contains(geom_attr)
            ))
            spatial_elt = None
            if geom and not isinstance(geom, geojson.geometry.Default):
                shape = asShape(geom)
                spatial_elt = from_shape(shape, srid=srid)
                allowed = allowed.filter(or_(
                    RestrictionArea.area.is_(None),
                    RestrictionArea.area.ST_Contains(spatial_elt)
                ))
            if allowed.scalar() == 0:
                raise HTTPForbidden()

            # check is geometry is valid
            self._validate_geometry(spatial_elt)

    def read_one(self):
        set_common_headers(self.request, "layers", NO_CACHE, add_cors=True)

        layer = self._get_layer_for_request()
        protocol = self._get_protocol_for_layer(layer)
        feature_id = self.request.matchdict.get("feature_id", None)
        feature = protocol.read(self.request, id=feature_id)
        if not isinstance(feature, Feature):
            return feature
        if layer.public:
            return feature
        if self.request.user is None:
            raise HTTPForbidden()
        geom = feature.geometry
        if not geom or isinstance(geom, geojson.geometry.Default):  # pragma: no cover
            return feature
        shape = asShape(geom)
        srid = self._get_geom_col_info(layer)[1]
        spatial_elt = from_shape(shape, srid=srid)
        allowed = DBSession.query(func.count(RestrictionArea.id))
        allowed = allowed.join(RestrictionArea.roles)
        allowed = allowed.join(RestrictionArea.layers)
        allowed = allowed.filter(Role.id == self.request.user.role.id)
        allowed = allowed.filter(Layer.id == layer.id)
        allowed = allowed.filter(or_(
            RestrictionArea.area.is_(None),
            RestrictionArea.area.ST_Contains(spatial_elt)
        ))
        if allowed.scalar() == 0:
            raise HTTPForbidden()

        return feature