How do you select custom colours to fill regions of a Voronoi diagram using scipy.spatial's Voronoi package?

I have a code that reads in a csv file containing places, their latitude and longitude, and whether I've been there (marked by 'Y' or 'N'). I want to present this in a Voronoi diagram, but crucially I want to colour the regions according to the following criteria:

• If I've not been there ('N') mark region as white
• If I have been there ('Y') mark region as a random colour that is unique to any other colours on the map

I have two issues with my code: I can't seem to match the regions produced by Voronoi(coords) to my original list of places (so I can't readily correspond the region to the colour I want to fill that region, because the order in which the regions are defined and filled is not the order in which they appear in the dataframes), and secondly, how to ensure each random colour is "unique", i.e. not the same random colour appearing in each region.

The second issue I am fairly certain I can find a solution for, it's the first issue that I am struggling with. Any help would be much appreciated. Here is my code:

import pandas as pd
from scipy.spatial import Voronoi, voronoi_plot_2d
import numpy as np

# read places, with lat and lon
places = pd.read_csv("places.csv")

# convert subset to numpy array
coords = places[['Longitude','Latitude']].values

# add 4 distant dummy points
coords = np.append(coords, [[999,999], [-999,999], [999,-999], [-999,-999]], axis = 0)
colours = np.append(places[['Been']].values, [['N'], ['N'], ['N'], ['N']], axis = 0).flatten()
colours = np.insert(colours, 0, ['N','N'])

# assign a random colour to the array if visited, leave white if not
colours[colours == 'N'] = 'w'
import random
r = lambda: random.randint(0,255)
colours[colours == 'Y'] = '#%02X%02X%02X' % (r(),r(),r())
print(places['Place'])
print(colours)
print(coords)
# compute voronoi tesselation
vor = Voronoi(coords)

# plot voronoi diagram
import matplotlib.pyplot as plt
fig = voronoi_plot_2d(vor, show_vertices = False)

j = -1
# colourise the regions
for region in vor.regions:
    j = j+1
    print(region)
    if not -1 in region:
        polygon = [vor.vertices[i] for i in region]
        plt.fill(*zip(*polygon), colours[j])

# fix the range of axes, plot locations
plt.plot(coords[:,0], coords[:,1], 'ko')
plt.xlim([places['Longitude'].min() - 0.6, places['Longitude'].max() + 0.6]), plt.ylim([places['Latitude'].min() - 0.6, places['Latitude'].max() + 0.6])

# annotate each point with the place name
[plt.annotate(places['Place'][i], (coords[i,0], coords[i,1]), xytext=(coords[i,0]-0.2, coords[i,1]+0.2)) for i in range(len(places))]
plt.show()

I've produced this in Jupyter Notebooks. My csv file looks like this:

Place,Latitude,Longitude,Been
Bern,46.948,7.4474,N
Juras,47.0086,6.7856,N
Lake Como,46.016,9.2572,N
Lyon,45.764,4.8357,N
Marseille,43.2965,5.3698,N
Milan,45.4642,9.19,N
Monaco,43.7384,7.4246,N
Mont Blanc,45.8326,6.8652,N
Mont Saleve,46.0942,6.1403,Y
Munich,48.1351,11.582,N
Turin,45.0703,7.6869,N
Zurich,47.3769,8.5417,N

I've left in print statements, if that helps, and some commented lines of code from my tests. I currently colour the region marked by Lyon, which demonstrates the issue I have of not knowing how to use region in vor.regions to retrieve the coordinates it corresponds to. Thanks in advance.