crop circles

Abstract: The features of the R raster package are demonstrated using a Landsat ETM+ image of Kansas farmland. Basic file manipulation, image display and raster calculations are performed. Linear and histogram equalization contrast enhancements are implemented for visualization improvement. Measuring distances and extracting pixel spectra are also demonstrated

Preliminaries, Contrast Enhancement, Annotation, Interacting with Images, Summary

Cimarron, KS from 705 km

Center pivot irrigation makes farming arid regions possible, and produces a beautifully alien landscape when viewed from 705 km above. Green vegetation appears bright red in the false color composite image because the Landsat near-infrared band is substituted for red in the RGB composite. Healthy chloroplasts reflect most of the near-infrared radiation incident upon them.

All of the code used to generate the figures on this page are included, although particularly long sections may be collapsed by default:

print("Hello World")

Preliminaries
GloVis provides GUI browsing and downloading of images from NASA and USGS administered missions. Data is available from Landsat, EO-1, Aster, MODIS and other sensors. This example uses a Landsat ETM+ image of Kansas (Path:30 Row:34 Date:2003/04/11) downloaded from GloVis. ETM+ has 8 spectral bands: band 6 is a 60 m thermal infrared channel, band 8 is a 15 m panchromatic channel and the remainder are 30 m resolution channels: blue, green, red, nir, swir1, & swir2 (full data details in the Landsat Handbook). This demonstration uses only the 30 m channels and a 52×62 km subset of the full scene. Additionally, we will use the satellite Digital Numbers, DN, instead of converting to radiance or reflectance.

GDAL Utilities are available if GDAL was installed with Python support, and offer an efficient set of tools for basic raster preprocessing such as subsetting and reprojecting. Here’s how I stacked the 30 m multispectral bands and subset in 2 lines from the bash shell:

gdal_merge.py -separate -o kansas -of ENVI `ls|grep B[^68]`
gdalwarp -te 324063 4140016 386707 4192167 -of ENVI kansas kansas_sub

There are many R packages for manipulating geographic data. The sp package offers spatial data classes for R, and the raster package has been developed to offer raster manipulation using sp classes. First we load the raster library, read the ENVI file kansas_sub into a RasterStack and display some image summary information:

library(raster)
etm<-stack("kansas_sub") #read the ENVI file into a RasterStack
etm #print information about etm

The raster package provides several tools for viewing rasters including plot, plotRGB, and image. First, plotRGB is used to display RGB false-color composites for two different image extents. extent is used in two different ways to generate identically sized areas: by min/max lines and samples, and by min/max x,y coordinates. The project function from the rgdal library is a useful tool to reproject coordinate pairs. R can plot to many devices, but we use the jpeg device in order to produce images with file sizes suitable for web display. Inner and outer margins are set to 0 and the device background color is set to black for added contrast. The first extent is plotted as a RGB:421 combination and the second as a RGB:543 composite.

ETM+ 421, and 543 RGB Composites at two different extents

jpeg("SpEx_Extents.jpg", height=6.2, width=12, units="in", res=120, qual=85) #for plotting to a jpg file
par(mfrow=c(1,2), mar=rep(0,4), oma=rep(0,4), bg="black")
extent1<-extent(etm,1,650,1,650) #define extent via max/min lines/samples
plotRGB(etm,4,2,1, extent=extent1, main="RGB:421", axes=F) #plot as a RGB:421 composite

#a different way to subset using extent
ulsub<-c(-100.51,37.86) #define the UL coord using geog coords (x,y not Lat,Long)
dimsub<-c(19500,19500) #the x and y dim of the subset (in m, not pixels)
ulsub.utm<-project(matrix(ulsub,nrow=1), proj=projection(etm)) #project UL in geo coords to dataset projection
extent2<-extent(c(ulsub.utm[1], ulsub.utm[1]+dimsub[1], ulsub.utm[2]-dimsub[2], ulsub.utm[2])) #define extent
plotRGB(etm,5,4,3, extent=extent2, main="RGB:543", axes=F)
dev.off()

Contrast Enhancement
plotRGB does not perform any contrast enhancements, and this results in a muted image with low color range. Digital numbers, DN, are mapped directly to an RGB intensity between 0-255. Visualization of imagery is typically enhanced by stretching the pixel values across the range 0-255, thereby using the full range of intensity. Typical linear stretches include the min/max stretch and the 2% stretch. An R function to perform a linear stretch on a RasterLayer may be written:

# Linear stretch between minmax=c(min,max)
linstretch<-function(img,minmax=NA){
	if(is.na(minmax)) minmax<-c(min(getValues(img),na.rm=T),max(getValues(img),na.rm=T))
	temp<-calc(img,fun=function(x) (255*(x-minmax[1]))/(minmax[2]-minmax[1]))
	#set all values above or below minmax to 0 or 255
	temp[temp<0]<-0;temp[temp>255]<-255;
	return(temp)
}

Another common contrast enhancement is histogram equalization. This is can be easily accomplished utilizing the cumulative distribution function. In R the ecdf function may be used:

# Histogram equalization stretch
eqstretch<-function(img){
	ecdf<-ecdf(getValues(img))
	return(calc(img,fun=function(x) ecdf(x)*255))
}

The effect of alternative contrast enhancements is most clearly demonstrated using a single image band. The stretching functions are used to visualize band 4 subject to four different contrasting techniques: no stretch, linear min/max, linear 2% and histogram equalization. Histograms also reveal the unique qualities of each enhancement:

Comparison of different contrast enhancements on ETM+ band 4

#use crop() to subset etm to extent1
etm.sub<-crop(etm,extent1) #create a RasterStack from the subset
etm.sub.b4<-raster(etm.sub,layer=4) #get *RasterLayer

# 8 panel plot comparing 4 contrast stretches w/ histograms
# uses B4 of etm.sub, custom color ramp
jpeg("Spex_ContrastEnhancement.jpg", height=5.75, width=12, units="in", res=170,qual=85)
nf<-layout(matrix(seq(1,8), byrow=T, nrow=2))
par(mar=c(3,3,1,1),oma=rep(0,4))
myramp<-colorRampPalette(brewer.pal(5,"PiYG")) #custom ramp from RColorBrewer

# no stretch
# BAD hack in order to get plot(*Raster) to look like it didn't apply a Min/Max stretch
temp<-etm.sub.b4; temp[1]<-0; temp[2]<-255;
plot(temp, col=myramp(255)) #plot with diff color map
hist(getValues(etm.sub.b4), breaks=seq(0,255), freq=F, main="0-255")
rm(temp)

#min/max stretch and histogram
b4.minmax<-linstretch(etm.sub.b4)
plot(b4.minmax, col=myramp(255))
hist(getValues(b4.minmax), breaks=seq(0,255), xlim=c(0,255), freq=F, main="Linear Min/Max")
rm(b4.minmax)

#histogram equalization
plot(eqstretch(etm.sub.b4), col=myramp(255))
b4.ecdf<-ecdf(getValues(etm.sub.b4))
hist(getValues(calc(etm.sub.b4, fun=function(x) b4.ecdf(x)*255)), breaks=seq(0,255), xlim=c(0,255), freq=F, main="Histogram Equalization")

#2% linear stretch
b4quant<-quantile(getValues(etm.sub.b4), c(0.02,0.98)) # 2% cutoffs
b4.linstretch<-linstretch(etm.sub.b4, minmax=c(b4quant[1],b4quant[2]))
plot(b4.linstretch, col=myramp(255))
hist(getValues(b4.linstretch), breaks=seq(0,255), xlim=c(0,255), freq=F, main="Linear 2%")
dev.off()

Note that because plot(Raster) automatically performs a min/max contrast enhancement, we must set the value of the first two pixels to 0 and 255, respectively. This gives the appearance of a direct 0-255 mapping without any contrast stretch. The plotRGB function does not behave this way, instead performing no contrast adjustment at all.

Annotation
It is straightforward to apply annotation to the plot using map coordinates. For example, a simple scale-bar can be implemented:

# function to draw a simple scalebar in the lower left
# plot corner that is dim map units long.
scalebar<-function(dim,coord=NA,padding=c(5,5),labels=paste(dim),...){
	#defaults to a lower left hand position
	padding<-padding/100
	parrange<-c(par()$usr[2]-par()$usr[1],par()$usr[4]-par()$usr[3])
	if(is.na(coord)) coord<-c(par()$usr[1]+(padding[1]*parrange[1]),par()$usr[3]+(padding[2]*parrange[2]))
	lines(matrix(c(coord[1],coord[2],coord[1]+dim,coord[2]),byrow=T,nrow=2),...)
	text(coord[1]+(0.5*dim),coord[2],labels=labels,adj=c(0.5,-0.2),...)
}

This function draws a horizontal line starting from coord that is dim map units long. Default behavior is to put the scalebar in the lower left hand corner with padding from the coordinate specified as a percentage of the total plot extent.

Finally, we can use plotRGB to visualize two different contrast enhancements on a false-color composite image, display a simple scale-bar and annotate the plot with text, and axes:

Effect of linear min/max and histogram equalization contrast enhancements on a RGB:421 false-color composite

# Display two false-color composites with different contrast enhancements
# plot a scale bar and annotate each plot
jpeg("SpEx_RGB.jpg",height=6.5,width=12,units="in",res=120,qual=85) #for plotting to a jpg file
# quartz(height=6.5,width=12)
par(mfrow=c(1,2), mar=rep(0,4), oma=c(3,3,3,3), bg="black", col.axis="white", col.lab="white", col.main="white", col.sub="white", col="white")

plotRGB(stack(linstretch(raster(etm.sub,layer=4)), linstretch(raster(etm.sub,layer=2)), linstretch(raster(etm.sub,layer=1))), axes=F);
# add a scalebar
scalebar(5000, padding=c(8,6), labels="5 km", col="white", lwd=1.2, cex=0.8)
# axes and a box
box(col="white")
axis(1, col="white", col.ticks="white"); axis(2, col="white", col.ticks="white")
axis(3, col="white", col.ticks="white", labels=NA); axis(4, col="white", col.ticks="white", labels=NA)
# plot annotation
parrange<-c(par()$usr[2]-par()$usr[1], par()$usr[4]-par()$usr[3])
text(par()$usr[2]-0.01*parrange[1], par()$usr[4]-0.01*parrange[2], labels="Cimarron, KS. ETM+ RGB:421 Max/Min Stretch", adj=c(1,1),col="white", cex=0.65)

plotRGB(stack(eqstretch(raster(etm.sub,layer=4)), eqstretch(raster(etm.sub,layer=2)), eqstretch(raster(etm.sub,layer=1))),axes=F)
scalebar(5000, padding=c(8,6), labels="5 km", col="white", lwd=1.2, cex=0.8)
box(col="white")
axis(1, col="white", col.ticks="white", labels=NA); axis(2, col="white", col.ticks="white", labels=NA)
axis(3, col="white", col.ticks="white"); axis(4, col="white", col.ticks="white")
parrange<-c(par()$usr[2]-par()$usr[1], par()$usr[4]-par()$usr[3])
text(par()$usr[2]-0.01*parrange[1], par()$usr[4]-0.01*parrange[2], labels="UTM14N WGS84", adj=c(1,1),col="white", cex=0.65)
dev.off()

Interacting with Images
raster provides some tools for interacting with a plotted RasterLayer on the screen. For instance, click returns the raster values, and optionally x,y locations. This can be used in conjunction with pointDistance to measure image distances on the screen. For example, we can plot a false-color image with histogram equalization stretch, measure the distance between two points, and extract and plot spectra from 4 points:

Measuring distance and extracting spectra from RasterLayers plotted to the screen

extent3<-extent(etm,300,600,300,600) #define extent via max/min lines/samples
etm.sub3<-crop(etm,extent3)
quartz(height=5, width=10)
par(mfrow=c(1,2), mar=rep(0,4), oma=rep(0,4))
plotRGB(stack(eqstretch(raster(etm.sub3, layer=4)), eqstretch(raster(etm.sub3, layer=2)), eqstretch(raster(etm.sub3, layer=1))), axes=F)

#Find the distance between two points supplied by user
clicks<-click(etm.sub3, n=2, xy=T, type="o") #get 2 clicks in plot() window
dist<-pointDistance(clicks[1,1:2], clicks[2,1:2]) #calculate distance
text(mean(c(clicks[1,1], clicks[2,1])), mean(c(clicks[1,2], clicks[2,2])), labels=round(dist), cex=0.8, adj=c(0.5,1)) #annotate line with distance

numclicks<-4 #number of spectra to query and plot
clicks<-click(etm.sub3, n=numclicks, type="p", xy=T) #get clicks
points(clicks[,1], clicks[,2], pch=21, lwd=2, cex=1.2, col=1, bg="white", pt.cex=0.6) #add points to screen
text(clicks[,1], clicks[,2], labels=1:dim(clicks)[1], lwd=2, cex=1.3, col="black", adj=c(0.5,-0.5)) #label points
bandCenterNM<-c(mean(c(450,515)), mean(525,605), mean(630,690), mean(775,900), mean(1550,1750), mean(2090,2350)) #etm band centers (nm)
par(mar=c(3,4,1,1), oma=rep(0,4), new=T)
plot(bandCenterNM,rep(255,  length(bandCenterNM)), type="n", ylim=c(min(clicks[,4:dim(clicks)[2]])*0.98, max(clicks[,4:dim(clicks)[2]])*1.02), xlab="Wavelength (nm)", ylab="DN") #plot the spectra
for(i in 1:dim(clicks)[1]) lines(bandCenterNM, clicks[i,3:dim(clicks)[2]], type="b", col=i, lty=i)
legend("topleft", legend=1:dim(clicks)[1], bty="n", col=1:dim(clicks)[1], lty=1:dim(clicks)[1])

Summary: The raster library brings basic remote sensing functionality to R. Reading and displaying images from a variety of file formats is straightforward. So is extracting data and performing calculations and analyses. Plotting is a bit clunky and there are no built-in contrast enhancements, though greater efficiency can likely be achieved using the maxpixels argument to plot and plotRGB, and contrast enhancements are simple to implement. We demonstrated simple linear and histogram equalization stretches. It could be argued that a color palette remapping would be preferable to an image value remapping, but the increase in efficiency is likely small. With respect to maxpixels it should probably be set based on the par() settings and the device resolution. Annotation of raster plots is accomplished using the map coordinates, and this can be used to add points, lines, and information such as scalebars. Functions for interacting with an image plotted to the screen can be used to extract pixel values and locations. Overall, raster provides a solid set of foundational remote sensing tools that advanced practitioners can use to accomplish most remote sensing tasks. However, the visualization capabilities offered by plot, plotRGB, and image are not likely to replace GUI tools like ENVI, Imagine or ArcMap for data exploration.

Tagged: GDAL, landsat, R, raster, remote sensing

Posted in: R, remote sensing, tutorial