Reports

For anyone who might have the same issue, here is how I solved it:
#1) Augment the outerplanar graph to a biconnected outerplanar graph:

def make_biconnected(graph):
    # Find cut vertices
    cut_vertices = list(nx.articulation_points(graph))

    for v in cut_vertices:
        neighbors = list(graph.neighbors(v))

        # Partition neighbors into blocks based on biconnected components
        subgraph = nx.Graph(graph)
        subgraph.remove_node(v)
        blocks = []
        for component in nx.connected_components(subgraph):
            block_neighbors = [n for n in neighbors if n in component]
            if block_neighbors:
                blocks.append(block_neighbors)

        # Add edges between blocks to make the graph biconnected
        for j in range(len(blocks) - 1):
            u = blocks[j][-1]
            w = blocks[j + 1][0]
            if not graph.has_edge(u, w):
                graph.add_edge(u, w)

                # If the edge breaks outerplanarity, revert and try other combinations
                if not is_outerplanar(graph):
                    graph.remove_edge(u, w)
                    for u_alt in blocks[j]:
                        for w_alt in blocks[j + 1]:
                            if not graph.has_edge(u_alt, w_alt):
                                graph.add_edge(u_alt, w_alt)
                                if is_outerplanar(graph):
                                    break
                                graph.remove_edge(u_alt, w_alt)

    return graph

#2) Find the largest simple cycle (this should be the Hamiltonian cycle or the outer face) and use the node ordering to create a circular layout

def generate_pos(graph):
    cycles = list(nx.simple_cycles(graph))
    maxLength = max( len(l) for l in cycles )
    path = list(l for l in cycles if len(l) == maxLength)
    pos = nx.circular_layout(path[0])
    return pos

#3) With the positions in #2 and the edges from the original outerplanar graph, use some math to sort the neighbor list of each node in ccw order and add the edges into an embedding.

def convert_to_outerplanar_embedding(graph, pos):
    # Step 1: Create a PlanarEmbedding object
    planar_embedding = nx.PlanarEmbedding()

    # Step 2: Assign the cyclic order of edges around each node based on positions
    for node in graph.nodes:
        neighbors = list(graph.neighbors(node))
        
        # Sort neighbors counterclockwise based on angle relative to the node's position
        neighbors.sort(key=lambda n: math.atan2(pos[n][1] - pos[node][1], pos[n][0] - pos[node][0]))

        # Step 3. Add edges in sorted order to the embedding
        planar_embedding.add_half_edge(node, neighbors[0])
        for i in range(1, len(neighbors)):
            u = neighbors[i-1]
            v = neighbors[i] 
            planar_embedding.add_half_edge(node, v, ccw = u)
    return planar_embedding

The result can be plotted using nx.draw_planar and should look like an outerplanar graph. I haven't tested it extensively, but it works so far for my use case.

79281382