R/SpatialGrid-methods.R

Defines functions plot.SpatialPixels plot.SpatialGrid length.SpatialGrid length.SpatialPixels IDvaluesSpatialPixels as.SpatialPolygons.SpatialPixels print.SpatialGrid print.SpatialPixels as.data.frame.SpatialPixels gridIndex2nb rcFromGridIndex getGridIndex gridparameters areaSpatialGrid getGridTopology coordnamesSG row.names.SpatialGrid row.names.SpatialPixels SpatialGrid SpatialPixels

Documented in areaSpatialGrid as.data.frame.SpatialPixels as.SpatialPolygons.SpatialPixels getGridIndex getGridTopology gridIndex2nb gridparameters IDvaluesSpatialPixels plot.SpatialGrid SpatialGrid SpatialGrid SpatialPixels

SpatialPixels = function(points, tolerance = sqrt(.Machine$double.eps), 
		proj4string = CRS(as.character(NA)), round = NULL, grid = NULL) {
	if (!is(points, "SpatialPoints"))
		stop("points should be of class or extending SpatialPoints")
	is.gridded = gridded(points)
	if (is.na(slot(slot(points, "proj4string"), "projargs"))) 
		proj4string(points) = proj4string
	if (is.null(grid))
		grid = points2grid(points, tolerance, round)
	if (!is.gridded) {
		points@bbox[,1] = points@bbox[,1] - 0.5 * grid@cellsize
		points@bbox[,2] = points@bbox[,2] + 0.5 * grid@cellsize
	}
	new("SpatialPixels", points, grid = grid, 
		grid.index = getGridIndex(coordinates(points), grid))
}

SpatialGrid = function(grid, proj4string = CRS(as.character(NA))) {
	stopifnot(is(grid, "GridTopology"))
	if (is.character(proj4string))
		proj4string = CRS(proj4string)
	cc = rbind(grid@cellcentre.offset - 0.5 * grid@cellsize,
			grid@cellcentre.offset + (grid@cells.dim - 0.5) * grid@cellsize)
	bb = .bboxCoords(cc)
	new("SpatialGrid", grid = grid, bbox = bb, proj4string = proj4string)
}

setAs("SpatialGrid", "SpatialPixels", 
	function(from) { 
		pts = as(from, "SpatialPoints")
		pts@bbox[,1] = pts@bbox[,1] - 0.5 * from@grid@cellsize
		pts@bbox[,2] = pts@bbox[,2] + 0.5 * from@grid@cellsize
		new("SpatialPixels", pts, grid = from@grid, grid.index = 1:length(pts))
	}
)

setAs("SpatialGrid", "SpatialPoints", 
	function(from) SpatialPoints(coordinates(from), from@proj4string))

setMethod("coordinates", "SpatialPixels", function(obj) obj@coords)

row.names.SpatialPixels <- function(x) {
    ret = dimnames(slot(x, "coords"))[[1]]
	if (is.null(ret))
		ret = seq_len(nrow(x@coords))
	ret
}

row.names.SpatialGrid <- function(x) 1:prod(x@grid@cells.dim)

setMethod("coordinates", "SpatialGrid", function(obj) coordinates(obj@grid))

setMethod("plot", signature(x = "SpatialGrid", y = "missing"),
    function(x, y, ..., grid = TRUE) {
		if (grid)
			plot.SpatialGrid(x, ...)
		else
			plot(as(x, "SpatialPoints"), ...)
	}
)

setMethod("plot", signature(x = "SpatialPixels", y = "missing"),
    function(x, y, ..., grid = TRUE) {
		if (grid)
			plot.SpatialPixels(x, ...)
		else
			plot(as(x, "SpatialPoints"), ...)
	}
)

coordnamesSG = function(x, value) {
	dimnames(x@bbox)[[1]] = value
	if (is(x, "SpatialPixels"))
		dimnames(x@coords)[[2]] = value
	coordnames(x@grid) = value
	x
}

setReplaceMethod("coordnames", 
	signature(x = "SpatialGrid", value = "character"), coordnamesSG)
setReplaceMethod("coordnames", 
	signature(x = "SpatialPixels", value = "character"), coordnamesSG)

getGridTopology = function(obj) {
	if (!(is(obj, "SpatialPixels") || is(obj, "SpatialGrid")))
		stop("object is not or does not extend class SpatialPixels or SpatialGrid")
	obj@grid
}

areaSpatialGrid = function(obj) {
	cellarea = prod(obj@grid@cellsize)
	if (is(obj, "SpatialGrid"))
		return(prod(obj@grid@cells.dim) * cellarea)
	else # take number of cells:
		length(obj@grid.index) * cellarea
}

gridparameters = function(obj) { 
	if (is(obj, "SpatialPixels") || is(obj, "SpatialGrid"))
		obj = obj@grid
	if (is(obj, "GridTopology"))
		return(data.frame(
			cellcentre.offset = obj@cellcentre.offset,
			cellsize = obj@cellsize,
			cells.dim = obj@cells.dim))
	return(numeric(0))
}

getGridIndex = function(cc, grid, all.inside = TRUE) {
	n = ncol(cc)
	idx = rep(1, nrow(cc))
	cumprod = 1
	for (i in 1:n) {
		this.idx = round((cc[,i] - grid@cellcentre.offset[i])/grid@cellsize[i])
		if (i == 2)
			this.idx = grid@cells.dim[2] - (this.idx + 1)
		outside = this.idx >= grid@cells.dim[i] | this.idx < 0
		if (any(outside)) {
			if (all.inside) {
				print(summary(this.idx))
				stop("this.idx out of range")
			} else
				this.idx[outside] = NA
		}
		idx = idx + this.idx * cumprod
		cumprod = cumprod * grid@cells.dim[i]
	}
	outside = idx < 1 | idx > .NumberOfCells(grid)
	if (any(na.omit(outside))) {
		print(summary(idx))
		stop("index outside boundaries")
	}
	as.integer(round(idx))
}

rcFromGridIndex = function(obj) { # returns <x,y> entries: col,row rather than row,col
	obj = as(obj, "SpatialPixels")
	gi = obj@grid.index
	grid = obj@grid
	stopifnot(ncol(coordinates(obj)) == 2)
	xi = ((gi - 1) %% grid@cells.dim[1]) + 1
	yi = grid@cells.dim[2] - ((gi - 1) %/% grid@cells.dim[1])
	cbind(xi,yi)
}

gridIndex2nb = function(obj, maxdist = sqrt(2), fullMat = TRUE, ...) {
	xy = rcFromGridIndex(obj)
	if (fullMat)
		lst = apply(as.matrix(dist(xy, ...)), 1, function(x) which(x <= maxdist))
	else {
		dst = function(X, Xi) apply(X, 1, function(Y) sqrt(sum((Y-Xi)^2)))
		n = nrow(xy)
		lst = vector(mode = "list", length = n)
		for (i in 1:n) { # avoid the n x n matrix construction:
			d = dst(xy, xy[i,])
			lst[[i]] = which(d <= maxdist)
		}
	}
	class(lst) = c("nb", "list")
	lst
}

setMethod("[", "SpatialPixels",
	function(x, i, j, ..., drop = FALSE) {
		if (!missing(j))
			stop("can only select pixels with a single index")
		if (missing(i))
			return(x)
		if (is(i, "Spatial"))
			i = !is.na(over(x, geometry(i)))
		if (drop) { # if FALSE: adjust bbox and grid
			res = as(x, "SpatialPoints")[i]
			tolerance = list(...)$tolerance
			if (!is.null(tolerance))
				res = SpatialPixels(res, tolerance = tolerance)
			else
				gridded(res) = TRUE
			res
		} else # default: don't adjust bbox and grid
			new("SpatialPixels", bbox = x@bbox, 
				proj4string = x@proj4string,	
				coords = x@coords[i, , drop = FALSE], grid = x@grid, 
				grid.index = x@grid.index[i])
	}
)

setMethod("[", "SpatialGrid",
	function(x, i, j, ..., drop = TRUE) {
		drop <- FALSE
#		if (!missing(drop))
#			stop("don't supply drop: it needs to be FALSE anyway")
		gr = x@grid
		if (missing(i))
			rows = 1:gr@cells.dim[2]
		else {
			if (is(i, "Spatial"))
				stop("area selection only makes sense for objects of class SpatialPixels or SpatialGridDataFrame; for object of class SpatialGrid you can only select x[rows,cols]")
			rows = i
		}
		if (missing(j))
			cols = 1:gr@cells.dim[1]
		else
			cols = j
		idx = 1:prod(gr@cells.dim[1:2])
		if (any(is.na(rows)) || any(is.na(cols))) 
			stop("NAs not permitted in indices")
		m = matrix(idx, gr@cells.dim[2], gr@cells.dim[1], byrow = TRUE)[rows,cols]
		idx = as.vector(m) # t(m)?
		cc = SpatialPixels(SpatialPoints(coordinates(x)[idx,,drop=FALSE], rebuild_CRS(slot(x, "proj4string"))))
		cc = as(cc, "SpatialGrid")
		cc
	}
)

setAs("SpatialPixels", "SpatialGrid", function(from) SpatialGrid(from@grid, from@proj4string))
#setAs("SpatialGrid", "SpatialPixels", 
#	function(from) {
#		pts = new("SpatialPoints", coords = coordinates(from),
#			bbox = from@bbox, proj4string = from@proj4string)
#		new("SpatialPixels", pts, grid = from@grid, grid.index = 1:NROW(cc))
#})

as.data.frame.SpatialPixels = function(x, row.names, optional, ...)
	as.data.frame(coordinates(x))

as.data.frame.SpatialGrid = as.data.frame.SpatialPixels

setAs("SpatialPixels", "data.frame", 
	function(from) as.data.frame.SpatialPixels(from))
setAs("SpatialGrid", "data.frame", 
	function(from) as.data.frame.SpatialGrid(from))

print.SpatialPixels = function(x, ...) {
	cat("Object of class SpatialPixels\n")
	print(summary(x@grid))
	cat("SpatialPoints:\n")
	print(coordinates(x))
	pst <- paste(strwrap(paste(
		"Coordinate Reference System (CRS) arguments:", 
		proj4string(x))), collapse="\n")
	cat(pst, "\n")
	invisible(x)
}
setMethod("show", "SpatialPixels", function(object) print.SpatialPixels(object))

print.SpatialGrid = function(x, ...) {
	cat("Object of class SpatialGrid\n")
	print(summary(x@grid))
	cat("SpatialPoints:\n")
	print(coordinates(x))
	pst <- paste(strwrap(paste(
		"Coordinate Reference System (CRS) arguments:", 
		proj4string(x))), collapse="\n")
	cat(pst, "\n")
	invisible(x)
}
setMethod("show", "SpatialGrid", function(object) print.SpatialGrid(object))

# make a SpatialPolygons from a SpatialPixels - Kohris Sahlen workshop
as.SpatialPolygons.SpatialPixels <- function(obj) {
	obj_crds <- coordinates(obj)
	IDs <- IDvaluesSpatialPixels(obj)
	nPolygons <- nrow(obj_crds)
	cS <- slot(slot(obj, "grid"), "cellsize")
	cS2 <- cS/2
	cS2x <- cS2[1]
	cS2y <- cS2[2]
	Srl <- vector(mode="list", length=nPolygons)
	for (i in 1:nPolygons) {
		xi <- obj_crds[i,1]
		yi <- obj_crds[i,2]
		x <- c(xi-cS2x, xi-cS2x, xi+cS2x, xi+cS2x, xi-cS2x)
		y <- c(yi-cS2y, yi+cS2y, yi+cS2y, yi-cS2y, yi-cS2y)
		Srl[[i]] <- Polygons(list(Polygon(coords=cbind(x, y))), ID=IDs[i])
                comment(Srl[[i]]) <- "0"
	}
	res <- SpatialPolygons(Srl, proj4string=rebuild_CRS(slot(obj, "proj4string")))
	res
}
setAs("SpatialPixels", "SpatialPolygons", 
	function(from) as.SpatialPolygons.SpatialPixels(from))

IDvaluesSpatialPixels <- function(obj) {
	if (!is(obj, "SpatialPixels"))
		stop("function only works for objects of class or extending SpatialPixels")

	cc <- slot(obj, "grid.index")
	res <- as.matrix(sapply(cc, as.integer))
	paste("g", res, sep="")
}

length.SpatialPixels = function(x) { nrow(x@coords) }
length.SpatialGrid = function(x) { .NumberOfCells(x@grid) }

setAs("SpatialGrid", "SpatialPolygons", function(from) {
		ret = as.SpatialPolygons.GridTopology(from@grid)
		proj4string(ret) = proj4string(from)
		ret
	}
)
setMethod("coordnames", signature(x = "SpatialGrid"),
	function(x) dimnames(bbox(x))[[1]])

setReplaceMethod("coordnames", signature(x = "SpatialGrid", 
	value =  "character"),
    function(x, value) {
		dimnames(x@bbox)[[1]] = value
		coordnames(x@grid) = value
		x
	}
)

setAs("SpatialGrid", "GridTopology", function(from) getGridTopology(from))
setAs("SpatialPixels", "GridTopology", function(from) getGridTopology(from))

plot.SpatialGrid = function(obj, ..., col = par("fg"), 
		lty = par("lty"), lwd = par("lwd"), add = FALSE) {

  if (! add)
  	plot(as(obj, "Spatial"), ...)
  gr = obj@grid
  # Don MacQueen, Feb 12 2015
  csiz <- gr@cellsize
  ncells <- gr@cells.dim
  nbounds <- ncells+1
  
  ## first get and sort the cell centers
  cv <- coordinatevalues(gr)
  cv[[1]] <- sort(cv[[1]])
  cv[[2]] <- sort(cv[[2]])

  ## calculate cell boundaries
  cv[[1]] <- c(cv[[1]][1] - csiz[1]/2,  cv[[1]] + csiz[1]/2 )
  cv[[2]] <- c(cv[[2]][1] - csiz[2]/2,  cv[[2]] + csiz[2]/2 )
  
  ## construct endpoints of cell boundary lines
  ## vertical lines
  vfrom <- cbind(cv[[1]], cv[[2]][1])
  vto   <- cbind(cv[[1]], cv[[2]][nbounds[2]])

  ## horizontal lines
  hfrom <- cbind(cv[[1]][1], cv[[2]])
  hto   <- cbind(cv[[1]][nbounds[1]], cv[[2]])
    
  ## add to plot
  segments(vfrom[,1], vfrom[,2] , vto[,1], vto[,2], 
  	col = col, lty = lty, lwd = lwd)
  segments(hfrom[,1], hfrom[,2] , hto[,1], hto[,2],
  	col = col, lty = lty, lwd = lwd)
}

plot.SpatialPixels = function(obj, ..., col = par("fg"), 
		lty = par("lty"), lwd = par("lwd"), add = FALSE) {

	# based on plot.SpatialGrid:
	if (! add)
		plot(as(obj, "Spatial"), ...)

	gr = obj@grid
	csiz <- gr@cellsize
	ncells <- gr@cells.dim
	nbounds <- ncells + 1

	m = matrix(FALSE, ncells[2], ncells[1])
	rc = rcFromGridIndex(obj)[,c(2,1)]
	m[rc] = TRUE # the pattern
  
	## first get and sort the cell centers
	cv <- coordinatevalues(gr)
	cv[[1]] <- sort(cv[[1]])
	cv[[2]] <- sort(cv[[2]])

	## calculate cell boundaries
	x <- c(cv[[1]][1] - csiz[1]/2,  cv[[1]] + csiz[1]/2 )
	y <- c(cv[[2]][1] - csiz[2]/2,  cv[[2]] + csiz[2]/2 )
  
	# horizontal lines:
  	p = do.call(rbind, lapply(1:nbounds[2], function(i) {
			if (i == 1) # bottom line
				cells = m[1,]
			else if (i == nbounds[2]) # top line
				cells = m[nbounds[2]-1,]
			else # in-between, draw for:
				cells = m[i-1,] | m[i,]
			r = rle(cells) # figure out line pieces
			if (any(r$values)) {
				wr = which(r$values) # where to draw/end
				cs0 = c(0, cumsum(r$lengths)) # all start/end indices
				cbind(x[cs0[wr] + 1], y[i], x[cs0[wr+1] + 1], y[i])
			}
		}
	))
	segments(p[,1], p[,2], p[,3], p[,4], col = col, lty = lty, lwd = lwd)
	# vertical lines:
  	p = do.call(rbind, lapply(1:nbounds[1], function(i) {
			if (i == 1) # left boundary
				cells = m[,1]
			else if (i == nbounds[1]) # right boundary
				cells = m[,nbounds[1]-1]
			else # non-boundary lines, draw for:
				cells = m[,i-1] | m[,i]
			r = rle(cells)
			if (any(r$values)) {
				wr = which(r$values)
				cs0 = c(0, cumsum(r$lengths))
				cbind(x[i], y[cs0[wr]+1], x[i], y[cs0[wr+1]+1])
			}
		}
	))
	segments(p[,1], p[,2], p[,3], p[,4], col = col, lty = lty, lwd = lwd)
}

Try the sp package in your browser

Any scripts or data that you put into this service are public.

sp documentation built on May 29, 2024, 5:33 a.m.