# Volumetric analysis from topological maps

I want to calculate the volume of fill required to level an area about the size of three foodball fields. The Department of Public Works in my city can provide me with topo data for the area. Their request form lists ArcView shapefile, AutoCAD.dwg, Geodatabase, and "other" data file formats, and JPEG, TIFF, and PDF map file formats. What kind(s) of files should I ask for? How do I calculate the volume of data required for an area I'm interested in? I don't know anything about this stuff.

I would also like to estimate the number of trees, associated carbon dioxide sequestration, etc. for the area, and other things related to the wildlife in the area (based on some assumptions relevant to my area and the type of vegetation in the area). This additional question may not really be a Mathematica question, but if there is something relevant in Mathematica's curated data, I would appreciate suggestions.

-
How big is a foodball field? And what's the ball made of? –  Jens Jan 27 '13 at 18:10
Foodball is related to Calvin Ball, and the ball can be made of anything I want. But if you are going to be so fussy, I'll stipulate that foodball is played on a football field. –  George Wolfe Jan 27 '13 at 18:26
You might want to ask these questions instead on the GIS site: that community is aware of software designed specifically to manage such data and perform these computations. If you choose to do this, please flag this question for moderator attention and ask them to migrate it. –  whuber Jan 27 '13 at 19:10

Mathematica is not really a reference in digital terrain models, and, there are very powerful software packages to deal with geographical information.

But... where would be the fun...

If the information in the files includes a digital terrain model (DTM), with the surface of the terrain defined by triangular faces (regular or not), then you can easily calculate the volume contained in the triangular projection down to a certain defined level.

Obviously, only CAD files have vector information (and the PDF format, but I wouldn't recomend it), and so, I believe that their option of JPG and TIFF is a no go, since it is probably just the Rasterize[CAD].

Since we are talking of such a small area, I believe there's no need for projection corrections (the fact that the earth is round kind of stuff).

To facilitate the algorithm, you first find the volume to a level below all levels of all the triangles, and then you add the box volume up to your desired level. This avoids mathematically dealing with triangles that cross the desired terrain level (that have a negative and positive filling).

Most of the trouble you will have is dealing with the information format.

If the triangles are defined just by individual lines, first you will need to join them together back into triangles:

lines = Import["dtm.dxf", "LineObjects"];


Something like:

triangleBuilder[lines_] :=
Module[{triangles, firstSelection, secondSelection, aux},
triangles = {};
Do[
(*first, for each line,
select the other ones that share the same vertice*)
firstSelection =
Select[lines, Length@Intersection[lines[[i]], #] == 1 &];
(*Then, from those,
we select the one that, for which, there's another line that links \
the vertices of it to the one of the originaly selected line*)
secondSelection =
Select[firstSelection, (aux =
Complement[Join[#, lines[[i]]], Intersection[#, lines[[i]]]];
MemberQ[firstSelection, aux] ||
MemberQ[firstSelection, Reverse@aux]) &];
(*clean duplicates, add to the triangle list*)
Do[
AppendTo[triangles,
DeleteDuplicates[Join[secondSelection[[j]], lines[[i]]]]],
{j, Length[secondSelection]}];
, {i, Length[lines]}];
(*clean duplicate triangles*)
DeleteDuplicates[
Map[Sort[#,
Function[{a, b},
If[a[[1]] != b[[1]], a[[1]] > b[[1]], a[[2]] > b[[2]]]]] &,
triangles]]
]


And you obtain this:

Once you have the list of triangles, you can easily calculate the volume contained below them. (from its name, this information is probably directly available from the ArcView shapefile, if you can access it, but I'm no expert...)

But if the information you have are just points, or isolines (line having the same height), that is, if you don't have the triangles, things get a little more complicated to do it right.

points = Import["map.dxf", "VertexData"];

ListPointPlot3D[points, ColorFunction -> "Topographic",
PlotStyle -> PointSize[0.002]]


You can ask Mathematica for different representation:

ListPlot3D[points, ColorFunction -> "Topographic",
MeshFunctions -> {#3 &}, Filling -> Bottom, FillingStyle -> Brown]


ListPlot3D[points, ColorFunction -> "Topographic", Mesh -> All,
Filling -> Bottom, FillingStyle -> Brown]


And then you could use Mathematica's triangles to do the volume calculation (extracting them from the Graphics element).

Or I imagine you could build a 3D Interpolation function from the points, and work directly with it (probably, NIntegrate it to your desired level).

But I don't know if you can then explain how the mesh was calculated (since it belongs to Mathematica internals).

So, I would recommend picking up a secific algorithm, like a reciprocal distance algorithm (see: Computational Geosciences with Mathematica By William C. Haneberg; if you are lucky, Google Books shows you page 290; I will not copy the code here, since it probably breaks some copyrights), to obtain a rectangular point map.

With a rectangular grid, you can easily calculate the volume by the triangles...

Long time ago, GIS also used raster information as data source, but I don't think this is what they are talking about when they say JPG, TIFF, etc.

You can see that, as long as you have the scale information, you can calculate the volume from each pixel color/tone:

In what respects other elements of the GIS, well, things can get even fuzzier without more specifics...

After this introduction, take also a look at: ArcGRID and GeoTIFF, and probably try to find examples (over the internet) of the formats they can supply, to see if you can work with them.

-
You can try Non-grid interpolation package, if you can't find Haneberg's book –  Tuku Jan 28 '13 at 21:23