diff --git a/DESCRIPTION b/DESCRIPTION
index 8d64e20..7d0ba68 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,6 +1,6 @@
 Package: GLCMTextures
 Title: GLCM Textures of Raster Layers
-Version: 0.4
+Version: 0.4.1
 Authors@R: 
     person(given = "Alexander",
            family = "Ilich",
@@ -20,7 +20,10 @@ LinkingTo:
     RcppArmadillo
 Imports: 
     Rcpp,
-    raster
+    raster,
+    future,
+    future.apply,
+    dplyr
 Suggests:
     knitr,
     rmarkdown
diff --git a/NAMESPACE b/NAMESPACE
index d3c1101..30204a8 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -6,6 +6,11 @@ export(make_glcm)
 export(quantize_raster)
 import(terra)
 importFrom(Rcpp,sourceCpp)
+importFrom(dplyr,"%>%")
+importFrom(dplyr,lead)
+importFrom(dplyr,mutate)
+importFrom(future,plan)
+importFrom(future.apply,future_lapply)
 importFrom(raster,raster)
 importFrom(raster,stack)
 importFrom(raster,writeRaster)
diff --git a/R/glcm_textures.R b/R/glcm_textures.R
index 2793d1e..ed3d97f 100644
--- a/R/glcm_textures.R
+++ b/R/glcm_textures.R
@@ -12,6 +12,7 @@
 #' @param maxcell positive integer used to take a regular sample for quantization if "equal prob" is used (default is Inf)
 #' @param na.rm a logical value indicating whether NA values should be stripped before the computation proceeds (default=FALSE)
 #' @param include_scale Logical indicating whether to append window size to the layer names (default = FALSE).
+#' @param ncores Number of cores to use for parallel processing (default is 1 aka serial processing).
 #' @param filename character Output filename. Can be a single filename, or as many filenames as there are layers to write a file for each layer
 #' @param overwrite logical. If TRUE, filename is overwritten (default is FALSE).
 #' @param wopt list with named options for writing files as in writeRaster
@@ -27,14 +28,19 @@
 #' @importFrom raster raster
 #' @importFrom raster stack
 #' @importFrom raster writeRaster
-
+#' @importFrom future plan
+#' @importFrom future.apply future_lapply
 #' @references
 #' Hall-Beyer, M., 2017. GLCM Texture: A Tutorial v. 3.0. University of Calgary, Alberta, Canada.
 #'
 #' Haralick, R.M., Shanmugam, K., Dinstein, I., 1973. Textural features for image classification. IEEE Transactions on Systems, Man, and Cybernetics 610–621. https://doi.org/10.1109/TSMC.1973.4309314
 #' @export
 #'
-glcm_textures<- function(r, w = c(3,3), n_levels, shift=list(c(1,0), c(1,1), c(0,1), c(-1,1)), metrics= c("glcm_contrast", "glcm_dissimilarity", "glcm_homogeneity", "glcm_ASM", "glcm_entropy", "glcm_mean", "glcm_variance", "glcm_correlation", "glcm_SA"), quantization, min_val=NULL, max_val=NULL, maxcell=Inf, na.rm=FALSE, include_scale=FALSE, filename=NULL, overwrite=FALSE, wopt=list()){
+glcm_textures<- function(r, w = c(3,3), n_levels, shift=list(c(1,0), c(1,1), c(0,1), c(-1,1)), metrics= c("glcm_contrast", "glcm_dissimilarity", "glcm_homogeneity", "glcm_ASM", "glcm_entropy", "glcm_mean", "glcm_variance", "glcm_correlation", "glcm_SA"), quantization, min_val=NULL, max_val=NULL, maxcell=Inf, na.rm=FALSE, include_scale=FALSE, ncores=1, filename=NULL, overwrite=FALSE, wopt=list()){
+
+  oplan<- future::plan()
+  on.exit(future::plan(oplan)) #restore original parallelization plan on exit of function
+
   og_class<- class(r)[1]
   if(og_class=="RasterLayer"){
     r<- terra::rast(r) #Convert to SpatRaster
@@ -93,14 +99,39 @@ glcm_textures<- function(r, w = c(3,3), n_levels, shift=list(c(1,0), c(1,1), c(0
     stop("Error: raster must have values between 0 and n_levels-1")}
 
   out_list<- vector(mode = "list", length=length(shift))
-  for (k in 1:length(shift)) {
-    out_list[[k]]<- terra::focalCpp(r, w=w, fun = C_glcm_textures_helper, w2=w, n_levels= n_levels, shift = shift[[k]], metrics = needed_metrics, na_rm=na.rm, fillvalue=NA, wopt=wopt)
-    out_list[[k]]<- terra::subset(out_list[[k]], metrics, wopt=wopt) #Subset from needed to requested metrics
+
+  if(ncores==1){
+    for (k in 1:length(shift)) {
+      out_list[[k]]<- terra::focalCpp(r, w=w, fun = C_glcm_textures_helper, w2=w, n_levels= n_levels, shift = shift[[k]], metrics = needed_metrics, na_rm=na.rm, fillvalue=NA, wopt=wopt)
+      out_list[[k]]<- terra::subset(out_list[[k]], metrics, wopt=wopt) #Subset from needed to requested metrics
+    }} else{
+      buffer<- (w[1]-1)/2
+      r_list<- chunk_raster(r, n_chunks=ncores, buffer=buffer) #Break raster into smaller chunks
+      breaks_df<- r_list[[1]]
+      r_list<- r_list[[2]]
+      nchunks<- length(r_list)
+      for (k in 1:length(shift)) {
+      future::plan(strategy = "multisession", workers=nchunks) #Set up parallel
+      out_list[[k]]<- future.apply::future_lapply(r_list, FUN = focalCpp_parallel,
+                               w=w,
+                               fun = C_glcm_textures_helper,
+                               w2=w,
+                               n_levels= n_levels,
+                               shift = shift[[k]],
+                               metrics = needed_metrics,
+                               na_rm=na.rm,
+                               fillvalue=NA,
+                               wopt=wopt)
+      plan(oplan) #Go back to original plan
+      out_list[[k]]<- lapply(out_list[[k]], terra::unwrap)
+      out_list[[k]]<- combine_raster_chunks(out_list[[k]], breaks_df=breaks_df)
+      out_list[[k]]<- terra::subset(out_list[[k]], metrics, wopt=wopt) #Subset from needed to requested metrics
+      }
     }
 
   n_layers<- length(metrics)
   if(length(shift) > 1){
-    output<- terra::app(terra::sds(out_list), fun=mean, na.rm=na.rm, wopt=wopt)
+    output<- terra::app(terra::sds(out_list), fun=mean, na.rm=na.rm, wopt=wopt) #Could try making this parallel with future package.
     } else{
       output<- out_list[[1]]
     }
diff --git a/R/parallel_code.R b/R/parallel_code.R
new file mode 100644
index 0000000..934afab
--- /dev/null
+++ b/R/parallel_code.R
@@ -0,0 +1,78 @@
+#' FocalCpp for use in parallel code
+#'
+#' FocalCpp for use in parallel code
+#' @param x a packed SpatRaster
+#' @param ... Other arguments for terra::focalCpp
+#' @return a packed SpatRaster
+#' @import terra
+#'
+focalCpp_parallel<- function(x,...){
+  terra::wrap(terra::focalCpp(terra::unwrap(x),...))
+}
+
+#' subset metrics and average across shifts
+#'
+#' subset metrics and average across shifts
+#' @param x a packed SpatRaster
+#' @param shift a list of shifts from glcm_textures
+#' @param metrics Requested GLCM textures metrics to be returned
+#' @param na.rm logical indicating if NA values should be removed from calculations
+#' @param wopt Additional options
+#' @return a packed SpatRaster
+#' @import terra
+#'
+avg_subset_parallel<- function(x,shift,metrics, na.rm, wopt){
+  x<- terra::subset(unwrap(x), metrics, wopt=wopt) #Subset from needed to requested metrics
+  if(length(shift) > 1){
+    output<- terra::app(terra::sds(x), fun=mean, na.rm=na.rm, wopt=wopt)
+  } else{
+    output<- out_list[[1]]
+  }
+  return(terra::wrap(output))
+} #Not yet used. Function in progress for parallelizing averaging of shifts
+
+#' Break raster into smaller chunks
+#'
+#' Break raster into smaller chunks
+#' @param r SpatRaster
+#' @param n_chunks Desired number of chunks to break raster into
+#' @param buffer the number of rows above/below the cell value that the calculation needs access to. For cell by cell calculations this should be zero and for standard focal operations this would be (w-1)/2 where w is the number of rows in the focal window.
+#' @return a list containing a dataframe specifying how rasters are chunked and a list of the chunked SpatRasters
+#' @import terra
+#' @importFrom dplyr mutate
+#' @importFrom dplyr lead
+#' @importFrom dplyr %>%
+
+chunk_raster<- function(r, n_chunks, buffer){
+  breaks<- data.frame(write_start = ceiling(seq(1, nrow(r)+1, length.out = n_chunks+1)))
+  breaks<- breaks[!duplicated(breaks$write_start), ,drop=FALSE] #remove duplicates
+  breaks<- breaks %>% dplyr::mutate(write_end=dplyr::lead(write_start, n=1)-1)
+  breaks<- breaks %>% dplyr::mutate(chunk_start=write_start - buffer)
+  breaks<- breaks %>% dplyr::mutate(chunk_end = write_end + buffer)
+  breaks<- breaks[-nrow(breaks),]
+  breaks$chunk_end[breaks$chunk_end > nrow(r)]<- nrow(r)
+  breaks$chunk_start[breaks$chunk_start < 1]<- 1
+
+  r_list<- vector(mode = "list", length = nrow(breaks))
+  for (i in 1:nrow(breaks)) {
+    r_list[[i]]<- terra::wrap(r[breaks$chunk_start[i]:breaks$chunk_end[i], ,drop=FALSE])
+  }
+  return(list(breaks, r_list))
+}
+
+#' Merge raster chunks from a list into a single raster
+#'
+#' Merge raster chunks from a list into a single raster
+#' @param x A list of SpatRasters
+#' @param breaks_df Dataframe output from chunk_raster
+#' @return a SpatRaster
+#' @import terra
+
+combine_raster_chunks<- function(x, breaks_df){
+  for (i in 1:length(x)) {
+    st<- 1 + breaks_df$write_start[i] - breaks_df$chunk_start[i]
+    end<- nrow(x[[i]]) - (breaks_df$chunk_end[i] - breaks_df$write_end[i])
+    x[[i]]<- x[[i]][st:end, , , drop=FALSE]  # Trim chunk
+  }
+  out<- do.call(terra::merge, x)
+}
diff --git a/man/avg_subset_parallel.Rd b/man/avg_subset_parallel.Rd
new file mode 100644
index 0000000..9e9fb59
--- /dev/null
+++ b/man/avg_subset_parallel.Rd
@@ -0,0 +1,25 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/parallel_code.R
+\name{avg_subset_parallel}
+\alias{avg_subset_parallel}
+\title{subset metrics and average across shifts}
+\usage{
+avg_subset_parallel(x, shift, metrics, na.rm, wopt)
+}
+\arguments{
+\item{x}{a packed SpatRaster}
+
+\item{shift}{a list of shifts from glcm_textures}
+
+\item{metrics}{Requested GLCM textures metrics to be returned}
+
+\item{na.rm}{logical indicating if NA values should be removed from calculations}
+
+\item{wopt}{Additional options}
+}
+\value{
+a packed SpatRaster
+}
+\description{
+subset metrics and average across shifts
+}
diff --git a/man/chunk_raster.Rd b/man/chunk_raster.Rd
new file mode 100644
index 0000000..591bc72
--- /dev/null
+++ b/man/chunk_raster.Rd
@@ -0,0 +1,21 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/parallel_code.R
+\name{chunk_raster}
+\alias{chunk_raster}
+\title{Break raster into smaller chunks}
+\usage{
+chunk_raster(r, n_chunks, buffer)
+}
+\arguments{
+\item{r}{SpatRaster}
+
+\item{n_chunks}{Desired number of chunks to break raster into}
+
+\item{buffer}{the number of rows above/below the cell value that the calculation needs access to. For cell by cell calculations this should be zero and for standard focal operations this would be (w-1)/2 where w is the number of rows in the focal window.}
+}
+\value{
+a list containing a dataframe specifying how rasters are chunked and a list of the chunked SpatRasters
+}
+\description{
+Break raster into smaller chunks
+}
diff --git a/man/combine_raster_chunks.Rd b/man/combine_raster_chunks.Rd
new file mode 100644
index 0000000..e4b13cd
--- /dev/null
+++ b/man/combine_raster_chunks.Rd
@@ -0,0 +1,19 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/parallel_code.R
+\name{combine_raster_chunks}
+\alias{combine_raster_chunks}
+\title{Merge raster chunks from a list into a single raster}
+\usage{
+combine_raster_chunks(x, breaks_df)
+}
+\arguments{
+\item{x}{A list of SpatRasters}
+
+\item{breaks_df}{Dataframe output from chunk_raster}
+}
+\value{
+a SpatRaster
+}
+\description{
+Merge raster chunks from a list into a single raster
+}
diff --git a/man/focalCpp_parallel.Rd b/man/focalCpp_parallel.Rd
new file mode 100644
index 0000000..ae7ccf8
--- /dev/null
+++ b/man/focalCpp_parallel.Rd
@@ -0,0 +1,19 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/parallel_code.R
+\name{focalCpp_parallel}
+\alias{focalCpp_parallel}
+\title{FocalCpp for use in parallel code}
+\usage{
+focalCpp_parallel(x, ...)
+}
+\arguments{
+\item{x}{a packed SpatRaster}
+
+\item{...}{Other arguments for terra::focalCpp}
+}
+\value{
+a packed SpatRaster
+}
+\description{
+FocalCpp for use in parallel code
+}
diff --git a/man/glcm_textures.Rd b/man/glcm_textures.Rd
index 4a0d3ef..135abd0 100644
--- a/man/glcm_textures.Rd
+++ b/man/glcm_textures.Rd
@@ -17,6 +17,7 @@ glcm_textures(
   maxcell = Inf,
   na.rm = FALSE,
   include_scale = FALSE,
+  ncores = 1,
   filename = NULL,
   overwrite = FALSE,
   wopt = list()
@@ -45,6 +46,8 @@ glcm_textures(
 
 \item{include_scale}{Logical indicating whether to append window size to the layer names (default = FALSE).}
 
+\item{ncores}{Number of cores to use for parallel processing (default is 1 aka serial processing).}
+
 \item{filename}{character Output filename. Can be a single filename, or as many filenames as there are layers to write a file for each layer}
 
 \item{overwrite}{logical. If TRUE, filename is overwritten (default is FALSE).}