How do I actually load these massive satellite images without crashing my R session?

Look, I've spent way too many hours staring at a frozen screen because I tried to load a multi-gigabyte GeoTIFF directly into memory. The secret is to stop using the old raster package and move entirely to terra. It's significantly faster and, more importantly, it doesn't actually load the whole image into your RAM immediately; it just creates a pointer to the file on your disk.

library(terra)

# Instead of loading everything, terra just reads the metadata
# Let's assume we have a Sentinel-2 composite image
img <- rast("sentinel2_composite.tif")

# You can check the dimensions and resolution without a memory spike
print(img)

If you're working with a "stack" of bands (like a typical remote sensing product), terra treats this as a SpatRaster object where each layer is a different wavelength. Just be careful when you start doing heavy calculations; that's when R will actually start pulling data into memory.

How do I calculate a vegetation index like NDVI using specific bands?
Once you have your raster stack, you're basically doing matrix algebra across pixels. For the Normalized Difference Vegetation Index (NDVI), we need the Red and Near-Infrared (NIR) bands. In a standard Sentinel-2 image, that's usually Band 4 (Red) and Band 8 (NIR). 
The beauty of terra is that it supports vectorized arithmetic. You don't need to write a loop (please, for the love of your CPU, don't write a loop over pixels in R).
# Extract the specific bands
# Assuming Band 4 is Red and Band 8 is NIR
red <- img[[4]]
nir <- img[[8]]

# The NDVI formula: (NIR - Red) / (NIR + Red)
ndvi <- (nir - red) / (nir + red)

# Let's give it a proper name so we don't lose track
names(ndvi) <- "NDVI"

# Quick check to see the result
plot(ndvi, col=rev(terrain.colors(10)), main="Vegetation Index")

I usually suggest checking the range of your results immediately. NDVI should be between -1 and 1. If you see values outside that, you've probably got some "no-data" pixels or outliers that need masking.

My image is huge, but I only care about a small study area. How do I clip it?
You don't want to run analysis on the entire tile if you only need a few hectares. You'll need an sf object (a shapefile or GeoJSON) that defines your area of interest (AOI). There are two steps here: crop, which cuts the rectangular bounding box, and mask, which sets everything outside your actual polygon to NA.
library(sf)

# Load your study area boundary
aoi <- st_read("study_area.shp")

# Ensure the CRS matches! This is where 90% of remote sensing errors happen.
# We project the AOI to match the image's coordinate system.
aoi <- st_transform(aoi, crs(img))

# 1. Crop to the bounding box (fast)
cropped_img <- crop(img, aoi)

# 2. Mask to the exact polygon shape (precise)
final_img <- mask(cropped_img, aoi)

plot(final_img[[1]])
plot(st_geometry(aoi), add=TRUE, border="red")

I always perform the crop before the mask. Cropping reduces the number of pixels the computer has to process during the masking phase, which saves you a lot of time when working with high-resolution imagery.

📋 Practical Task

Exercise: Processing a Drought Stress Map for the Napa Valley Vineyard

You have been provided with a multi-band GeoTIFF napa_vineyard.tif containing 10 spectral bands and a shapefile vineyard_boundary.shp outlining the specific vineyard plots.

Write a script to perform the following:

• Load the raster and the shapefile using terra and sf.
• Ensure the shapefile is projected to the same CRS as the raster.
• Crop and mask the raster to the vineyard boundary.
• Calculate the NDVI using Band 4 (Red) and Band 8 (NIR).
• Identify the "stressed" areas by creating a new binary raster where pixels with an NDVI value below 0.3 are marked as 1 (stressed) and all others are 0.
• Save the final stressed-area map as a new GeoTIFF named vineyard_stress_map.tif.