def main(infile, outfile, driver):

    with fio.open(infile) as src:
        meta = src.meta
        meta['driver'] = driver

    with fio.open(infile) as src, fio.open(outfile, 'w', **meta) as dst:
        with click.progressbar(src) as features:
            for feat in features:

                east = deepcopy(feat)
                west = deepcopy(feat)

                east_geom = shape(east['geometry'])
                west_geom = shape(west['geometry'])

                # if 'Point' not in asShape(feat['geometry']).type:
                #     east_geom = east_geom.simplify(0.0001).buffer(0)
                #     west_geom = west_geom.simplify(0.0001).buffer(0)

                east_geom = translate(east_geom, xoff=180)
                west_geom = translate(west_geom, xoff=-180)

                if not east_geom.is_empty:
                    east['geometry'] = mapping(east_geom)
                    dst.write(east)

                if not west_geom.is_empty:
                    west['geometry'] = mapping(west_geom)
                    dst.write(west)

def map_to_polygon(shape,origin,position,rotation,size):
    geom_obj=Polygon(shape)
    geom_obj=affinity.translate(geom_obj, -origin[0],-origin[1],0)
    geom_obj=affinity.scale(geom_obj,size[0],size[1],origin=(0,0))
    geom_obj=affinity.rotate(geom_obj,int(rotation*360),origin=(0,0))
    geom_obj=affinity.translate(geom_obj, position[0],position[1],0)
    return geom_obj

    def _refine_predug(self, candidate):
        """ uses the color information directly to specify the predug """
        # determine a bounding rectangle
        region = candidate.bounds
        size = max(region.width, region.height)
        region.buffer(0.5 * size)  # < increase by 50% in each direction

        # extract the region from the image
        slice_x, slice_y = region.slices
        img = self.image[slice_y, slice_x].astype(np.uint8, copy=True)

        # build the estimate polygon
        poly_p = affinity.translate(candidate.polygon, -region.x, -region.y)
        poly_s = affinity.translate(self.ground.get_sky_polygon(), -region.x, -region.y)

        def fill_mask(color, margin=0):
            """ fills the mask with the buffered regions """
            for poly in (poly_p, poly_s):
                pts = np.array(poly.buffer(margin).boundary.coords, np.int32)
                cv2.fillPoly(mask, [pts], color)

        # prepare the mask for the grabCut algorithm
        burrow_width = self.params["burrows/width"]
        mask = np.full_like(img, cv2.GC_BGD, dtype=np.uint8)  # < sure background
        fill_mask(cv2.GC_PR_BGD, 0.25 * burrow_width)  # < possible background
        fill_mask(cv2.GC_PR_FGD, 0)  # < possible foreground
        fill_mask(cv2.GC_FGD, -0.25 * burrow_width)  # < sure foreground

        # run GrabCut algorithm
        img = cv2.cvtColor(img, cv2.COLOR_GRAY2RGB)
        bgdmodel = np.zeros((1, 65), np.float64)
        fgdmodel = np.zeros((1, 65), np.float64)
        try:
            cv2.grabCut(img, mask, (0, 0, 1, 1), bgdmodel, fgdmodel, 2, cv2.GC_INIT_WITH_MASK)
        except:
            # any error in the GrabCut algorithm makes the whole function useless
            logging.warn("GrabCut algorithm failed for predug")
            return candidate

        # turn the sky into background
        pts = np.array(poly_s.boundary.coords, np.int32)
        cv2.fillPoly(mask, [pts], cv2.GC_BGD)

        # extract a binary mask determining the predug
        predug_mask = (mask == cv2.GC_FGD) | (mask == cv2.GC_PR_FGD)
        predug_mask = predug_mask.astype(np.uint8)

        # simplify the mask using binary operations
        w = int(0.5 * burrow_width)
        kernel = cv2.getStructuringElement(cv2.MORPH_ELLIPSE, (w, w))
        predug_mask = cv2.morphologyEx(predug_mask, cv2.MORPH_OPEN, kernel)

        # extract the outline of the predug
        contour = regions.get_contour_from_largest_region(predug_mask)

        # translate curves back into global coordinate system
        contour = curves.translate_points(contour, region.x, region.y)

        return shapes.Polygon(contour)

def crop_line(line,pts_array):
    line_start,line_end=list(line.coords)
    line_start=Point(line_start); line_end=Point(line_end)
    x=(line_start.x-line_end.x)/line.length
    y=(line_start.y-line_end.y)/line.length
    center=line.centroid
    pts=pts_array-np.array([center.x,center.y])
    pts=pts[:,0]*x+pts[:,1]*y
    return pts2line([translate(center,xoff=(np.min(pts)*x), yoff=np.min(pts)*y),translate(center,xoff=(np.max(pts)*x), yoff=np.max(pts)*y)])

def fit_shape(shape, width, height, padding=0):
    width -= padding * 2
    height -= padding * 2
    x1, y1, x2, y2 = shape.bounds
    w, h = x2 - x1, y2 - y1
    s = min(width / w, height / h)
    shape = translate(shape, -x1, -y1)
    shape = scale(shape, s, s, origin=(0, 0, 0))
    shape = translate(shape, padding, padding)
    return shape

def findIntersect( ):
    time, max_area, max_time = 0.0, 0.0, 0.0
    while time < 5:
        time+=.001
        temp = asteroid1.asteroid.intersection(asteroid2.asteroid).area
        if temp > max_area:
            max_area = temp
            max_time = time
        asteroid1.asteroid = translate(asteroid1.asteroid, asteroid1.velocity[0]*.001, asteroid1.velocity[1]*.001)
        asteroid2.asteroid = translate(asteroid2.asteroid, asteroid2.velocity[0]*.001, asteroid2.velocity[1]*.001)
    return max_time

 def wrap_americas(self, wrapping_point=-50):
     new_mp = []
     try:
         for i, p in enumerate(self.patch):
             if min(p.envelope.exterior.xy[0]) < wrapping_point or\
                     max(p.envelope.exterior.xy[0]) < wrapping_point:
                 new_mp.append(translate(p, xoff=360.))
             else:
                 new_mp.append(p)
         self.patch = MultiPolygon(new_mp)
     except TypeError:
         if self.get_maxx() < wrapping_point or\
                 self.get_minx() < wrapping_point:
             self.patch = translate(self.patch, xoff=360.)

def clip_geometry_to_srs_bounds(geometry, pyramid, multipart=False):
    """
    Clip input geometry to SRS bounds of given TilePyramid.

    If geometry passes the antimeridian, it will be split up in a multipart
    geometry and shifted to within the SRS boundaries.
    Note: geometry SRS must be the TilePyramid SRS!

    - geometry: any shapely geometry
    - pyramid: a TilePyramid object
    - multipart: return list of geometries instead of a GeometryCollection
    """
    try:
        assert geometry.is_valid
    except AssertionError:
        raise ValueError("invalid geometry given")
    try:
        assert isinstance(pyramid, TilePyramid or MetaTilePyramid)
    except AssertionError:
        raise ValueError("not a TilePyramid object")
    pyramid_bbox = box(
        pyramid.left, pyramid.bottom, pyramid.right, pyramid.top)

    # Special case for global tile pyramids if geometry extends over tile
    # pyramid boundaries (such as the antimeridian).
    if pyramid.is_global and not geometry.within(pyramid_bbox):
        inside_geom = geometry.intersection(pyramid_bbox)
        outside_geom = geometry.difference(pyramid_bbox)
        # shift outside geometry so it lies within SRS bounds
        if isinstance(outside_geom, Polygon):
            outside_geom = [outside_geom]
        all_geoms = [inside_geom]
        for geom in outside_geom:
            geom_left = geom.bounds[0]
            geom_right = geom.bounds[2]
            if geom_left < pyramid.left:
                geom = translate(geom, xoff=2*pyramid.right)
            elif geom_right > pyramid.right:
                geom = translate(geom, xoff=-2*pyramid.right)
            all_geoms.append(geom)
        if multipart:
            return all_geoms
        else:
            return GeometryCollection(all_geoms)

    else:
        if multipart:
            return [geometry]
        else:
            return geometry

def load_letters(letters):
    x = 0.0
    s = 1.82748538
    p = 0.125
    polygons = []
    for letter in letters:
        polygon = load_letter(letter)
        polygon = scale(polygon, s, s)
        x1, y1, x2, y2 = polygon.bounds
        polygon = translate(polygon, -x1, -y1)
        polygon = translate(polygon, x)
        x += (x2 - x1) + p
        polygons.append(polygon)
        print polygon.bounds
    return MultiPolygon(polygons)

    def tile(self, tile_builder, offset, max_width, max_height, polygon_mask=None):
        """Extract a tile from the image

        Parameters
        ----------
        tile_builder: TileBuilder
            A tile builder for constructing the Tile object
        offset: (int, int)
            The (x, y) coordinates of the pixel at the origin point of the tile in the parent image
        max_width: int
            The maximum width of the tile
        max_height: int
            The maximum height of the tile
        polygon_mask: Polygon (optional, default: None)
            The polygon representing the alpha mask to apply to the image window. The polygon must be referenced
            to the full image top-left pixel.

        Returns
        -------
        tile: Tile
            The extracted tile

        Raises
        ------
        IndexError: if the offset is not inside the image
        TileExtractionException: if the tile cannot be extracted
        """
        if not self._check_tile_offset(offset):
            raise IndexError("Offset {} is out of the image.".format(offset))
        width = min(max_width, self.width - offset[0])
        height = min(max_height, self.height - offset[1])
        translated_polygon = translate(polygon_mask, -offset[0], -offset[1]) if polygon_mask is not None else None
        return tile_builder.build(self, offset, width, height, polygon_mask=translated_polygon)

def pack(items):
    shifted = [ ]
    for i, item in enumerate(items):
        shifted.append(
            affinity.translate(item, (i%6)*15, (i//6)*15)
            )
    return ops.cascaded_union(shifted)

def module_third(wafer_size, ncells, module_center=Point((0,0)), module_id=0, triggercell_size=1):
    # Compute the cell grid length for 1/3 of a module
    grid_size = int(math.sqrt(ncells/3))
    # This geometry is of the rotated type
    # The wafer size below is the vertex to vertex distance
    # The cell size is the edge to edge distance
    cell_size = wafer_size/grid_size/2.
    # Create grid of cells along the usual hexagon axes (60deg rotated axes)
    grid_generator = GridGenerator('diamond', grid_size)
    reference_position = translate(module_center,
            xoff=-cell_size*(grid_size-1), 
            yoff=cell_size*tan30)
    cell_centers = grid_generator(reference_position, cell_size)
    # Create cells corresponding to 1/3 of a module
    hex_generator = HexagonGenerator(cell_size*tan30)
    cell_transform = CellTransform(cell_size, grid_size)
    cell_vertices = [cell_transform(i)(hex_generator(point)) for i,point in enumerate(cell_centers)]
    # Merge cells in trigger cells if requested
    if triggercell_size>1:
        cell_vertices = trigger_cells(cell_vertices, size=triggercell_size)
        cell_centers = [c.centroid for c in cell_vertices]
    cells = []
    for i,(vertices,center) in enumerate(zip(cell_vertices,cell_centers)):
        cells.append(Cell(
            id=compute_id(module=module_id, third=0, cell=i),
            layer=1,
            zside=1,
            subdet=3,
            module=module_id,
            center=center,
            vertices=vertices
            ))
    return cells

def gen_feature(feature_name, model_number, image_path, shape, transl, scaling, scn):
    print("Called gen_feature @ %s" % feature_name)
    # let's first translate our feature.
    ip = Polygon(list(shape))  # map(lambda x: (x[0], 4x[1]), shape)))
    p = translate(ip, xoff=transl[0], yoff=transl[1])
    if(feature_name == "Mountain"):
        center_z = 0
        center_pos = p.centroid.coords[0]
        rd = int((max([dist(x, center_pos) for x in p.exterior.coords]) / 2) ** 0.5)
        print("Radius = %d" % rd)
        print("Center = %d, %d" % (center_pos[0], center_pos[1]))
        return Mountain(rd, center_z, center_pos)
    elif(feature_name == "MountainImg"):
        center_z = 0
        center_pos = p.bounds[0:2]
        print("Center = %d, %d" % (center_pos[0], center_pos[1]))
        return MountainImg(p, center=center_pos)
    elif(feature_name == "Roads"):
        pass
    elif(feature_name == "Image"):
        f = ImageFeature(image_path)
        f.shape = p
        return f
    elif(feature_name == "Vegetation"):
        for a in p.exterior.coords:
            print(a)
        return Vegetation(p, model=AbstractModel(scn["model_path"], 0.02, (0, 0)), tree_number=model_number)
    elif(feature_name == "Urban"):
        pass
    elif(feature_name == "WaterArea"):
        pass
    elif(feature_name == "River"):
        pass

    def window(self, offset, max_width, max_height, polygon_mask=None):
        """Build an image object representing a window of the image

        Parameters
        ----------
        offset: (int, int)
            The (x, y) coordinates of the pixel at the origin point of the window in the parent image
        max_width: int
            The maximum width of the window
        max_height: int
            The maximum height of the window
        polygon_mask: Polygon (optional, default: None)
            The polygon representing the alpha mask to apply to the image window. The polygon must be referenced
            to the full image top-left pixel.

        Returns
        -------
        window: ImageWindow
            The resulting image window
        """
        # width are bound to the current window size, not the parent one
        width = min(max_width, self.width - offset[0])
        height = min(max_height, self.height - offset[1])
        translated_polygon = translate(polygon_mask, -offset[0], -offset[1]) if polygon_mask is not None else None
        return ImageWindow(self, offset, width, height, polygon_mask=translated_polygon)

 def test_translate(self):
     ls = load_wkt('LINESTRING(240 400 10, 240 300 30, 300 300 20)')
     # test default offset of 0.0
     tls = affinity.translate(ls)
     self.assertTrue(tls.equals(ls))
     # test all offsets
     tls = affinity.translate(ls, 100, 400, -10)
     els = load_wkt('LINESTRING(340 800 0, 340 700 20, 400 700 10)')
     self.assertTrue(tls.equals(els))
     # Do explicit 3D check of coordinate values
     for a, b in zip(tls.coords, els.coords):
         for ap, bp in zip(a, b):
             self.assertEqual(ap, bp)
     # retest with named parameters for the same result
     tls = affinity.translate(geom=ls, xoff=100, yoff=400, zoff=-10)
     self.assertTrue(tls.equals(els))

    def valid_edge(self, state, primitive, plot = True):

        bounding_poly = affinity.rotate(primitive.bounding_poly, state[2], origin = (0.0, 0.0), use_radians = True)
        bounding_poly = affinity.translate(bounding_poly, state[0], state[1])

        #Drawing Primitive-TAKE OUT TODO
        if plot:
            if bounding_poly.intersects(self.world_polys):
                color = 'r'
            else:
                color = 'b'

            fig = plt.figure(2)
            fig.clear()
            plt.ion()
            ax = fig.add_subplot(111, aspect = 'equal')
            for poly in self.world_polys:
                P = PolygonPatch(poly, fc = 'k', zorder = 2)
                ax.add_patch(P)
            P = PolygonPatch(bounding_poly, fc = color, ec = '#999999', alpha = 1, zorder = 1)
            polyPatch = ax.add_patch(P)
            ax.set_xlim(0,50)
            ax.set_ylim(0,50)
            fig.show()
            plt.pause(0.1)
            #pdb.set_trace()
        
        if bounding_poly.intersects(self.world_polys):

            return False
        return True

def saveR(rectangles, A, cc, n):
    zona = takeZona(n)
    polis =[]
    for r in rectangles:
        polis.append(r[0])

    union = affinity.rotate(cascaded_union(polis), -A, origin=cc)
    dx = union.centroid.x-cc.x
    dy = union.centroid.y-cc.y
    print 'translate : ',dx, dy
    data2save=()
    for r in rectangles:
        rotated=affinity.rotate(r[0], -A, origin=cc)
        translated = affinity.translate(rotated, -dx, -dy)
        #verificar si interseca
        print zona.intersects(translated)
        if zona.intersects(translated):
            data = (n, A, r[1], r[2], "SRID=25831;"+str(translated))
            data2save += (data,)
        #print data
    conn = db.connect(rootData)
    c = conn.cursor()
    c.executemany('''insert into elementsdiv ( name, ang, x, y, geometry ) values ( ?, ?, ?,?, GeomFromEWKT( ? ) )''', data2save )
    conn.commit()
    conn.close()

    return

    def _load_shapefile(self, shp, index_field, convert_coordinates, remove_offset, simplify):

        df = shp2df(shp)

        if index_field is not None:
            df.index = df[index_field]

        proj4 = get_proj4(shp)

        if proj4 != self.proj4:
            df['geometry'] = projectdf(df, proj4, self.proj4)

        # convert projected coordinate units and/or get rid z values if the shapefile has them
        if convert_coordinates != 1 or df.iloc[0]['geometry'].has_z:
            df['geometry'] = [transform(lambda x, y, z=None: (x * convert_coordinates,
                                                              y * convert_coordinates), g)
                              for g in df.geometry]

        # remove model offset from projected coordinates (llcorner = 0,0)
        if remove_offset:
            df['geometry'] = [translate(g,
                                        -1 * self.extent_proj[0],
                                        -1 * self.extent_proj[1]) for g in df.geometry]

        if simplify > 0:
            df['geometry'] = [g.simplify(simplify) for g in df.geometry]
        return df

 def make(self):
     # Create new geometry
     dx = self.destination[0] - self.origin[0]
     dy = self.destination[1] - self.origin[1]
     self.geometry = [DrawToolShape(affinity.translate(geom.geo, xoff=dx, yoff=dy))
                      for geom in self.draw_app.get_selected()]
     self.complete = True

def main(opts):
    pattern = loads(open(opts.input, "r").read())
    extent = loads(open(opts.extent, "r").read())

    if not contains.matches(extent.relate(pattern)):
        print "ERROR: pattern must be contained within the extent"
        return

    c = pattern.centroid
    (xs, ys) = extent.boundary.xy
    (minx, maxx, miny, maxy) = (min(xs) - c.x, max(xs) - c.x, min(ys) - c.y, max(ys) - c.y)

    outputFile = open(opts.output, "w")

    geoms = []

    while len(geoms) < opts.number:
        dx = random.uniform(minx, maxx)
        dy = random.uniform(miny, maxy)

        geom = translate(pattern, xoff=dx, yoff=dy)

        if contains.matches(extent.relate(geom)):
            # Check that it is within the extent
            overlap = False
            for g in geoms:
                if intersects.matches(g.relate(geom)):
                    overlap = True
            if overlap == False:
                geoms.append(geom)

    for geom in geoms:
        outputFile.write(dumps(geom) + "\n")
    outputFile.close()