Creating a single map composed of three separate shapefiles

I have the following shapefiles:

qmasp <- st_read("Rock_Lobster_(Spiny_Red)_QMAs.shp", quiet = TRUE)
statShp <- st_read("Rock_Lobster_Statistical_Areas.shp", quiet = TRUE)
NE_countries <- st_read("ne_110m_admin_0_countries", quiet = TRUE)

I want them to appear together as a single map. Each shapefile looks like this (see images):

and I am using the following code:

transform QMA shapefile to an sf object

qmasf=st_as_sf(qmasp)
qmasf4326=st_transform(qmasf, 4326)

then NZ UTM for plotting

qmasf2134=st_transform(qmasf, 2193)#2134
qmasf2134 %>% 
  ggplot()+
  geom_sf(aes(fill=FishstockC))

transform Stat Area Shapefile to sf object

statsf=st_as_sf(statShp)
statsf4326=st_transform(statsf, 4326)

then NZ UTM for plotting

statsf2134=st_transform(statsf, 2134)#check if NZ2000 is better
stats<-statsf2134 %>%
 ggplot()+
  geom_sf(fill= 'White')+
  geom_sf_text(
    aes(label = Statistica),
    size = 4 / .pt
  )
stats

Plot Map

theme_set(theme_bw())

newzealand <- ne_countries(country="new zealand", type="countries", scale = 'large', returnclass = "sf")

ggplot(data = newzealand) +
  geom_sf() +
  coord_sf(crs = 2193)

I can then plot data points using coordinate data onto specific polygons from the QMA shapefile on the map with this code:

to prepare data for plotting on map

samples_sf <- sf::st_as_sf(d2, coords=c("Longitude","Latitude"), crs=4326)

Subset CRA3 polygon from the QMA shapefile

CRA3_shapefile<-qmasf2134 %>%
  filter(FishstockC %in% c("CRA2","CRA3"))
CRA3_shapefile

which looks like this (see below):

or this depending on which code I use:

My questions are, (i) How do I produce a map that has all three shape files on the one image, essentially layered on top of one another? (ii) Why is my map of New Zealand so small and how do I make it the same size as the other shapefiles, without losing Chatham Island (the third island of New Zealand)? This island is to the East of the north and south island of New Zealand.

Any help will be greatly appreciated!

Tēnā koe e hoa, here's a repex given we don't have access to your data. You can think of coord_sf() as a zoom function. To get the extent of your data, use st_bbox() and work forwards from there. With some trial and error, you'll achieve the desired result. The example image uses the xlim/ylim (long/lat) as stated. For each layer, add a new geom_sf() as shown below. You have to explicitly state the data = each time, but I prefer this approach as I find it easier to follow. The geom_sf() are plotted in order, so your background layers go first:

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

newzealand <- ne_countries(country = "new zealand", 
                           type = "countries", 
                           scale = 'large', 
                           returnclass = "sf") %>%
  st_transform(2193)

# Dummy CRA3_shapefile
CRA3_shapefile <- data.frame(cra = rep(c("CRA2", "CRA3"), each = 4),
                             lat = c(5700000, 5800000, 5800000, 5700000,
                                     5800000, 5900000, 5900000, 5800000),
                             long = c(2000000, 2000000, 2200000, 2200000)) %>%
  st_as_sf(coords = c("long","lat")) %>%
  group_by(cra) %>%
  summarise(geometry = st_combine(geometry)) %>%
  st_cast("POLYGON") %>%
  st_set_crs(2193)

# Dummy point data. set.seed() for replicable result, only needed for dummy point data
set.seed(123)
sf_points <- data.frame(lat = sample(5700000:5900000, 20),
                        long = sample(2000000:2200000, 20)) %>%
  st_as_sf(coords = c("long","lat")) %>%
  st_cast("POINT") %>%
  st_set_crs(2193)

# Get extent of newzealand data as basis for setting coord_sf limits
st_bbox(newzealand)
xmin    ymin    xmax    ymax 
1092701 4165390 3357853 9024853 

ggplot() +
  geom_sf(data = CRA3_shapefile,
          aes(fill = cra)) +
  geom_sf(data = newzealand,
          colour = "black") +
  geom_sf(data = sf_points,
          colour = "grey40",
          size = 1.5) +
  coord_sf(xlim = c(1002701, 2557853),
           ylim = c(4165390, 6204853))

Update based on OP's comment

Often it is more desirable to filter out the unnecessary geometries rather than using coord_sf(). It's trickier in this case as your NZ data are MULTIPOLYGON so every individual polygon will share the same non-geometry values (attributes in ArcGIS parlance). In other words, there's no way of knowing which polygons represent Avaiki Nui/Cook Islands other than by their long/lat values.

In your case, the simplest approach is to 'borrow' the max ylim value from the previous step e.g 6204853 and to filter out everything north of 6204853. This approach involves breaking your MULTIPOLYGON into individual polygons using st_cast():

# Convert MULTIPOLYGON TO POLYGON (you can safely ignore the resultant warning)
nz1 <- newzealand %>%
  st_cast("POLYGON") %>%
  mutate(id = 1:n()) # This column added as a key for filtering

# 1. Extract coordinates from nz1
# 2. Filter out polygon ids north of max ylim value
# 3. left_join() polygon ids to nz1
# 4. Filter out unwanted polygons
# NOTE: the Y and L2 columns referenced here are 
#       default names generated by data.frame(st_coordinates(nz1))
nz2 <- data.frame(st_coordinates(nz1)) %>%
  filter(Y <= 6204853) %>%
  select(L2) %>%
  distinct() %>%
  left_join(nz1, ., by = join_by(id == L2), keep = TRUE) %>%
  filter(!is.na(L2))