R - ggplot mapping data on international country

I am trying to map a zoomed-in view along the Mwanihana forest in Tanzania.

I have looked into Plotting latitude and longitude points in r, but can't get it to work for my dataset.

Data <-

date     forestsite        primaryinterface secondaryinterface                   
<chr>    <chr>             <chr>            <chr>
8.27.13  UMNP-HQ	       ecotourism	    forest use/PD
8.28.13  UMNP-campsite3	   ecotourism	    forest use
8.28.13  UMNP-campsite3	   ecotourism	    forest use
8.28.13  UMNP-hondohondoa  ecotourism	    forest use/PD
8.28.13	 UMNP-hondohondob  ecotourism       forest use/PD/crop-raiding
8.29.13	 UMNP-njokamoni	   ecotourism	    forest use
8.29.13	 UMNP-mangabey	   ecotourism	    forest use   

anthrochange                     species       geometry
<chr>                            <chr>         <S3: sfc_POINT>      
small trails, snares	         BWC           <S3: sfc_POINT>
roads, buildings	             RC            <S3: sfc_POINT>
roads, campsite	                 RC            <S3: sfc_POINT>
roads, campsite	                 BWC           <S3: sfc_POINT>
none	                         SY            <S3: sfc_POINT>
clearcut road, houses, crops     YB            <S3: sfc_POINT>
trail to njokamoni	             RC            <S3: sfc_POINT>
small trails, habituated troop	 SM            <S3: sfc_POINT>

Map data and ggplot attempt

library(elevatr)
library(rgdal)

se.elevations <- get_elev_raster(locations = world, z = 7, clip = "locations") 
se.elevation.df <- raster::as.data.frame(se.elevations, xy = TRUE) %>%
  rename(elevation =  3) %>%
  na.omit() %>%
  filter(elevation >= 0)
xlimit <- c(36.85, 37)
ylimit <- c(-7.76, -7.88)

ggplot()+
  geom_tile(data = se.elevation.df, aes(x = x, y = y, fill = elevation))+
  scale_fill_gradientn(colours= c("gray", "black"))+
  geom_sf(data = world, fill = NA)+
  xlim(xlimit)+
  ylim(ylimit)+
  geom_sf(data = data, aes(color = forestsite, shape = primaryinterface), size = 3) +
    coord_sf(xlim = c(36.85, 37), ylim = c(-7.76, -7.88), expand = FALSE)

I am able to get this map -

from this code

ggplot() +
  geom_sf(data = data, aes(color = species, shape = primaryinterface), size = 3) +
  coord_sf(xlim = c(36.85, 37), ylim = c(-7.76, -7.88), expand = FALSE) -> ggplottt

but anytime I try to add the mymap data to include the Tanzanian lat, long, I get various error messages.

Any advice is greatly appreciated.

Deput

structure(list(date = c("8.27.13", "8.28.13", "8.28.13", "8.28.13", 
"8.28.13", "8.29.13", "8.29.13"), forestsite = c("UMNP-HQ", "UMNP-campsite3", 
"UMNP-campsite3", "UMNP-hondohondoa", "UMNP-hondohondob", "UMNP-njokamoni", 
"UMNP-mangabey"), primaryinterface = c("ecotourism", "ecotourism", 
"ecotourism", "ecotourism", "ecotourism", "ecotourism", "ecotourism"
), secondaryinterface = c("forest use/PD", "forest use", "forest use", 
"forest use/PD", "forest use/PD/crop-raiding", "forest use", 
"forest use"), anthrochange = c("roads, buildings", "roads, campsite", 
"roads, campsite", "none", "clearcut road, houses, crops", "trail to njokamoni", 
"small trails, habituated troop"), species = c("RC", "RC", "BWC", 
"SY", "YB", "RC", "SM"), geometry = structure(list(structure(c(36.883988, 
-7.844929), class = c("XY", "POINT", "sfg")), structure(c(36.884545, 
-7.849439), class = c("XY", "POINT", "sfg")), structure(c(36.884545, 
-7.849439), class = c("XY", "POINT", "sfg")), structure(c(36.887065, 
-7.833501), class = c("XY", "POINT", "sfg")), structure(c(36.891369, 
-7.832091), class = c("XY", "POINT", "sfg")), structure(c(36.878388, 
-7.82738), class = c("XY", "POINT", "sfg")), structure(c(36.87753, 
-7.8232), class = c("XY", "POINT", "sfg"))), class = c("sfc_POINT", 
"sfc"), precision = 0, bbox = structure(c(xmin = 36.87753, ymin = -7.849439, 
xmax = 36.891369, ymax = -7.8232), class = "bbox"), crs = structure(list(
    input = "EPSG:4326", wkt = "GEOGCRS[\"WGS 84\",\n    ENSEMBLE[\"World Geodetic System 1984 ensemble\",\n        MEMBER[\"World Geodetic System 1984 (Transit)\"],\n        MEMBER[\"World Geodetic System 1984 (G730)\"],\n        MEMBER[\"World Geodetic System 1984 (G873)\"],\n        MEMBER[\"World Geodetic System 1984 (G1150)\"],\n        MEMBER[\"World Geodetic System 1984 (G1674)\"],\n        MEMBER[\"World Geodetic System 1984 (G1762)\"],\n        MEMBER[\"World Geodetic System 1984 (G2139)\"],\n        ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n            LENGTHUNIT[\"metre\",1]],\n        ENSEMBLEACCURACY[2.0]],\n    PRIMEM[\"Greenwich\",0,\n        ANGLEUNIT[\"degree\",0.0174532925199433]],\n    CS[ellipsoidal,2],\n        AXIS[\"geodetic latitude (Lat)\",north,\n            ORDER[1],\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        AXIS[\"geodetic longitude (Lon)\",east,\n            ORDER[2],\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n    USAGE[\n        SCOPE[\"Horizontal component of 3D system.\"],\n        AREA[\"World.\"],\n        BBOX[-90,-180,90,180]],\n    ID[\"EPSG\",4326]]"), class = "crs"), n_empty = 0L)), row.names = c(NA, 
-7L), sf_column = "geometry", agr = structure(c(date = NA_integer_, 
forestsite = NA_integer_, primaryinterface = NA_integer_, secondaryinterface = NA_integer_, 
anthrochange = NA_integer_, species = NA_integer_), .Label = c("constant", 
"aggregate", "identity"), class = "factor"), class = c("sf", 
"tbl_df", "tbl", "data.frame"))

I think that is better an interactive map. The data is df. You need get the coordinate in a good format because your object is a sfc_POINT.

I'm work with you dput()

library(sf)
library(leaflet)
library(raster)
library(sp)

# Extract coordinates
df$latitudes <- st_coordinates(df$geometry)[, 2]
df$longitudes <- st_coordinates(df$geometry)[, 1]

adm <- getData('GADM', country='TZA', level=1)

popup <- paste0("<strong>POPULATION: </strong>", 
                adm$NAME_1)

# Tanzania interactive map
TZA_MAP <- leaflet() %>% 
  addTiles() %>%
  addPolygons(data=adm, weight = 2, col = 'blue', fillColor = "#FEF800", popup=popup,
              fillOpacity = 0.1) |> 
  addCircles(data = df ,lat= df$latitudes , lng= df$longitudes,
             label = c(df$primaryinterface, df$secondaryinterface), labelOptions = 
  labelOptions(noHide = T, direction = "bottom"),
             group = c(df$primaryinterface, df$secondaryinterface)) |> 
  addLayersControl(overlayGroup= c(df$primaryinterface, df$secondaryinterface),
                       options = layersControlOptions(collapsed=F))

TZA_MAP 

You could select for filter primaryinterface and secondaryinterface, but you can change in function what you want to do.