Friday, 7 December 2012

R and the SGeMS blockdata format 7.12.12

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))
} 
5 Bart Rogiers: R and the SGeMS blockdata format The popular geostatistical software SGeMS  has some options for working with non-point support ...

No comments:

Post a comment:

< >