The maptools package has a pruneMap() function t0 crop map objects in R. In practice, the function extracts data from SpatialPolygon or SpatialLine objects given a boundary box or specific area of interest. Unfortunately, there is no equivalent function for high resolution, large data, raster images, which are common in many Earth Science applications. The following post defines a custom function to crop raster images in R and to extract data from SpatialGridDataFrames. The function is tested using a raster image from the Shuttle Radar Topography Mission (SRTM; shown at left). The resulting data is then mapped using the image() function in R.
Function Data, Use and Definition
Digitial elevation data from the SRTM is packaged by the U.S. Geological Survey as GOTOPO30 data. The pruneDEM() function is used to ingest GOTOPO30 data into R given a path and file name. The function then extracts data using a boundary box defined by min/max coordinates. The result is a SpatialGridDataFrame. that can be used as a land surface definition for GIS mapping, engineering design, or spatial weather modeling.
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 |
pruneDEM <- function(x, xlim = NULL, ylim = NULL, ...) { require(rgdal) # Basic data object diagnostics GTOPOinfo <- GDALinfo(x) DEMoffset <- GTOPOinfo[c("ll.x", "ll.y")] DEMscale <- GTOPOinfo[c("res.x", "res.y")] DEMdim <- GTOPOinfo[c("columns", "rows")] # Define area of interest and boundary box Xseq <- seq(DEMoffset[1], by = DEMscale[1], length = DEMdim[1]) + DEMscale[1]/2 Yseq <- rev(seq(DEMoffset[2], by = DEMscale[2], length = DEMdim[2]) + DEMscale[2]/2) # xlim condition testing if (!is.null(xlim)) { if (!is.numeric(xlim)) stop("xlim must be numeric") if (!length(xlim) == 2) stop("xlim must be of length 2") if (!diff(xlim) > 0) stop("xlim[1] must be less than xlim[2]") Xindex <- which(Xseq >= xlim[1] & Xseq <= xlim[2]) } # ylim condition testing if (!is.null(ylim)) { if (!is.numeric(ylim)) stop("ylim must be numeric") if (!length(ylim) == 2) stop("ylim must be of length 2") if (!diff(ylim) > 0) stop("ylim[1] must be less than ylim[2]") Yindex <- which(Yseq >= ylim[1] & Yseq <= ylim[2]) } # Confirm boundary box bbox <- matrix(c(Xseq[Xindex[1]], Yseq[Yindex[1]], Xseq[Xindex[length(Xindex)]], Yseq[Yindex[length(Yindex)]]), 2, 2) # Finalize readGDAL pruning parameters rgdal.offset <- rev(c(min(Xindex), min(Yindex))) rgdal.dim <- rev(c(length(Xindex)-1, length(Yindex))) # Ingest and prune DEM data given xlim and ylim readGDAL(x, offset = rgdal.offset, region.dim = rgdal.dim) } |
Working with Raster Images in R
The following example ingests a large tile of digital elevation data, which can be downloaded below. pruneDEM() will ingest the raster file and extract data given a boundary box. Function output is then plotted with custom colors and axes legends.
[sdm-download id=”3769″ fancy=”0″]
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 |
# Boundary box coordinates lon <- c(28, 60) lat <- c(14, 35) lon <- c(28, 60) lat <- c(14, 35) bbox <- matrix(c(lon[1], lat[1], lon[2], lat[2]), 2, 2, dimnames = list(c("lon", "lat"), c("min", "max"))) # Prune digital elevation map dem.data <- "/Users/bjh/@Data/SpatialData/GTOPO30/e020n40/E020N40.DEM" map.tile <- pruneDEM(dem.data, xlim=lon, ylim=lat) # Plot DEM M0title <- "Digital Elevation Map" M1title <- "Shuttle Radar Topography Mission" S0title <- "Fig. 1 Custom boundary box defined by pruneDEM()" S1title <- "GOTOPO30 Data; USGS Earth Resources Observation and Science Center" dem.colors <- colorRampPalette(c("papayawhip", "bisque3", "orchid2", "red4")) par(mar = c(5, 4, 5, 4) + 0.1, pin=c(6.5, 5.6)) image(map.tile, col = dem.colors(12)) axis(1, at=c(28 + 0:8 * 4), pos=lat[1], las=1, cex.axis = 0.5, col = "grey45", col.axis = "grey45", labels=parse(text=sp:::degreeLabelsEW(c(28 + 0:8 * 4)))) axis(2, at=c(14 + 0:7 * 3), pos=lon[1], las=2, cex.axis = 0.5, col = "grey45", col.axis = "grey45", labels=parse(text=sp:::degreeLabelsNS(c(14 + 0:7 * 3)))) axis(3, at=c(28 + 0:8 * 4), pos=lat[2], las=1, cex.axis = 0.5, col = "grey45", col.axis = "grey45", labels=parse(text=sp:::degreeLabelsEW(c(28 + 0:8 * 4)))) axis(4, at=c(14 + 0:7 * 3), pos=lon[2], las=2, cex.axis = 0.5, col = "grey45", col.axis = "grey45", labels=parse(text=sp:::degreeLabelsNS(c(14 + 0:7 * 3)))) title(main=M0title, font.main = 1, line = 2) title(main=M1title, cex.main = 0.85, font.main = 1.75, line = 0.75) title(sub=S0title, cex.sub = 0.85, line = 0.75) title(sub=S1title, cex.sub = 0.7,font.sub = 1.3, line = 2.25) |
