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.

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.

Never miss an update!
Subscribe to R-bloggers to receive
e-mails with the latest R posts.
(You will not see this message again.)

Click here to close (This popup will not appear again)