Mapping Open Map Raods to Terrain lat/long DEM - Problems with CRS transform Maybe

[rhgof]

Thanks for the Rayshader packege. Very clever and simple to use.
I am following the Bryce tutorial that adds Open Street Map data to terrain.
The Rayshader tutorial uses USM metres - but I am rendering in lat/long

I am new to the GIS domain - and still learning.
I am posting here as could not find a forum for support.

Issue
The problem is that when I overlay the transformed roads - they are overlayed incorrectly.
The roads are scale/down/shortened along the latitude (vertically). Longitudinal scale looks about right (I think).

I think this is due to the different CRS used by Open Street Map to CRS terrain DEM file.
But Iam not sure.

My guess is with the transform of OPenStreetMap Data below is not working.

roadPath = st_transform(roads$osm_lines, crs=crs(localtif))

Basically the output roadPath is not transformed at all.
See at bottom of code where I plot the untransformed roads$osm_lines with identical results.

The osm_lines comes back as OpenSteetMap lat/long which I think is a projected lat/long in EPSG:3875 that Open Street Map uses - not true lat/long.

See this good discussion https://gis.stackexchange.com/questions/48949/epsg-3857-or-4326-for-googlemaps-openstreetmap-and-leaflet?rq=1

I have tried multiple underlying CRS's in the terrain DEM file eg:
EPSG 4326,
EPSG:7844,
EPSG:4383
with same effect.

Reprex below
(NB - warning from curl that I cannot understand due to retrieving the remote file)
But the example works

Many thanks..

#!/usr/local/bin/Rscript
library(tidyverse)
library(ggplot2)

options(rgl.useNULL = TRUE)
library(rgl)
library(rayshader)
library(magick)
library(terra)
library(osmdata)
library(sf)

tifFile = "https://github.com/rhgof/sealevels/raw/main/PalmBeachNSW_Marine_5m_TopoBathy_DEM.tif"

# helpers
DEBUG=TRUE

epsgID="+init=epsg:4326"
#epsgID="+init=epsg:3857"

convert_coords = function(lat,long, from = CRS(epsgID), to) {
  data = data.frame(long=long, lat=lat)
  coordinates(data) <- ~ long+lat
  proj4string(data) = from
  #Convert to coordinate system specified by EPSG code
  xy = data.frame(sp::spTransform(data, to))
  colnames(xy) = c("x","y")
  return(unlist(xy))
}


# get the file

remotefile = tempfile() 
download.file(tifFile, remotefile)
localtif = raster::raster(remotefile)
elmat = raster_to_matrix(localtif,verbose = interactive())
unlink(remotefile)


# fix NA in the DEM for missing values in water
elmat[is.na(elmat)]<- -20.0

if (DEBUG) crs(localtif)

mapExt = extent(localtif)
if (DEBUG) mapExt

lat_range = c(ymin(mapExt),ymax(mapExt))
long_range = c(xmin(mapExt), xmax(mapExt))

utm_bbox = convert_coords(lat = lat_range, long=long_range, from = epsgID,  to = crs(localtif))
if (DEBUG) utm_bbox

extentzoom = extent(utm_bbox[1], utm_bbox[2], utm_bbox[3], utm_bbox[4])
if (DEBUG) extentzoom

# get roads from OpenStreet Map

osmap_bbox = c(xmin(mapExt),ymin(mapExt),xmax(mapExt),ymax(mapExt))
if (DEBUG) osmap_bbox

roads = opq(osmap_bbox) %>% 
  add_osm_feature("highway") %>% 
  osmdata_sf() 

roadPath = st_transform(roads$osm_lines, crs=crs(localtif))

# - check if we got the roads

ggplot(roadPath,aes(color=osm_id)) + 
  geom_sf() +
  theme(legend.position = "none") +
  labs(title = "Roads")

# - plot terrain

water_palette = colorRampPalette(c("darkblue", "dodgerblue", "lightblue"))(200)
bathy_hs = height_shade(elmat, texture = water_palette)

depth=0
Zscale=5  #5m grid

  elmat %>%
    sphere_shade(sunangle=45,zscale=Zscale, texture = "imhof3") %>%
    add_overlay(generate_altitude_overlay(bathy_hs, elmat, 0, depth))  %>%
    add_overlay(generate_line_overlay(roadPath,extent = extentzoom,linewidth = 2, heightmap = elmat)) %>% 
    add_shadow(ray_shade(elmat), 0.5) %>%
    plot_map()

#plot untransformed roads
  elmat %>%
    sphere_shade(sunangle=45,zscale=Zscale, texture = "imhof3") %>%
    add_overlay(generate_altitude_overlay(bathy_hs, elmat, 0, depth))  %>%
    add_overlay(generate_line_overlay(roads$osm_lines,extent = extentzoom,linewidth = 2, heightmap = elmat)) %>% 
    add_shadow(ray_shade(elmat), 0.5) %>%
    plot_map()