95m Difference between R .laz file spTransform() and the Result from epsg.io - why?

In comparing R’s result using spTransform() to reproject an .laz file, I find that its result in one case is 95m off that given by epsg.io for the same reprojection.

I have found that results from R and LAStools (run via QGIS) with each other and that results from epsg.io and from https://mygeodata.cloud/cs2cs/ agree with each other, so something is systematically amiss with what is happening on my PC vs what is happening with the online transformations.

• Can you help me understand why this is happening? Do I misunderstand transformation as being the same thing as reprojection?

Full example below, if helpful:

> # Test to compare R's projection transform with that of epsg.io
> 
> #rm(list = ls()); gc()
> 
> library(lidR) # also loads raster and sp packages
> library(sf)
> 
> #setwd("C:/MappingProjects/Project_20210209_CeriseParkMontrose/Lidar/")
> # download .laz file from TNM USGS Download, using "-107.882239°, 38.460569°"
> # setwd() here
> las.in1 <- readLAS("USGS_LPC_CO_WestCentral_2019_A19_w1024n1778.laz")
> wkt(las.in1)
[1] "COMPD_CS["NAD83(2011) / Conus Albers + NAVD88 height - Geoid12B (m)",PROJCS["NAD83(2011) / Conus Albers",GEOGCS["NAD83(2011)",DATUM["NAD83 (National Spatial Reference System 2011)",SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],AUTHORITY["EPSG","1116"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","6318"]],PROJECTION["Albers_Conic_Equal_Area"],PARAMETER["standard_parallel_1",29.5],PARAMETER["standard_parallel_2",45.5],PARAMETER["latitude_of_center",23],PARAMETER["longitude_of_center",-96],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["meter",1,AUTHORITY["EPSG","9001"]],AXIS["X",EAST],AXIS["Y",NORTH],AUTHORITY["EPSG","6350"]],VERT_CS["NAVD88 height - Geoid12B (m)",VERT_DATUM["North American Vertical Datum 1988",2005,AUTHORITY["EPSG","5103"]],UNIT["meter",1,AUTHORITY["EPSG","9001"]],AXIS["Up",UP],AUTHORITY["EPSG","5703"]]]"
> crs(las.in)
CRS arguments:
 +proj=aea +lat_0=23 +lon_0=-96 +lat_1=29.5 +lat_2=45.5 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +vunits=m +no_defs 
> 
> # attempt to transform using OGSWKT
> # las.out <- spTransform(las.in1, CRS('PROJCS["WGS 84 / UTM zone 13N",GEOGCS["WGS 84",
> # DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],
> # AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],
> # UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],
> # PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],
> # PARAMETER["central_meridian",-105],PARAMETER["scale_factor",0.9996],
> # PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,
> # AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],
> # AUTHORITY["EPSG","32613"]]'))
> 
> # attempt to transform using EPSG code (same result as with OGSWKT above):
> las.out <- spTransform(las.in1, CRS("+init=epsg:32613"))
> 
> crs(las.out)
CRS arguments:
 +proj=utm +zone=13 +datum=WGS84 +units=m +no_defs 
> las.in1@bbox
       min      max
x -1024000 -1023000
y  1778000  1779000
> las.out@bbox
        min       max
x  248879.7  249980.1
y 4260326.2 4261404.4
> 
> # per epsg.io (5070, 5071, 5072, 6350)
> # https://epsg.io/transform#s_srs=6350&t_srs=32613&x=-1024000.0000000&y=1778000.0000000
> minX = 248974.68   # 5069 = 248941.74
> minY = 4260326.01   # 5069 = 4260385.05
> # same output as with my GEoCloud: https://mygeodata.cloud/cs2cs/
> 
> # Error:
> X.err <- las.out@bbox[1, 1] - minX
> X.err
[1] -95.02
> # Y Error
> Y.err <- las.out@bbox[2, 1] - minY
> Y.err
[1] 0.14
> sqrt((X.err)^2 + (Y.err)^2)
[1] 95.0201
> 
> # compare to reprojection in LAStools - agrees with R, not with epsg.io
> las.lastools <- readLAS("n1024_1778_reproject.laz")
> las.lastools@bbox
        min       max
x  248879.7  249980.1
y 4260326.2 4261404.4