How to correct erroneous CRS in the function coord_sf() in ggplot after merging shapefiles together to zoom into a specific region in R

Question

Issue

I have three shapefiles that exhibit our study area and contain the "EPSG - 4326". After researching, the exact EPSG for our study area is "3355 with transformation: 15846" for the Suez Canal and Red Belt regions. I need to zoom in just to our study area, which is the Northern Reefs in The Red Sea in Egypt as we study a marine species within this region.

I merged the shapefiles together using the package ggplot2, but the outcome was a map of the world (see below).

I played around with both ESPG numbers "3355" and "4326" in the function sf_coord() in the sf package to try and plot the shapefiles zoomed in, but it doesn't make a difference with the output.

I thought the function geom_sf() automatically took the CRS from the first layer and plotted the shapefiles based on the information contained in the shapefile.

I can see that the shapefiles have been plotted, but I'm not sure how to zoom in with R to just our study area. Please see the screenshot below which shows a zoomed out and zoomed in screenshot of the three shapefiles plotted in QGIS.

The typical coordinates for our study area are 27 degrees latitude and 33 degrees longitude.

Does anyone know how to do this? Do you need to specify exact latitude and longitude coordinates?

I'm a novice with GIS analysis in R.

Many thanks if anyone can help.

R-Code

#Load the shape file for our study area in The Red Sea
Map_Study_Area<-read_sf("/Users/GIS_Shapefiles/Map_RedSea_middleEast.shp")
Map_Study_Area

#Load the shape file of the reefs in our study area 
Red_N_Reefs<-read_sf("/UsersGIS_Shapefiles/Reef_Northen_Islands.shp")
Red_N_Reefs

#Load the shape for the Northern Red Sea Islands 
R_N_Islands<-read_sf("/Users/GIS_Shapefiles/Islands_Egypt.shp")
R_N_Islands

#Explore the CRS of the study area 
st_crs(Map_Study_Area)

Coordinate Reference System:
  User input: WGS 84 
  wkt:
GEOGCRS["WGS 84",
    DATUM["World Geodetic System 1984",
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["latitude",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["longitude",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    ID["EPSG",4326]]

**#Egypt Gulf of Suez S-650 TL / Red Belt
#EPSG:3355 with transformation: 15846

#Plot the shape files together using ggplot
Study_Area_Shp<-ggplot() + 
                      geom_sf(data = Map_Study_Area, color = "gray", fill = "grey") +
                      geom_sf(data = R_N_Islands, color = "green", fill = "green") +
                      geom_sf(data = Red_N_Reefs, color = "red", fill = "red") +
                      ggtitle("The Northern Reefs of The Red Sea") + 
                      coord_sf(crs = st_crs(3355))

Output

Zoomed out in QGIS

Zoomed in: Desired Output