Mathematica Stack Exchange is a question and answer site for users of Mathematica. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

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.

share|improve this question
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
up vote 12 down vote accepted

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"];

enter image description here

Something like:

triangleBuilder[lines_] := 
 Module[{triangles, firstSelection, secondSelection, aux},
  triangles = {};
   (*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*)
     DeleteDuplicates[Join[secondSelection[[j]], lines[[i]]]]],
    {j, Length[secondSelection]}];
   , {i, Length[lines]}];
  (*clean duplicate triangles*)
      Function[{a, b}, 
       If[a[[1]] != b[[1]], a[[1]] > b[[1]], a[[2]] > b[[2]]]]] &, 

And you obtain this:

enter image description here

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]]

enter image description here

You can ask Mathematica for different representation:

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

enter image description here

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

enter image description here

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.

enter image description here

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:

enter image description here

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.

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

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.