How Do I Convert from Map Coordinates from SWEREF99 TM to WGS84 Using R?

I am trying to use the sf package for R to convert the coordinates North 6855790, East 557541 from coordinate reference sysstem SWEREF99 TM to WGS84 and expecting the output North 61.830576°, East 16.092562° (in decimal, or North 61° 49′ 50,07″, East 16° 5′ 33,22″ in degrees, minutes, and seconds, although I would prefer decimal), as I get when using this website. I am, however, getting North 64.54803°, East 3.262454°.

The EPSG codes for SWEREF99 TM and WGS84 should be 3006 and 4326 respectively, as far as I understand.

For context: I plan to apply this conversion to a large dataset of coordinates using dplyr, once it works, so using the conversion website mentioned above is not practical. Alternative solutions, not using the sf package, would also work.

What am I doing wrong?

I have tried:

# This part’s here for reproducibility.
# Ignore it if you already have sf installed
chooseCRANmirror(ind=1)
install.packages('sf')

# The actual code.
library(sf)
data.frame(x = 6855790, y = 557541) |>
  st_as_sf(coords = c('x', 'y')) |>
  st_set_crs(3006) |>
  st_transform(4326)

Which gives me the following output:

## Simple feature collection with 1 feature and 0 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 64.54803 ymin: 3.262454 xmax: 64.54803 ymax: 3.262454
## Geodetic CRS:  WGS 84
##                    geometry
## 1 POINT (64.54803 3.262454)

I searched the forums before posting. For example this post is similar but uses the package rdal, which is no longer on CRAN as of October 16, 2023.

Switch values

library(sf)
#> Linking to GEOS 3.11.0, GDAL 3.5.3, PROJ 9.1.0; sf_use_s2() is TRUE
data.frame(x = 557541, y = 6855790) |>
  st_as_sf(coords = c('x', 'y')) |>
  st_set_crs(3006) |>
  st_transform(4326)
#> Simple feature collection with 1 feature and 0 fields
#> Geometry type: POINT
#> Dimension:     XY
#> Bounding box:  xmin: 16.09256 ymin: 61.83058 xmax: 16.09256 ymax: 61.83058
#> Geodetic CRS:  WGS 84
#>                    geometry
#> 1 POINT (16.09256 61.83058)

That is c("longitude", "latitude"). Looks promising

library(rnaturalearth)
library(ggplot2)
library(sf)

data.frame(x = 557541, y = 6855790) |>
  st_as_sf(coords = c('x', 'y')) |>
  st_set_crs(3006) |>
  st_transform(4326) -> pt

swe <- ne_countries(country = "sweden", scale = "large", returnclass = "sf")

ggplot(data = swe) +
  geom_sf() +
  geom_sf(data = pt)

^{Created on 2023-10-20 with reprex v2.0.2}