1

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")

NDVI values

3
  • They are the same image, there are not any buildings or human made infrastructure. The problem is that for the same image, I get two different results
    – Ric_02
    Nov 5, 2020 at 7:18
  • Can you provide a smaller reprex tidyverse.org/help? I think a good step would be to check the CRS and projection of your rasters in R. You mention that R threw an error there, so that is where I would start.
    – GISHuman
    Nov 6, 2020 at 15:17
  • Originally, the rasters are in EPSG 3301 spatialreference.org/ref/epsg/3301; Perhaps this projection is not recognised by R
    – Ric_02
    Nov 6, 2020 at 15:28

1 Answer 1

0

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.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.