Calculating volume of archaeological excavation levels from points

1944
6
12-11-2019 09:44 AM
TimDennehy
New Contributor II

Hi all,

I've been trying to perform what seems like it should be a fairly simple task, but I've been unable to get it to work using the simple "dummy" versions of the data I have, so I figured I'd ask here.

Essentially, I need to know the volume of sediment excavated in each level of several archaeological units. The units themselves are either square or rectangular pits, but the tricky part is that each level has irregular starting and ending depths. I have XYZ coordinates for the four corners and center point of each level. While the XY coordinates of those points remain essentially unchanged throughout the unit, the Z coordinates are different even for the top or bottom of the same level (because unfortunately we're less-than-perfect excavators, as it turns out). I'm trying to quickly turn these points into 3D representations in ArcGIS so I can calculate their volumes. I've found several posts about this but can't seem to get the procedures that were suggested to work. For example, I tried two workflows from this post: 

https://community.esri.com/message/369720?commentID=369720#comment-369720 

... but can't seem to get a simple "dummy" level (i.e. a cube 100 m on a side) to look right in ArcGlobe or to have the correct volume. Here's a stripped down version of workflow 1: 

Save coordinates in XYZ format as a csv

ASCII 3D to Feature Class (3D Analyst)

Create TIN (3D Analyst)

Surface Volume (3D Analyst)

A different version of a similar workflow suggested in that post is below:

Add Data (import csv file to map)

Display XY Data

Export Data as Feature Class

Create TIN from that Feature Class (3D Analyst)

Surface Volume (3D Analyst)

In both of the above workflows, the point features look correct in ArcGlobe (i.e. I get 8 points floating in space at the corners of a cube), but the TIN's produced always look non-cube-like, and the Surface Volume tool doesn't yield the correct calculations when applied to them. 

I also found a potential solution involving 3D polygons in another post:

https://community.esri.com/message/452823?commentID=452823#comment-452823 

I was able to create a cube-shaped 3D polygon, but couldn't figure out how to calculate volume measurements  because both Surface Volume and Polygon Volume tools require a TIN for the whole 3D shape.

Creating a TIN seems to be where I run into problems with every workflow I've tried. I'm not sure if this is because the point data aren't ordered correctly - I read somewhere that polygons needs to be "oriented clockwise", but I don't know what that means for a 3D shape. Or it may be because the xyz data isn't "closed" in the sense that the first and last vertices are different. I've tried making them the same in my dummy data by adding a ninth point that is the same as the first, but that yields an odd looking TIN with an edge stretching from diagonally from top to bottom. 

Any help would be greatly appreciated! Xander Bakker‌, I saw you had written a blog post back in 2015 dealing with some of the data from the 2nd post above, but couldn't figure out how to get volume calculations from it (and I'm not proficient in python - yet! - so wasn't able to follow some of what you write anyway). 

EDIT: I'm attaching both the dummy data (TextXYZz1) in .csv format and a sample of the actual data that I'm using in .xlsx format.

0 Kudos
6 Replies
XanderBakker
Esri Esteemed Contributor

Hi Tim Dennehy , 

Would it be possible to share a sample set of data? That would help me to understand your specific requirements and would speed up the process to find a solution. 

0 Kudos
TimDennehy
New Contributor II

Xander Bakker‌, I added the dummy and sample data to the original post, but in doing so, I somehow removed your reply. Thanks in advance for any help you can give!

XanderBakker
Esri Esteemed Contributor

Hi Tim Dennehy ,

I'm not sure if the data that you shared is representative for what you will work on, but is see some strange things:

In the data send in the Excel file I see some nodes that coincide. In the image I created some 3D triangle polygons based on the list of coordinates including the central point (see test code below). I exaggerated the X, Y and Z values to be able to see something:

def main():
    import arcpy
    import os

    arcpy.env.overwriteOutput = True

    # create spatial reference (WGS_1984_Web_Mercator_Auxiliary_Sphere)
    sr = arcpy.SpatialReference(3857)

    # list of x, y, z tuples with coordinates
    lst_crds = [(0, 0, 700),
                (0, 1000, 770),
                (1000, 1000, 760),
                (1000, 0, 740),
                (500, 500, 730)]

##    lst_crds = [(0, 0, 720),
##                (0, 1000, 770),
##                (1000, 1000, 760),
##                (1000, 0, 780),
##                (500, 500, 730)]

    # create a point for each x, y, z tuple in the list and append it to the list of points
    dct_i = {0: [0, 1, 4, 0],
             1: [1, 2, 4, 1],
             2: [2, 3, 4, 2],
             3: [3, 0, 4, 3]}
    lst_pols = []
    for i in range(4):
        lst_pnts = []
        lst_j = dct_i[i]
        for j in lst_j:
            crds = lst_crds[j]
            pnt = arcpy.Point(crds[0], crds[1], crds[2])
            lst_pnts.append(pnt)

        # create polygon with Z, but without M
        polygon = arcpy.Polygon(arcpy.Array(lst_pnts), sr, True, False)
        lst_pols.append(polygon)

    # create a featureclass with the polygon Z geometry
    ws = r'C:\GeoNet\3DVolume\data.gdb'
    fc_name = 'polygons_b02'
    fc = os.path.join(ws, fc_name)
    arcpy.CopyFeatures_management(lst_pols, fc)

if __name__ == '__main__':
    main()‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

After this, I created TINs and converted them to raster and did a CutFill analysis to determine the volume.

0 Kudos
XanderBakker
Esri Esteemed Contributor

Hi tdennehy , 

I also saw this publication on GeoNet this week, which might have relevant content: https://community.esri.com/community/learn-arcgis/blog/2019/12/12/underground-exploration-and-more-w... 

0 Kudos
TimDennehy
New Contributor II

Thanks very much Xander! Yeah, the data are a little strange... some corners (and sometimes the center) of the excavation unit were left unexcavated in a certain level while other corners were excavated. So co-inciding data points do occur sometimes. 

Just last night, I found the Cut Fill tool and figured out how to go from points to TINS to rasters to volumes. I noticed that the volumes calculated by Cut Fill were different from those calculated by Surface Volume or Polygon Volume using the TINS. Not sure why that would be the case, but I think Cut Fill is more straightforward so I'm going to go with that. 

One additional question - what did you use to create the image you posted, is that ArcScene? I'm having a lot of trouble getting my points and TINS to show up in either ArcGlobe or ArcScene at a zoom level where I can actually examine them. I've tried playing around with the visibility in the properties of each, but nothing seems to help - ArcGlobe only shows them when zoomed way out, and they don't even appear in ArcScene at all. Is there anything I'm missing to make them appear only at far-out zoom levels?

THANKS again for your help!

-Tim

0 Kudos
PeterKnoop
MVP Regular Contributor

Another 3D Analyst tool you might consider for tasks like this is Minimum Bounding Volume.

I find Minimum Bounding Volume to be very convenient for converting points to volumes in the context of archaeological excavations. I use it to turn the points that define the often uneven boundaries of levels (or contexts, or stratagraphic units, or whatever your preferred terminology is...) into closed multipatch features (i.e., 3D volume representations.) I use ArcGIS Pro, though if you are a Desktop user, it also has a version of this tool, Minimum Bounding Volume - Desktop.

As you've already done, the first step is loading the data from your Excel spreadsheet into a point feature class. You did not provide the coordinate system for your data; however, since it looks like an arbitrary example with units of meters, picking a coordinate system, like UTM, should do for this example.

I arbitrarily choose to work with your data in my local UTM Zone. Otherwise, by default, ArcGIS may assume a coordinate system that treats your X and Y values as degrees longitude and latitude, and your Z value as meters. As you can probably imagine, it would then take a lot of vertical exaggeration for a few centimeters of elevation difference to show up in a 3D view of your Level, when it is being treated as one-degree in areal size. Not to mention that you might get some very unexpected volume estimates as well!)

XY Table to Point

(I notice that your data seems to go up in elevation from your Start of the Level to the End of the Level, for point pairs with different elevations, rather than down, as I would expect for excavation. Your data will still work, however, it places your Start points below your End points, when you visualize it.)

Sample Level 1 Points

The second step is to generate multipatch features from each Levels' points. In your sample data you only have a single Level. Assuming your full data set contains more than one Level, then you would want to use the Level column to indicate which points should be grouped together to form an individual Level's volume. 

Minimum Bounding Volume

Because you have duplicate points for the same Level, you will see a warning message about such points being ignored in creating the multipatch feature. (When you have more than one Level in your data, then you should have points with duplicate coordinates, however, they will be associated with different Levels, so that won't interfere with each other; one Level's ending point is another's starting point, until you hit the bottom.)

Warning Message

when running the tool, checking the box for "Add geometry..." will add an attribute, MBV_Volume, which will contain the volume of your Level. This is not an auto-calculating attribute, however, so if you edit your multipatch feature, then you will want to re-calculate this field (e.g., Update Feature Z.)

Why would you edit your multipatch features? If you surveyed your "duplicate" points multiple times -- maybe it took more than a day to excavate a level, then the inaccuracies of the survey method will likely produce slightly different coordinates each time. So you might edit the multi-patch feature's vertices to average the difference, or to ensure one Level's vertices snap to an adjacent Level's. 

Also, because of those duplicate points, you will get an area that may not agree with how you are thinking of the area of this feature. Keep in mind that the whole "square" area is defined by the points you provided, however, more than half of the area the tool identified also has zero-thickness, and is, therefore, not part of the actual multipatch feature.

Sample Level 1 Volume

If you visualize the resulting points and multipatch features in a Local Scene in Pro, then you can see the 3D shape of your Level. Something like this:

Local Scene of Level 1 Sample Data 

In practice, levels/contexts/stratagraphic units can end up with some fairly convoluted shapes, so I would strongly recommend taking the time to visualize and explore the resulting multipatch features, before trusting their volumes (or areas). The "Concave Hull" method usually does the right thing, however, occasionally you may find that you need to manually drop some points that are duplicates or superfluous, which are confusing the algorithm relative to your expectations for the shape of the volume.

Hope that helps!

0 Kudos