Trouble visualizing gridded Copernicus Marine data on a map after transforming it to a GeoDataFrame

Question

I'm working with copernicus marine data for the school and I have a problem:

I took my inspiration from this code.

We have copernicus marine service data in netcdf4 format, these data have as variables environmental data, depth, time, latitude and longitude. I transform these ".nc" files into a dataframe and then select only one depth. Once all the data has been processed I transform the dataframe into a geodataframe, to obtain geometry. Once we've obtained the geodataframe, we want to generate a grid of the space under study, based on the coordinates present in the geodataframe. We then merge the geodataframe data (index and coordinates) with the grid.

We then want to obtain a grid with the corresponding parameter value for each cell in the grid. Finally, we want to visualize the results. The problem is that when I try to visualize after having done the steps described above, I get a map with a grid on top that overflows onto the ground.

I think the problem comes from the transformation in geodataframe and the geometry point. Because I try to plot the gdf and it gives me bad results.

This is my code:

import matplotlib.pyplot as plt
import contextily as cx



# Open the NetCDF file for reading
fn = 'file.nc'
ds_2 = nc.Dataset(fn)
ds_2.variables.keys()


var = ds_2.variables['var']

# Get the dimension coordinates

time = ds_2.variables['time'][:]
lat = ds_2.variables['latitude'][:]
lon = ds_2.variables['longitude'][:]
depth = ds_2.variables['depth'][:]

# Create a meshgrid of time, lat, and lon coordinates
time_mesh, lat_mesh, lon_mesh , depth_mesh= np.meshgrid(time, lat, lon,depth, indexing='ij')

# Reshape the arrays into 1D arrays
time_1d = time_mesh.flatten()
lat_1d = lat_mesh.flatten()
lon_1d = lon_mesh.flatten()
depth_1d = depth_mesh.flatten()
var_1d = var[:].flatten()
# Create a pandas DataFrame

data = {
    'time': time_1d,
    'latitude': lat_1d,
    'longitude': lon_1d,
    'depth': depth_1d,
    'var':var_1d
}

df = pd.DataFrame(data)

def generate_grid(total_bounds, n_cells, cell_size):

    xmin, ymin, xmax, ymax = total_bounds
    crs = 3857
    grid_cells = []

    for x0 in np.arange(xmin, xmax+cell_size, cell_size):
        for y0 in np.arange(ymin, ymax+cell_size, cell_size):
            # bounds
            x1 = x0-cell_size
            y1 = y0+cell_size
            grid_cells.append(shapely.geometry.box(x0, y0, x1, y1))
    cell = gpd.GeoDataFrame(grid_cells, columns=['geometry'], crs=crs)
    return cell

gdf = gpd.GeoDataFrame(df, geometry=gpd.points_from_xy(df['longitude'], df['latitude']), crs='EPSG:4326')

gdf['geometry'].apply(lambda x: WKTElement(x.wkt, srid=4326))

df_wm =gdf.to_crs(epsg=3857)

df_wm.plot()
plt.show()