El recorte mediante R y ArcMap no coinciden

Question

Estoy realizando un recorte Raster de un shapefile, utilizando su forma exacta.

He hecho un clip utilizando el Clip (Gestión de Datos) en ArcMap. La casilla "Usar características de entrada para la geometría de recorte (opcional)" está marcada mientras que la casilla "Mantener la extensión del recorte (opcional)" no lo está. A continuación, realizo un recorte utilizando la función mask() en R. Los valores de los NDVI son bastante diferentes. Por favor, vea la imagen adjunta.

Creo que puede ser un problema de proyección porque tengo un aviso en R de que no reconoce el EPS 3301. ¿Debo llevar el raster y el shapefile a WGS84?

Pongo el código R:

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

He utilizado las dos funciones siguientes para realizar el clip:

mask(crop())

Y entonces, todos los rastreos NDVI son exactamente iguales.