rcoordinatesgeospatialr-spadehabitathr

Why is my MCP giving me an area size that is too small?


This is my first time doing any spatial statistics and I am slowly learning. I have seen a few other people here with similar problems, however none of the answers seemed to work for me.

I am trying to do a MCP, but when I look at the results the area size is far smaller than it should be. However when I plot the polygons over a satellite map it looks correct.

My data is set up like this:

set.seed(1)
animal_sp <- data.frame(ID = rep(c("S1", "S2"), each = 10),
                         Long = runif(20, 34.2002, 34.55328),
                         Lat = runif(20, -19.11998, -18.82139))

And the code that I have used is as follows:

library(sp)

sp::coordinates(animal_sp) <- c("Long", "Lat") # define the variables

proj4string(animal_sp) <- CRS("+init=epsg:4326") # Set the projection




# Create the MCP using the adehabitatHR package:                                                

library(adehabitatHR)

animal_mcp <- mcp(animal_sp, 
                 percent = 100,
                 unout = "m2") # Make the MCP

animal_mcp # look at the results

The results that I get are:

Object of class "SpatialPolygonsDataFrame" (package sp):

Number of SpatialPolygons:  2

Variables measured:
   id       area
S1 S1 0.02919310
S2 S2 0.03671279

I should be getting an area that is ~260 km².

Any help would be greatly appreciated!


Solution

  • As @r2evans alluded to, you need to project your sp object to its corresponding UTM Zone CRS before running mcp(). Currently, mcp() is treating your coordinates as literal metres, hence the small area values.

    Here's how to correct your issue. Note that you are 'lucky' in that all of your data fall inside the same UTM Zone. If your data span multiple UTM Zones, you will need to search online for how to manage this issue/choose the correct CRS relative to the extent of your data.

    library(sp)
    library(adehabitatHR)
    
    set.seed(1)
    animal_sp <- data.frame(ID = rep(c("S1", "S2"), each = 10),
                            Long = runif(20, 34.2002, 34.55328),
                            Lat = runif(20, -19.11998, -18.82139))
    
    # Function to return UTM Zone/s
    utm_epsg <- \(lon, lat) {
    
      utm_zone <- floor((lon + 180) / 6) + 1
      if(lat >= 0) {
        epsg <- 32600 + utm_zone
      } else {
        epsg <- 32700 + utm_zone
      }
      
    }
    
    # Get UTM Zone/s for animal_sp
    unique(mapply(utm_epsg, animal_sp$Long, animal_sp$Lat))
    # [1] 32736
    
    # Define coordinates
    coordinates(animal_sp) <- c("Long", "Lat")
    
    # Set intial projection
    proj4string(animal_sp) <- CRS("EPSG:4326")
    
    # Project to relevant UTM Zone CRS
    animal_sp_utm <- spTransform(animal_sp, "EPSG:32736")
    
    # Create the MCP
    animal_mcp <- mcp(animal_sp_utm, 
                      percent = 100,
                      unout = "m2") # Make the MCP
    
    animal_mcp
    # Object of class "SpatialPolygonsDataFrame" (package sp):
    # 
    # Number of SpatialPolygons:  2
    # 
    # Variables measured:
    #    id      area
    # S1 S1 601430654
    # S2 S2 461876993