R and the SGeMS blockdata format
[This article was first published on Bart Rogiers - Sreigor, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.
The popular geostatistical software SGeMS has some options for working with non-point support (block) data through the BGeost set of algorithms by Yongshe Liu (see his PhD thesis), and published in Liu and Journel (2009). A specific but fairly simple data format is required to get your block data into SGeMS.
If you have a set of square data supports (either for conditioning data or interpolation/simulation targets) you can use the R function below with the square center coordinates (centersX, centersY), the side lengths (lengthX, lengthY), the parameter value (values), block error estimate (error), and the output filename (outfile). Rectangles are possible too by using different side lengths, but the orientation should always be parallel to the x- and y-axis. Adding a rotation angle (or third dimension) would be straightforward but I didn’t need it, so it is not there yet ;).
################################################################################
# FUNCTION - SGeMS write square blockdata ######################################
################################################################################
write.square.blockdata <- function(centersX, centersY, lengthX, lengthY, values, error, outfile)
{
blocknumbers <- length(centersX)
blockvalues <- values
blockwest <- centersX - (lengthX/2)
blockeast <- centersX + (lengthX/2)
blocknorth <- centersY + (lengthY/2)
blocksouth <- centersY - (lengthY/2)
blockerror <- error
blockdata <- data.frame(nr=blocknumbers, value=blockvalues, east=blockeast, west=blockwest, north=blocknorth, south=blocksouth, error=blockerror)
# This function writes an SGeMS blockdata file with rectangular blocks
# The header
cat('Blockdata file created in R\n', file=outfile)
# Number of blocks
cat(length(blockdata[,1]), file=outfile, append=TRUE)
cat('\n', file=outfile, append=TRUE)
# Block loop
for(i in 1:length(blockdata[,1]))
{
write(paste('block #', i, sep=''), file=outfile, append=TRUE) # block name
write(blockdata$value[i], file=outfile, append=TRUE) # block value
write(blockdata$error[i], file=outfile, append=TRUE) # block error
write(paste(blockdata$west[i], blockdata$south[i], '0', sep=' '), file=outfile, append=TRUE) # lower left corner
write(paste(blockdata$west[i], blockdata$north[i], '0', sep=' '), file=outfile, append=TRUE) # upper left corner
write(paste(blockdata$east[i], blockdata$north[i], '0', sep=' '), file=outfile, append=TRUE) # upper right corner
write(paste(blockdata$east[i], blockdata$south[i], '0', sep=' '), file=outfile, append=TRUE) # lower right corner
}
}In case you have to deal with circular spatial supports, the function below might be useful. The arguments are comparable to the ones above, except there is a radius now, and the number of sides for a polygon approximation of the circular support. These functions certainly made my life more easy, I hope they will be useful to some of you too!
################################################################################
# FUNCTION - SGeMS write circular blockdata ####################################
################################################################################
sgems.write.circular.blockdata <- function(centersX, centersY, radii, values, error, npolygon=20, outfile)
{
### This function writes a 2D SGeMS blockdata file with default icosagon approximation of circular support blocks
## Header
cat('Blockdata file created in R\n', file=outfile)
## Number of blocks
cat(length(centersX), file=outfile, append=TRUE)
cat('\n', file=outfile, append=TRUE)
## Block loop
for(i in 1:length(centersX))
{
write(paste('block #', i, sep=''), file=outfile, append=TRUE) # block name
write(values[i], file=outfile, append=TRUE) # block value
write(error[i], file=outfile, append=TRUE) # block error
circlepoints <- pointsOnCircle(centersX[i], centersY[i], radii[i], npolygon)
circlepoints <- data.frame(circlepoints, z=circlepoints$x*0)
write.table(circlepoints,file=outfile, append=TRUE, sep='\t', col.names=FALSE, row.names=FALSE)
}
}
pointsOnCircle <- function(centerX, centerY, radius, n)
{
# This function creates a data frame with n points equally spaced on a circle
alpha = pi * 2 / n
pointsX = c(rep(NA,n))
pointsY = c(rep(NA,n))
i = -1
while( i < n-1 )
{
theta = alpha * i
pointOnCircleX = (cos( theta ) * radius) + centerX
pointOnCircleY = (sin( theta ) * radius) + centerY
pointsX[ i+2 ] = pointOnCircleX
pointsY[ i+2 ] = pointOnCircleY
i <- i+1
}
return(data.frame(x=pointsX, y=pointsY))
}
To leave a comment for the author, please follow the link and comment on their blog: Bart Rogiers - Sreigor.
R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.