R/utilities.R

In aurelhy: Hydrometeorological interpolation

# Size of one degree in latitude/longitude, according to the WGS84 ellipsoid
"deg.lat" <- function (latitude) {
	# Latitude is in decimal degrees => phi in radians
	phi <- latitude / 180 * pi
	# See http://en.wikipedia.org/wiki/Latitude 
	111.132954 - 0.559822 * cos(2*phi) + 0.001175*cos(4*phi)
}

"deg.lon" <- function (latitude) {
	# Note that size of one degree in longitude depends only on the latitude!
	# Latitude is in decimal degrees => phi in radians
	phi <- latitude / 180 * pi
	# See http://en.wikipedia.org/wiki/Longitude#Length_of_a_degree_of_longitude
	cos(phi) * pi * 6378.137 / (180 * sqrt(1 - 0.00669437999014* sin(phi)^2))
}

# Given a point, calculate angle and distance from grid data (geomat object)
"polar.coords" <- function (geomat, x, y, maxdist)
{
	if (!inherits(geomat, "geomat"))
		stop("'geomat' must be a 'geomat' object")
	if (!missing(maxdist) && !is.null(maxdist)) {
		# maxdist is in km, but x and y are in decimal degrees
		# Calculation of the size of ane degree in latitude/longitude according to
		# central latitude in the considered geographical area
		meanlat <- mean(range(coords(geomat, type = "y")))
		lenx <- deg.lon(meanlat)
		maxdegx <- maxdist / lenx
		leny <- deg.lat(meanlat)
		maxdegy <- maxdist / leny
		# Filter out data contained in a square of maxdist * 1.05
		mx <- maxdegx * 1.05
		my <- maxdegy * 1.05
		xlim <- c(x - mx, x + mx)
		ylim <- c(y - my, y + my)
		# Take a window out of these data
		geomat <- window(geomat, xlim, ylim)
	}
	
	coords <- coords(geomat, "xy")
	# Recenter coordinates around x and y
	X <- coords$x - x
	Y <- coords$y - y
	# Note: distances are in km
	angles <- atan2(Y, X)	# Angles go from -pi to pi, and we want 0 to 2 * pi
	angles <- ifelse(angles > 0, angles, 2 * pi + angles)
	res <- data.frame(angle = angles, dist = sqrt((X*lenx)^2 + (Y*leny)^2))
	attr(res, "geomat") <- geomat
	return(res)
}

# match.coords matches numerical data with a tol value
# TODO: optimize this, considering we have a grid!
"match.coords" <- function (points, table, tol = 0.002)
{
	# Look if point is in the tol vicinity of one point in table
	match.one <- function (point, table, tol)
		any(point[1] - tol < table$x && point[1] + tol > table$x &&
			   point[2] - tol < table$y && point[2] + tol > table$y)
	res <- apply(points[, c("x", "y")], 1, match.one, table = table, tol = tol)
	return(res)
}

# A new coords method
"coords" <- function (x, ...)
	UseMethod("coords")

# New resample method
"resample" <- function (x, ...)
	UseMethod("resample")

# New add.points method
"add.points" <- function (x, ...)
	UseMethod("add.points")
	
# Augment a geomask with points near geopoints coordinates
"add.points.geomask" <- function (x, geopoints, ...)
{
	# Check arguments
	if (!inherits(x, "geomask"))
		stop("'x' must be a 'geomask' object")
	if (!inherits(geopoints, "geopoints"))
		stop("'geopoints' must be a 'geopoints' object")
	# Look for the points that are closest to the geomask grid in the 'geopoints' object
	# Quick calculation by transforming coordinates into their closest index
	# in the grid
	Coords <- coords(x)
	X = (geopoints$x - Coords["x"]) / Coords["size"]
	Y = (geopoints$y - Coords["y"]) / Coords["size"]
	# Round values for X and Y to match a grid point
	Xind = round(X) + 1
	Yind = round(Y) + 1
	# Eliminate points that are outside the grid
	In <- rep(TRUE, length(X))
	In[Xind < 1] <- FALSE
	In[Xind > nrow(x)] <- FALSE
	In[Yind < 1] <- FALSE
	In[Yind > ncol(x)] <- FALSE

	# Check which point is inside the grid, but not in the mask (avoid negative
	# indices... points eliminated anyway)
	toAdd <- !diag(x[abs(Xind), abs(Yind)]) & In
	# Add these points to the mask
	added <- data.frame(x = Xind[toAdd], y = Yind[toAdd])
	x[added$x, added$y] <- TRUE
	# ... and record the list of added points as "added" attribute
	added2 <- attr(x, "added")
	if (is.null(added2)) {
		attr(x, "added") <- added 
	} else {
		# Merge added points
		attr(x, "added") <- rbind(added2, added)
	}
	return(x)
}

# Distance to the sea from DEM
"dist2sea" <- function (geotm)
{
	# Edges: find sea pixels that have at least one terrestrian neighbour
	"externaledge" <- function (geotm, o) {
		edgematrix <- geotm
		edgematrix[edgematrix >= 0] <- 0
		for (i in 1:nrow(o)) {
		    if (o[i, 1] == 1) x1 <- o[i, 1] else x1 <- o[i, 1] - 1
		    if (o[i, 1] == nrow(geotm)) x2 <- o[i, 1] else x2 <- o[i, 1] + 1
		    if (o[i, 2] == 1) y1 <- o[i, 2] else y1 <- o[i, 2] - 1
		    if (o[i, 2] == ncol(geotm)) y2 <- o[i, 2] else y2 <- o[i, 2] + 1
		    if (min(is.na(geotm[x1:x2, y1:y2])) == 0)
				edgematrix[o[i, 1], o[i, 2]] <- 1
		}
		return(edgematrix)
	}

	dist2seaout <- geotm
    o <- which(is.na(dist2seaout), arr.ind = TRUE)
    edgematrix <- externaledge(dist2seaout, o)
    xyland <- which(edgematrix == 0,  arr.ind = TRUE)
    xyedge <- which(edgematrix == 1,  arr.ind = TRUE)
    cellsize <- coords(geotm)["size"]
    lat0 <- coords(geotm)["y1"]
    for (i in 1:nrow(xyland)) {
        dist2seaout[xyland[i, 1], xyland[i, 2]] <- min(sqrt(((xyland[i, 1] -
			xyedge[, 1]) * cellsize * 111.2 * cos((xyland[i,2] *
			cellsize + lat0) * pi / 180))^2 + ((xyland[i,2] - xyedge[,2]) *
			cellsize * 111.2)^2))
    }
    return(dist2seaout)
}