Here is a the faster version:

point_geometries = [ee.Geometry.Point(lon, lat) for lon, lat in zip(df['Lon'], df['Lat'])]

collection = (ee.ImageCollection('COPERNICUS/S2_SR_HARMONIZED')

.filterDate(time_maps[region]['start'], time_maps[region]['end']))

mean_values = (collection

.select(BANDS)

.filterBounds(ee.FeatureCollection(point_geometries))

.mean()

.reduceRegions(collection=ee.FeatureCollection(point_geometries), reducer=ee.Reducer.mean(), scale=10)

)

# Extract the mean values for the specified bands

results = []

for feature in tqdm(mean_values.getInfo()['features']):

values = [feature['properties'][band] for band in BANDS]

results.append(values)

df[BANDS] = results