<html style="direction: ltr;">
<head>
<meta http-equiv="Content-Type" content="text/html; charset=UTF-8">
<style id="bidiui-paragraph-margins" type="text/css">body p { margin-bottom: 0cm; margin-top: 0pt; } </style>
</head>
<body bidimailui-charset-is-forced="true" style="direction: ltr;">
<p><font size="4">Hello</font></p>
<p><font size="4"></font><br>
</p>
<div class="moz-cite-prefix">On 16/03/2022 12:20, Philip Monk wrote:<br>
</div>
<blockquote type="cite"
cite="mid:CAMVDa_7uC98Qs93-Z=W=m1_PXBK+CFdD9-Lu=dBcedeyAr65zg@mail.gmail.com">
<pre class="moz-quote-pre" wrap="">I would like to create two stacks of raster images, each containing 12
images in total (one for each month over a one-year time period), then
process each image to:
*) display pixels which are 1/2/3 standard deviations from the mean (I
would amend and run the code several times to plot different values -
not all in one go)
*) add a mask to indicate the aoi
*) crop the image to the aoi
*) save each raster as an individual file.
The images are geoTIFFs created in and exported from Google Earth
Engine (GEE). Each pixel is a calculated value of Land Surface
Temperature derived from Landsat 8 (Collection 2) scenes. They are
identical in extent et al.
The mask is a shapefile exported from GEE. it has the same extent but
a different crs. The code below addresses that, but I have also used
QGIS to dissolve its geometry to create an outline (for the purposes
of indicating the aoi on the mapped data as its structure is not
relevant for this analysis).
</pre>
</blockquote>
<p><br>
</p>
<p>It's not clear why you are not just completing the process in
GEE, but see below.</p>
<p><br>
</p>
<blockquote type="cite"
cite="mid:CAMVDa_7uC98Qs93-Z=W=m1_PXBK+CFdD9-Lu=dBcedeyAr65zg@mail.gmail.com">
<pre class="moz-quote-pre" wrap="">
At present I can process an individual image:
####
library(raster)
library(rgdal)</pre>
</blockquote>
<p><br>
</p>
<p>First, it's time to switch to terra</p>
<p><br>
</p>
<blockquote type="cite"
cite="mid:CAMVDa_7uC98Qs93-Z=W=m1_PXBK+CFdD9-Lu=dBcedeyAr65zg@mail.gmail.com">
<pre class="moz-quote-pre" wrap="">
LST <- raster("temperature.tif")</pre>
</blockquote>
<p><br>
</p>
<p>You can read all your *.tif into a SpatRaster (stacked) with the
terra package, by passing a list of file names to the rast()
function. I see you did this below (using the raster package)...<br>
</p>
<p><br>
</p>
<blockquote type="cite"
cite="mid:CAMVDa_7uC98Qs93-Z=W=m1_PXBK+CFdD9-Lu=dBcedeyAr65zg@mail.gmail.com">
<pre class="moz-quote-pre" wrap="">
Mask <- readOGR("mask.shp")
# For mask not modified in QGIS
crs(LST)
crs(Mask)
lst = projectRaster(LST, crs=crs(Mask))
plot(lst)
plot(Mask, add=TRUE)
masked <- mask(x = lst, mask = Mask, inverse=FALSE)
plot(masked)
cropped <- crop(masked, Mask)
plot(cropped)
</pre>
</blockquote>
<p><br>
</p>
<p>I think this is backwards: crop first to limit the bounding box
to the Mask polygon, then apply the mask function.</p>
<p><br>
</p>
<blockquote type="cite"
cite="mid:CAMVDa_7uC98Qs93-Z=W=m1_PXBK+CFdD9-Lu=dBcedeyAr65zg@mail.gmail.com">
<pre class="moz-quote-pre" wrap="">
# Calculate mean and sd
sdVal <- cellStats(cropped, "sd")
meanVal <- cellStats(cropped, "mean")
s <- cropped
s[which(s[] > (meanVal-sdVal)] <- NA # for example
plot(s)
plot(Mask, add=TRUE)
####
I can also create a stack:
####
#Load libraries
library(raster)
# List files
files <- list.files("...", pattern = "tif", full.names = TRUE)
# Create stack
siteStack <- raster::stack(files)</pre>
</blockquote>
<p>Prefer:<font size="5" face="monospace"> terra::rast(files)</font><br>
</p>
<p>You can now apply the crop and mask to the whole stack at one go.</p>
<p><br>
</p>
<blockquote type="cite"
cite="mid:CAMVDa_7uC98Qs93-Z=W=m1_PXBK+CFdD9-Lu=dBcedeyAr65zg@mail.gmail.com">
<pre class="moz-quote-pre" wrap="">
####
I cannot, however, work out how to iterate through the stack, and save
individual tifs.
One suggestion I've come across is to use a for loop:
####
# List all the files you want to work on
raster.files = dir(“/myfolder/my.raster.folder”)
for (i in 1:length(raster.files)){
map = raster(raster.files[i]
# do some other stuff …
# Save results as csv or write.raster() to save maps
}
####
I've seen elsewhere that you should avoid for loops and use lapply().</pre>
</blockquote>
<p><br>
</p>
<p>Here is how I would do it:</p>
<p>(Using some fictitious rasters, prepared for this example)</p>
<p><br>
</p>
<p><font size="5" face="monospace">library(terra)<br>
# Prepare a few fictitious rasters<br>
# Terra comes with volcano data, convert to SpatRaster<br>
r1 <- rast(volcano)<br>
# Second SpatRaster with normalized values<br>
r2 <- scale(r1)<br>
# Third as Kmeans cluster<br>
km <- kmeans(r1, centers=10, iter.max=100)<br>
r3 <- rast(r1, vals=km$cluster)<br>
<br>
# Stack<br>
rStk <- c(r1, r2, r3)<br>
nlyr(rStk)<br>
<br>
# Now use lapply to loop over layers<br>
rList <- lapply(1:nlyr(rStk), function(i) {<br>
r = rStk[[i]]<br>
# Get cutoff value<br>
sdVal <- as.numeric(global(r, "sd"))<br>
meanVal <- as.numeric(global(r, "mean"))<br>
lowIndex <- which(values(r) <= meanVal + sdVal)<br>
values(r)[lowIndex] <- NA<br>
file_name = paste0("above_std_", as.character(i), ".tif")<br>
writeRaster(r, file_name)<br>
return(r)<br>
})<br>
# Recreate the stack:<br>
rAboveSd = rast(rList)<br>
</font></p>
<p><br>
</p>
<p><br>
</p>
<p>HTH</p>
<p><br>
</p>
<blockquote type="cite"
cite="mid:CAMVDa_7uC98Qs93-Z=W=m1_PXBK+CFdD9-Lu=dBcedeyAr65zg@mail.gmail.com">
<pre class="moz-quote-pre" wrap="">
Either way, though I can find plenty of examples that process the
stack using raster::calc() to create a single output file, I can't get
anything to work that sequentially processes each image in isolation.
Stackoverflow, R4DS, the usual websites and blogs, even my
University's internal R help Slack group have all drawn a blank.
(Which I find a bit surprising - this doesn't feel like that
complicated or an unusual thing to do).
I've never written anything like this in R, so any help would be very welcome.
Thanks for your time,
Philip
Part-time PhD student, Lancaster Environment Centre, Lancaster University.
_______________________________________________
R-sig-Geo mailing list
<a class="moz-txt-link-abbreviated" href="mailto:R-sig-Geo@r-project.org">R-sig-Geo@r-project.org</a>
<a class="moz-txt-link-freetext" href="https://stat.ethz.ch/mailman/listinfo/r-sig-geo">https://stat.ethz.ch/mailman/listinfo/r-sig-geo</a>
</pre>
</blockquote>
<pre class="moz-signature" cols="72">--
Micha Silver
Ben Gurion Univ.
Sde Boker, Remote Sensing Lab
cell: +972-523-665918</pre>
</body>
</html>