An embedded and charset-unspecified text was scrubbed... Name: not available URL: <https://stat.ethz.ch/pipermail/r-sig-geo/attachments/20100413/4162feb4/attachment.pl>
Obtaining 3d surface area using polygon layer and DEM grid
8 messages · Sharon Baruch-Mordo, Edzer Pebesma, Barry Rowlingson
On Tue, Apr 13, 2010 at 3:16 PM, Sharon Baruch-Mordo
<sharonb_m at yahoo.com> wrote:
Hello list, Is there anyway to calculate in R?the 3D surface of a polygon spatial object using a DEM grid similar to the ArcView tool developed by Jenness (2004; http://www.jennessent.com/arcview/surface_tools.htm)?
Any way? Well yes, but that way may mean you implementing the algorithm described in the paper. It looks quite simple, and probably best written in C for speed. A pure R implementation might be very slow for large DEMs. I don't know of any R implementation of a surface area algorithm for R at the moment. Barry
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
Barry, I wouldn't object integrating it in sp. Can it cope with missing values in the grid? Would you be willing to write a short help file with some reference to the algorithm used? (I found this one: http://www.jennessent.com/downloads/WSB_32_3_Jenness.pdf )
Barry Rowlingson wrote:
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 :)
That's why SpatialGridDataFrames have an as.matrix method. Of course you loose the cell size.
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
Edzer Pebesma Institute for Geoinformatics (ifgi), University of M?nster Weseler Stra?e 253, 48151 M?nster, Germany. Phone: +49 251 8333081, Fax: +49 251 8339763 http://ifgi.uni-muenster.de http://www.52north.org/geostatistics e.pebesma at wwu.de
An embedded and charset-unspecified text was scrubbed... Name: not available URL: <https://stat.ethz.ch/pipermail/r-sig-geo/attachments/20100414/098644dd/attachment.pl>
On Wed, Apr 14, 2010 at 4:38 PM, Edzer Pebesma
<edzer.pebesma at uni-muenster.de> wrote:
Barry, I wouldn't object integrating it in sp. Can it cope with missing values in the grid? Would you be willing to write a short help file with some reference to the algorithm used? (I found this one: http://www.jennessent.com/downloads/WSB_32_3_Jenness.pdf )
Not sure what to do about missing values, but that would be the way to deal with non-rectangular areas. An NA cell would obviously contribute 0 to the total surface, but it also affects the calculation of surface areas for the eight adjacent cells. Currently it's just C code called with .C, so dealing with NA is inherently troublesome anyway! Maybe I can rewrite to use .Call, but there's a cricket match on... Barry
On Wed, Apr 14, 2010 at 5:00 PM, Sharon Baruch-Mordo
<sharonb_m at yahoo.com> wrote:
Thanks you Barry and Edzer for your replies! I should have mentioned that I'm a relatively?new to?spatial analysis in R,?and unfortunately I'm not familiar with programing in C. But I'm ready to learn! In trying to understand the code below, I see that?you are generating a 2000x2000 matrix with a random uniform number then?using the function sacalc, but how is it calculating the 3d surface (I tried the web for documentation on sacalc but didn't get too far... )??Would that be the equivalent of getting a DEM clipped by the polygon of interest and then running it?through sacalc?
Ah, sorry, sacalc() is the function I wrote, and it's not included. The line: b=matrix(runif(2000*2000),2000,2000) just creates a 2000x2000 matrix of random numbers (representing heights) between 0 and 1. It won't be a very smooth surface! The resulting output being about 2000x2000 shows I'm getting in the right ballpark. In fact if I change the values in "b" to be smaller, meaning a relatively flatter surface compared to the grid size, the result tends towards 2000x2000, which is comforting. sacalc just works on a 2d matrix of values, and treats them as heights in the DEM. It can't clip it to a polygon, but Edzer's comments on NA (missing data) values has given me some ideas... Barry
On Wed, Apr 14, 2010 at 5:03 PM, Barry Rowlingson
<b.rowlingson at lancaster.ac.uk> wrote:
?Not sure what to do about missing values, but that would be the way to deal with non-rectangular areas. An NA cell would obviously contribute 0 to the total surface, but it also affects the calculation of surface areas for the eight adjacent cells. ?Currently it's just C code called with .C, so dealing with NA is inherently troublesome anyway! Maybe I can rewrite to use .Call, but there's a cricket match on...
Right, well I didn't realise you could call ISNA(x) from .C code. I've just implemented something for handling NA... When considering a cell with an NA value, the area contribution is zero (whatever its neighbours are) When considering a non-NA cell with NA-valued neighbours, pretend the neighbour cell had the same value as the considered cell. I've implemented this, and it seems to be consistent with expectations... Will sort out docs tomorrow. Barry