Clipping using R and ArcMap do not match

Question

I am carrying out a Raster clipping from a shapefile, using its exact shape.

I have done a clip using the Clip (Data Management) in ArcMap. The "Using Input Features for Clipping Geometry (optional) checkbox is checked whereas the "Maintain Clipping Extent (optional)" is not. Then I perform a clip using mask() function in R. Values of NDVIs are quite different. Please see the image attached.

I think this might be a problem of projection because I have a warning in R that it doesn't recognise EPS 3301. Shall I take the raster and shapefile to WGS84?

I post the R code:

library(tools)
library(raster)
library(rgdal)
library(rgeos)
library(sp)
library(RColorBrewer)

Study_Areas= list("a", "b")##TEST
##System Location Paths
#Inputs
VIS = "route to where i store the VIs"
MultiSpectral = "route where I store the multispectral bands"
Mask = "Route where I store the polygons for clipping"
Shapefile <- list.files(file.path(Mask), pattern = "\\.shp$",full.names = TRUE)

for(study_A in Study_Areas){
  study_A_cap <- toTitleCase(study_A)
  
  
  ##Checking the study areas
  
  cat("---initialising If statement---\n")
  
  print(study_A_cap)
  
  allBands <- list.files(file.path(MultiSpectral,study_A),pattern = "tif$",full.names = TRUE)
  
  print(paste("The bands in the study area:", study_A_cap),sep=" ")
  print(allBands)
  
  ##Calculate NDVI
  
  #1. stack bands because you need to index them
  stackBands <- stack(allBands)
  brikBands_br <- brick(stackBands)
  #find the band indexes
  
  print(brikBands_br[[4]])##RED EDGRE
  print("\n")
  print(brikBands_br[[3]])## RED
  print("\n")
  print(brikBands_br[[2]])##NIR
  print("\n")
  print(brikBands_br[[1]])##GREEN
  ## Warning: Discarded datum Unknown based on GRS80 ellipsoid in CRS definition #-> anyway, I can add the bands on ArcMap. (EPSG: 3301)
  
  # I am creating an empty list here which is going to have all the vegetation indices 
  VI_list = list()
  
  
  #Create the vegetation indices and visualise them in separated windows
  
  ndvi <- (brikBands_br[[2]] - brikBands_br[[3]]) / (brikBands_br[[2]] + brikBands_br[[3]])
  ndwi <-(brikBands_br[[1]] - brikBands_br[[2]]) / (brikBands_br[[1]] + brikBands_br[[2]])
  
  windows()
  plot(ndvi, main= paste("NDVI",study_A_cap, sep=" of "))
  windows()
  plot(ndwi, main= paste("NDWI",study_A_cap, sep=" of "))
  
  
  
  
  ## append the vegetation indices to the empty list
  VI_list = append(VI_list,list(ndvi,ndwi))
  
  VI_list
  
  ## fill (MANUALLY) this list, in ORDER, 
  names_of_VI = list("ndvi","ndwi")
  
  names_of_VI
  cat("\nVEGETATION INDEX HAVE BEEN CALCULATED")
  cat("\n")
  
  
  
  ##from this loop, I want to get the names of the vegetation indices in order to store them
  #so, for each of them, I carry out the following analysis and their output (saved VI) will have the name
  for(names in names_of_VI){
    
    ## and now, process a mask extraction for each study area
    for(shp in Shapefile){
      
      theShp <- readOGR(shp) ##this allows us to read the shapefile format, using OGR
      
      cat("Printing the study area: for clipping \n")
      print(shp)
      cat("\n")
      
      Split_path<- strsplit(shp, split = "/") ##get the name of the study area ONLY
      
      StAreaName <- strsplit(Split_path[[1]][10],split = ".shp")
      
      print(paste("\nStudy area equal ?? ---", StAreaName==study_A,"---"),sep= "")
      
      ##show the study areas
      plot(theShp, main = paste("Study area:",StAreaName, sep= " "),
           axes = TRUE,
           border = "blue",
           na.rm=TRUE)
      ##performs mask extraction for the study area!
      if(StAreaName==study_A){
        
        print("se ejecuta StAreaName == study_A, BIEN! ")
        
        for(ras in VI_list){ ##this is each raster of VI index in the list created before
          
          cat("The VI raster:\n")  
          print(ras)
          cat("\n")  
          
          ## PERFORM THE CLIP  (EXTRACT BY MASK) !!!
          cat("Doing clip\n")
          raster_clip <- crop(ras, theShp)
          
          plot(raster_clip, main = "The VI cropped")
          
          cat("Clip is done\n")
          
          ##we are going to save the Vegetation index in this route:
          print(file.path(Image_Masks_R ,study_A_cap,paste(study_A_cap,ras),sep = "_"))
          cat("\n")
          
          print("WRITING RASTER !!!")
          writeRaster(raster_clip, filename=file.path(Image_Masks_R ,study_A_cap,paste(study_A_cap,names,".tif",sep = "_")),
                      options=c('TFW=YES','CPG = YES'),
                      ##format = "GTiff", ## as geotiff
                      ##dataType = 'INT25', ##USE ?datatype to see options,
                      overwrite = TRUE)
        }
        
        
      }else{
        print("Nothing to output")
      }

      
    }##en for(shp in shapefile)
    
  }## end for(names in names_of_VI)
  cat(paste("\nFINISHEDTHIS AREA:\n",study_A_cap ))

}##cat("\nFINISHED\n")

Answer 1

I have used the following two functions in order to perform the clip:

mask(crop())

And then, all the NDVI rasters are exactly the same.