[R-sig-Geo] Obtaining 3d surface area using polygon layer and DEM grid
Barry Rowlingson
b.rowlingson at lancaster.ac.uk
Wed Apr 14 17:12:55 CEST 2010
On Tue, Apr 13, 2010 at 5:46 PM, Barry Rowlingson
<b.rowlingson at lancaster.ac.uk> wrote:
> I don't know of any R implementation of a surface area algorithm for
> R at the moment.
Well, I do now because I've just written it. It's an implementation
of the Jenness algorithm. He doesn't say how he copes with the edges,
so what I do is extend the grid by 1 pixel all round and fill the edge
with the nearest value from the grid.
It's written in C and does a 2000x2000 grid in a couple of seconds.
> b=matrix(runif(2000*2000),2000,2000)
> system.time(print(sacalc(b)))
[1] 4592546
user system elapsed
1.868 0.076 1.941
now currently it doesn't take into account the grid spacing - it
assumes the cells are 1 unit apart and that's the same scale as the
values. Simple enough to add scaling but fiddly to integrate it into
SpatialGrid objects. Matrices are so much simpler :)
Also, it's one C file and one R file, no documentation, and I've not
tested it hardly, and I can't build it for Windows.
Barry
--
blog: http://geospaced.blogspot.com/
web: http://www.maths.lancs.ac.uk/~rowlings
web: http://www.rowlingson.com/
twitter: http://twitter.com/geospacedman
pics: http://www.flickr.com/photos/spacedman
More information about the R-sig-Geo
mailing list