I’m currently doing a project trying to determine the home ranges of sharks using minimum convex polygons as a measurement. So far I’ve managed to successfully make the minimum convex polygons (MCPs) which tell you how large the home range for each animal is, however I’m running into problems when trying to place them on a map, where the map is created but the MCPs are plotted in the completely incorrect area of the world due to issues with coordinate transformations.

I think I’ve narrowed down the problem to the part of my code when I convert the data from a UTM 40 zone to a WGS84 CRS (coordinate system) to make it compatible with Google Maps which uses the WGS84 CRS. The coordinates my MCPs should be in are approximately 71.96, -5.24 which are located in the Chagos Archipelago. However, these coordinates change to 52.51, -4.72e-05 when I convert the data to the new CRS for use with the Google Maps. These new coordinates are very small and not in the Chagos area at all.

The first couple lines of code are converting the MCP data (named SP_cat2 and SP_cat2MCP) into EPSG:32640 (named SF32map_cat2 / cat2MCP), which is the UTM zone 40 / Chagos. When I check the coordinates of these outputs, they are correct:

SF32map_cat2 <- st_as_sf(SP_cat2, crs = 32640) 
SF32map_cat2MCP<- st_as_sf(SP_cat2MCP, crs = 32640)

However, when I convert them to EPSG:4326 (to be used with Google Maps) and check the coordinates, they have changed to the very small coordinates. This is the code used:

SFmap_cat2 <- st_transform(SF32map_cat2, crs = 4326) 
SFmap_cat2MCP <- st_transform(SF32map_cat2MCP, crs = 4326)

I’ve also used st_as_sf() for the second lot of code, and this does successfully put the MCPs on a map, but again, with the wrong coordinates.

After converting to the new CRS, the coordinates are completely wrong and this messes up the map as it either plots the MCPs very close to the equator which isn’t where they should be at all, or plots just the Chagos Archipelago with none of my data on top. This is the mapping code in case this is problematic as well:

ggmap(chagos_map) + 
  geom_sf(data = SFmap_cat2, aes(color = as.factor(code)), size = 2, inherit.aes = FALSE) +
  geom_sf(data = SFmap_cat2MCP, fill = NA, color = "black", alpha = 0.5, inherit.aes = FALSE) +
  ggtitle("MCP Silvertip Cat 2 Chagos") +
  theme_minimal()

I’ve plotted the MCPs themselves without the map behind them and they do plot correctly, but they are located in the complete wrong coordinates very close to the equator. I’ve double checked all the classes and structures of the data and they are all correct. I can’t see where the problem is at all.

Here are the outputs of SP_cat2:

new("SpatialPointsDataFrame", data = structure(list(code =     c("12964", 
    "12964", "12964", "12964", "12964", "12964")), row.names =     c(NA, 
    6L), class = "data.frame"), coords.nrs = 2:3, coords =     structure(c(71.9684, 
    71.9684, 71.9684, 71.9684, 71.9684, 71.9698, -5.2417,     -5.2417, 
    -5.2417, -5.2417, -5.2417, -5.3302), dim = c(6L, 2L),     dimnames = list(
    NULL, c("receiver_lon", "receiver_lat"))), bbox =     structure(c(71.9684, 
    -5.3302, 71.9698, -5.2417), dim = c(2L, 2L), dimnames = list(
    c("receiver_lon", "receiver_lat"), c("min", "max"))),     proj4string = new("CRS", 
    projargs = "+proj=utm +zone=40 +datum=WGS84 +units=m         +no_defs"))

and SP32map_cat2MCP:

structure(list(code = c("12964", "12964", "12964", "12964", "12964", 
"12964"), geometry = structure(list(structure(c(71.9684, -5.2417
), class = c("XY", "POINT", "sfg")), structure(c(71.9684, -5.2417
), class = c("XY", "POINT", "sfg")), structure(c(71.9684, -5.2417
), class = c("XY", "POINT", "sfg")), structure(c(71.9684, -5.2417
), class = c("XY", "POINT", "sfg")), structure(c(71.9684, -5.2417
), class = c("XY", "POINT", "sfg")), structure(c(71.9698, -5.3302
), class = c("XY", "POINT", "sfg"))), class = c("sfc_POINT", 
"sfc"), precision = 0, bbox = structure(c(xmin = 71.9684, ymin = -5.3302, 
xmax = 71.9698, ymax = -5.2417), class = "bbox"), crs = structure(list(
    input = "+proj=utm +zone=40 +datum=WGS84 +units=m +no_defs", 
    wkt = "PROJCRS["unknown",n    BASEGEOGCRS["unknown",n        DATUM["World Geodetic System 1984",n            ELLIPSOID["WGS 84",6378137,298.257223563,n                LENGTHUNIT["metre",1]],n            ID["EPSG",6326]],n        PRIMEM["Greenwich",0,n            ANGLEUNIT["degree",0.0174532925199433],n            ID["EPSG",8901]]],n    CONVERSION["UTM zone 40N",n        METHOD["Transverse Mercator",n            ID["EPSG",9807]],n        PARAMETER["Latitude of natural origin",0,n            ANGLEUNIT["degree",0.0174532925199433],n            ID["EPSG",8801]],n        PARAMETER["Longitude of natural origin",57,n            ANGLEUNIT["degree",0.0174532925199433],n            ID["EPSG",8802]],n        PARAMETER["Scale factor at natural origin",0.9996,n            SCALEUNIT["unity",1],n            ID["EPSG",8805]],n        PARAMETER["False easting",500000,n            LENGTHUNIT["metre",1],n            ID["EPSG",8806]],n        PARAMETER["False northing",0,n            LENGTHUNIT["metre",1],n            ID["EPSG",8807]],n        ID["EPSG",16040]],n    CS[Cartesian,2],n        AXIS["(E)",east,n            ORDER[1],n            LENGTHUNIT["metre",1,n                ID["EPSG",9001]]],n        AXIS["(N)",north,n            ORDER[2],n            LENGTHUNIT["metre",1,n                ID["EPSG",9001]]]]"), class = "crs"), n_empty = 0L)), sf_column = "geometry", agr = structure(c(code = NA_integer_), levels = c("constant", 
"aggregate", "identity"), class = "factor"), row.names = c(NA, 
6L), class = c("sf", "data.frame"))