(This article was first published on

**Bart Rogiers - Sreigor**, and kindly contributed to R-bloggers)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 his blog:**Bart Rogiers - Sreigor**.R-bloggers.com offers

**daily e-mail updates**about R news and tutorials on topics such as: visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...