Calculating volume of archaeological excavation levels from points
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.
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!)
(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.)
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.
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.)
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.
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:
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!
