st_transform in R: coordinates after transformation from EPSG 3857 to 3035 do not match location from before

I am trying to transform some coordinates from WGS84 / Pseudo-Mercator (EPSG 3857) to ETRS89 / ETRS-LAEA (INSPIRE; EPSG 3035), but the location after the transformation does not match the location before.

Both example locations are in Germany, the output coordinates are not. Any idea what could be wrong here?

#create example df
lat <- c("48.109", "52.221")
long <- c("9.212", "13.338")
df <- data.frame(lat, long)
df$coords <- data.frame(longitude = df$long, latitude = df$lat)

df$coords_sf <- st_as_sf(df$coords, coords = c("longitude", "latitude"), crs = 3857)

df$transformed_coords_sf <- st_transform(df$coords_sf, 3035)

This line:

df$coords_sf <- st_as_sf(df$coords, coords = c("longitude", "latitude"), crs = 3857)

Is not doing what you think is doing. Actually, you used WGS84 (i.e. EPSG:4326) when defining your lat/long objects, so you should use st_as_sf(df$coords, coords = c("longitude", "latitude"), crs = 4326) in first place and then transform to whatever EPSG you want.

Check this example:

library(sf)
#> Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1; sf_use_s2() is TRUE

# Proof with map
library(ggplot2)
deu <- giscoR::gisco_get_countries(country = "Germany")
deu_gg <- ggplot(deu) +
  geom_sf()
# Map stored

# create example df
lat <- c("48.109", "52.221")
long <- c("9.212", "13.338")
df <- data.frame(lat, long)
df$coords <- data.frame(longitude = df$long, latitude = df$lat)

# First to EPSG:4326
df$coords_sf <- st_as_sf(df$coords,
  coords = c("longitude", "latitude"),
  crs = 4326
)

df$coords_sf
#> Simple feature collection with 2 features and 0 fields
#> Geometry type: POINT
#> Dimension:     XY
#> Bounding box:  xmin: 9.212 ymin: 48.109 xmax: 13.338 ymax: 52.221
#> Geodetic CRS:  WGS 84
#>                geometry
#> 1  POINT (9.212 48.109)
#> 2 POINT (13.338 52.221)

# Proof
deu_gg +
  geom_sf(data = df$coords_sf) +
  coord_sf(crs = 4326)

# Mercator
df$mercator <- st_transform(df$coords_sf, 3857)
df$mercator
#> Simple feature collection with 2 features and 0 fields
#> Geometry type: POINT
#> Dimension:     XY
#> Bounding box:  xmin: 1025475 ymin: 6125008 xmax: 1484779 ymax: 6840184
#> Projected CRS: WGS 84 / Pseudo-Mercator
#>                  geometry
#> 1 POINT (1025475 6125008)
#> 2 POINT (1484779 6840184)

# Proof
deu_gg +
  geom_sf(data = df$mercator) +
  coord_sf(crs = 3857)

# LAEA
df$eters_laea <- st_transform(df$coords_sf, 3035)
df$eters_laea
#> Simple feature collection with 2 features and 0 fields
#> Geometry type: POINT
#> Dimension:     XY
#> Bounding box:  xmin: 4262291 ymin: 2777584 xmax: 4549026 ymax: 3239817
#> Projected CRS: ETRS89-extended / LAEA Europe
#>                  geometry
#> 1 POINT (4262291 2777584)
#> 2 POINT (4549026 3239817)

# Proof
deu_gg +
  geom_sf(data = df$eters_laea) +
  coord_sf(crs = 3035)

^{Created on 2023-08-02 with reprex v2.0.2}