setwd("C:/Users/Public/Documents/amsantac/data")
library(raster)

img <- raster("c_2004-2005_30_classif_03_v3.tif")
shp <- shapefile("footprints_2000_v2.shp")
img         # let's print some info about the imported raster

# Crop the image with the polygon shapefile
crop.img <- crop(img, shp)
crop.img    # compare the extent of the cropped image vs. the original image

plot(crop.img)  # if needed, install the package called 'Cairo'

hist(crop.img, maxpixels = 50137035)

# Export the cropped raster:
writeRaster(crop.img, "cropimg.tif", datatype = 'INT1U')

# Import arcpy and define the workspace
import arcpy
arcpy.env.workspace="C:/Users/Public/Documents/amsantac/data"
from arcpy.sa import *

# List the raster files in the workspace 
rasters = arcpy.ListRasters("*", "ALL")
for raster in rasters:
    print(raster)

# Set local variables
inRaster = "cropimg.tif"

# Execute MajorityFilter
outMajFilt = MajorityFilter(inRaster, "EIGHT", "HALF")

# Save the output 
outMajFilt.save("C:/Users/Public/Documents/amsantac/data/majfilter.tif")

# Verify the new raster file was created
rasters = arcpy.ListRasters("*", "ALL")
for raster in rasters:
    print(raster)

setwd("C:/Users/Public/Documents/amsantac/data")
library(raster)

crop.img <- raster("cropimg.tif")
maj.img <- raster("majfilter.tif")
plot(maj.img)    # Plot the filtered image

mpix <- 50137035
data.frame(crop.img = hist(crop.img, breaks = 0:8, maxpixels = mpix, plot = FALSE)$counts, 
      maj.img = hist(maj.img, breaks = 0:8, maxpixels = mpix, plot = FALSE)$counts)