# R and the SGeMS blockdata format

[This article was first published on

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

**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.