I tried applying a Genetic Algorithm to the problem above. I made an initial guess that two additional nodes would connect all three disconnected components.

import networkx as nx
import pygad
import math
from libpysal import weights
import numpy as np

no_comps_target = 1                    # a connected graph has 1 component
max_dist_between_nodes = math.sqrt(2)  # in km
num_pts = 2                            # number of additional nodes to add
num_genes = num_pts*2                  # number of genes required with the GA. A gene each for x and y coordinates.

# Generate coordinate np array of existing components
centroids = np.array([(v) for k, v in pos.items()])

# Create indcies for new nodes within the GA solution list
y_ix = [x for x in range(num_genes) if x % 2 != 0] 
x_ix = [x for x in range(num_genes) if x % 2 == 0] 

# Define fitness function
def my_fitness_func(solution, solution_idx):
    
    # Select coordinates of GA solution
    xs = np.array(solution[x_ix])
    ys = np.array(solution[y_ix])
    new_pts = np.column_stack((xs,ys))
   
    # Create one set for all coordinates
    all_pts = np.append(centroids, new_pts,axis=0)
        
    # Calculate weights using a distance band equal to the max distance between nodes
    w = weights.DistanceBand.from_array(all_pts, 
                                        threshold=max_dist_between_nodes, 
                                        silence_warnings=True)

    # Convert to a networkx obejct
    G = w.to_networkx()
    
    # Calculate the number of graph components for G
    # Target is 1 - fully connected graph
    no_comps_solution = nx.number_connected_components(G)
         
    
    # Calculate solution fitness
    fitness = 1.0 / np.abs(no_comps_target - no_comps_solution + 0.000001)
        
        return fitness
# Set constraints on possible solution locations
x_max = 5
x_min = 0
y_max = 5
y_min = 1

ga_instance = pygad.GA(
    num_generations=20,
    num_parents_mating=2,
    fitness_func=my_fitness_func,
    sol_per_pop=10,
    num_genes=num_genes,
    gene_type= int,
    gene_space= [{'low': x_min, 'high': x_max, 'step': 1}, 
                 {'low': y_min, 'high': y_max, 'step': 1}] * num_pts,
    mutation_num_genes=1,
    parent_selection_type = "sss",
    keep_parents =2,
    stop_criteria="saturate_10" # Stop if no progress after 10 generations
)

ga_instance.run()

# If final reached a maximum we should expect a fitness of 100,000.
solution, solution_fitness, solution_idx = ga_instance.best_solution()
print("Parameters of the best solution : {solution}".format(solution=solution))
print("Fitness value of the best solution = {solution_fitness}".format(solution_fitness=solution_fitness))

This gives a valid solution:

Parameters of the best solution : [3 3 2 4]
Fitness value of the best solution = 1000000.0

And from running it multiple times, I get multiple valid solutions. Which I think makes sense. Also, how to determine the optimal number of additional nodes? Especially if this problem were to be much larger. I'd still like to know if there are other ways of solving this problem. Especially if they come with less code!

I am quite convinced that this problem is NP-hard. The closest problem I know is the geometric Steiner tree problem with octilinear metric. I have two, rather quick-and-dirty, suggestions. Both are heuristic.

1st idea: Formulate the problem as an Euclidean Steiner tree problem (https://en.wikipedia.org/wiki/Steiner_tree_problem#Euclidean_Steiner_tree), where you consider just the nodes of your problem and forget about the edges at first. Solve the problem by using GeoSteiner: http://www.geosteiner.com/ This should quickly give you a solution for problems with 10000 or more nodes (if you need to solve bigger problems, you can write the problem out with GeoSteiner after the full-Steiner tree generation and use https://scipjack.zib.de/). There is no Python interface, but just write your problem to a plain text file, the syntax is quite easy. Afterward, put additional nodes into the solution provided by GeoSteiner such that the \sqrt(2) condition is satisfied. Finally, you need to do some clean-up to get rid of redundant edges, because the solution will not take into account that you already have edges in your original problem. Take all the edges and nodes that you have computed so far and define a weighted graph by giving all of your original edges weight 0 and all of the newly added edges weight 1. Consider a Steiner tree problem in graphs (https://en.wikipedia.org/wiki/Steiner_tree_problem#Steiner_tree_in_graphs_and_variants) on this weighted graph, where the terminal set corresponds to your original nodes. Solve this Steiner tree problem with SCIP-Jack: https://scipjack.zib.de/.

2nd idea: Consider your problem directly as a Steiner tree problem in graphs as follows: Each of the original edges is assigned weight 0, consider all original nodes as terminals. Add additional nodes and edges at distance at most \sqrt(2) from each other. For example, you could put a big rectangle around all your connected components and from each node add recursively additional 8 nodes in an angle at degrees 0,45,90,... at a distance of sqrt(2) and with edge of weight 1 in the Steiner tree problem in graphs, as long as they are within the rectangle. If one of these nodes is within distance sqrt(2) of one of your original nodes, connect them directly with an edge of weight 1. Solve the corresponding Steiner tree problem in graphs with SCIP-Jack.