xander_bakker

Modeling Archaeological Layers

Blog Post created by xander_bakker on Feb 23, 2015

A few week ago my attention was drawn by the thread: https://community.esri.com/message/453147 

 

The basic idea was to determine the soil layer of an archaeological find. In this thread 3D Analyst tools were used to obtain the result, but many steps were required to get there. I suggested to share some data to see if it was possible to make this process a little less manual and Armando was so kind to share not only part of this classified data, but also shared some insights in the wonders of archaeology.

 

Since the information is classified, the thread was made private. However, I do want to share the code and process that we worked. Maybe parts of it will be helpful for you or maybe you have suggestions to enhance the process.

 

Input data

The input data consists of two tables. One table holds the finds, with XYZ information:

tbl_finds.png

... and the other table hold topographic measurements and information on the soil type (the upper part of a soil layer):

tbl_topo.png

 

Outline of the process

A short description of the steps involved is:

  • Explore the data and determine
    • the extent of the xy coords of the find and topographic  information
    • the number of soil layers (a list of the names of the levels: [OH9, OH10, OH11])
  • Create points XYZ of the finds
  • Create points (XYZ) of the topographic coords per soil layer
  • Interpolate the topograhic points for each soil layer and convert to raster (3D Analyst)
  • Perform a "Extract Multi Values to Points" (Spatial Analyst). This adds the values of the topographic rasters as attributes to the finds.
  • Perform a loop to evaluate the attributes to find at what layer the find is located

 

A challenge

There was a challenge:

  • The finds have an extent a little larger than the topographic information. The topographic information should be extended to match the extent of the finds

points.png

 

Extending the topographic information

To extent the topographic information, I decided to create additional points at the corners of the data (the corners of the union of the finds and topographic extents). These corner points were added a Z value based on a weighted average of the 3 nearest points. Although this is an interpolation method, the influence of this error on the result is minimal.

 

Topographic points with additional corner points and the interpolation result (TIN, converted to raster):

Interpolation2.png

 

 

The result

The result of the analysis is the finds table converted to a z-aware featureclass enriched with the soil layer it was found in:

LayersAndFinds.png

 

 

The python code

Below the code that was used:

 

#-------------------------------------------------------------------------------
# Name:        archeo3D.py
# Purpose:     analyze finds in layers
#
# Author:      Xander
#
# Created:     10-02-2015
#-------------------------------------------------------------------------------

import arcpy
import os
import math

def main():

    # input tables
    tbl_finds = r"C:\Forum\Archeology\gdb\data.gdb\tbl_Finds"
    tbl_topo = r"C:\Forum\Archeology\gdb\data.gdb\tbl_Topo"
    fc_finds_name = "Find02" # output name
    fld_out = "FindLayerOH"

    # spatial reference and overwrite setting (True)
    sr = arcpy.SpatialReference(23030) # ED_1950_UTM_Zone_30N
    arcpy.env.overwriteOutput = True

    # settings
    cellsize = "CELLSIZE 0,01"

    # workspace
    ws = r"C:\Forum\Archeology\gdb\work.gdb"
    ws_tin = r"C:\Forum\Archeology\TIN"

    # output finds featureclass
    fc_finds = os.path.join(ws, fc_finds_name)

    # relevant fields for finds
    # fld_find_layer = "Occupation"
    fld_find_id = "Inventory"
    fld_find_x = "X_UTM"
    fld_find_y = "Y_UTM"
    fld_find_z = "Zr_MASL"

    # relevant fields for topo
    fld_topo_id = "PointNumber"
    fld_topo_layer = "OccupationTop"
    fld_topo_x = "X_UTM"
    fld_topo_y = "Y_UTM"
    fld_topo_z = "Z_MASL"

    # Check out extensions
    arcpy.CheckOutExtension("Spatial")
    arcpy.CheckOutExtension("3D")

    # list with fields
    flds_finds = [fld_find_id, fld_find_x, fld_find_y, fld_find_z]
    flds_topo = [fld_topo_id, fld_topo_layer, fld_topo_x, fld_topo_y, fld_topo_z]

    # create dictionary of finds
    dct_finds = createDictionaryFromTable(tbl_finds, flds_finds)

    # create dictionary of topo
    dct_topo = createDictionaryFromTable(tbl_topo, flds_topo)

    # get extent finds and topo and union
    ext_finds = getExtentFromDct(dct_finds, flds_finds, fld_find_x, fld_find_y)
    ext_topo = getExtentFromDct(dct_topo, flds_topo, fld_topo_x, fld_topo_y)
    ext_tot = unionExtents(ext_finds, ext_topo)

    # get soil layer list
    lst_oh = getUniqueListValues(dct_topo, flds_topo.index(fld_topo_layer))

    # create xyz featureclass for each layer with topo data
    lst_xyz = []
    for oh in lst_oh:
        fc_xyz = os.path.join(ws, "xyz_{0}".format(oh))
        lst_xyz.append(fc_xyz)
        createXYZfeatureclass(fc_xyz, sr)
        addFields2XYZfeatureclass(fc_xyz, flds_topo, tbl_topo)
        fillXYZfeatureclassOH(fc_xyz, dct_topo, oh, fld_topo_layer, flds_topo,
                            tbl_topo, fld_topo_x, fld_topo_y, fld_topo_z)
        updateXYZwithExtentCornerPoints(fc_xyz, ext_tot)

    # create the TIN's and convert to raster
    lst_multiras = []
    for oh in lst_oh:
        fc_xyz = os.path.join(ws, "xyz_{0}".format(oh))
        tin_xyz = os.path.join(ws_tin, "tin{0}".format(oh))
        in_params = "{0} Shape.Z Mass_Points <None>".format(fc_xyz)
        print in_params
        arcpy.CreateTin_3d(out_tin=tin_xyz,spatial_reference=sr,
                           in_features=in_params,
                           constrained_delaunay="DELAUNAY")

        # convert TIN to raster
        ras_xyz = os.path.join(ws, "ras{0}".format(oh))
        arcpy.TinRaster_3d(in_tin=tin_xyz, out_raster=ras_xyz, data_type="FLOAT",
                           method="LINEAR", sample_distance=cellsize, z_factor="1")

        lst_multiras.append([ras_xyz, "ras{0}".format(oh)])

    # create xyz featureclass finds
    createXYZfeatureclass(fc_finds, sr)
    addFields2XYZfeatureclass(fc_finds, flds_finds, tbl_finds)
    fillFindsfeatureclass(fc_finds, dct_finds, flds_finds, tbl_finds,
                            fld_find_x, fld_find_y, fld_find_z)

    # extract multi values to points
    arcpy.gp.ExtractMultiValuesToPoints_sa(fc_finds, lst_multiras, "NONE")

    # analyze attributes of finds featureclass
    print "analyze attributes of finds featureclass"
    flds_ohs = []
    for lst in lst_multiras:
        flds_ohs.append(lst[1])

    # add output field
    arcpy.AddField_management(fc_finds, fld_out, "TEXT", field_length=50)

    # update the finds featureclass with the soil layer the find belongs to
    flds = [fld_out, fld_find_z]
    flds.extend(flds_ohs)
    cnt = 0
    with arcpy.da.UpdateCursor(fc_finds, flds) as curs:
        for row in curs:
            cnt += 1
            lst_vals = list(row[2:])
            lst_srt = sorted(lst_vals, reverse=True)
            if row[1] > lst_srt[0]:
                layer = "above upper layer"
                layfound = None
            else:
                layfound = lst_srt[0]
                for laytop in lst_srt:
                    if row[1] < laytop:
                        layfound = laytop

            if not layfound == None:
                i = lst_vals.index(layfound)
                layer = flds[i+2]

            row[0] = layer
            curs.updateRow(row)


    # Check in extensions
    arcpy.CheckOutExtension("Spatial")
    arcpy.CheckOutExtension("3D")


def updateXYZwithExtentCornerPoints(fc, ext):
    """update XYZ featureclass with corner points of extent
       use IDW to assign Z value to the corner points"""
    flds = ('OID@', 'SHAPE@')
    dct_xyz = {r[0]: r[1] for r in arcpy.da.SearchCursor(fc, flds)}
    sr = arcpy.Describe(fc).spatialReference

    # insert cursor
    with arcpy.da.InsertCursor(fc, ('SHAPE@')) as curs:
        # loop through corners
        props = ["lowerLeft", "upperLeft", "upperRight", "lowerRight"]
        for crnr in props:
            dct_dist = {}
            pnt = eval("ext.{0}".format(crnr))
            for oid, pntg in dct_xyz.items():
                dist = getDistance(pnt, pntg)
                dct_dist[oid] = dist

            # sort the dict on value
            cnt = 0
            lst_dist = []
            lst_z = []
            for oid, dist in sorted(dct_dist.items(), key=lambda x: x[1]):
                cnt += 1
                lst_dist.append(dist)
                lst_z.append(dct_xyz[oid].firstPoint.Z)
                if cnt > 3:
                    break

            # distance weighted height for corners
            sum_dist = sum(lst_dist)
            z_res = 0
            for i in range(len(lst_dist)):
                d = lst_dist[i]
                z = lst_z[i]
                z_res += (d / sum_dist) * z

            # create the point and insert it
            pntz = arcpy.Point(pnt.X, pnt.Y, z_res)
            curs.insertRow((pntz, ))


def getDistance(pnt, pntg):
    """Calculate the 2D distance between Point and PointGeometry"""
    return math.hypot(pnt.X - pntg.firstPoint.X, pnt.Y - pntg.firstPoint.Y)


def fillFindsfeatureclass(fc_finds, dct_finds, flds_finds, tbl_finds,
                          fld_find_x, fld_find_y, fld_find_z):
    """Fill the finds featureclass based on finds table, use Z aware points"""
    flds = ["SHAPE@"]
    flds.extend(flds_finds)
    with arcpy.da.InsertCursor(fc_finds, flds) as curs:
        for oid, lst_values in dct_finds.items():
            pnt = arcpy.Point(lst_values[flds_finds.index(fld_find_x)],
                              lst_values[flds_finds.index(fld_find_y)],
                              lst_values[flds_finds.index(fld_find_z)])
            lst_row = [pnt]
            lst_row.extend(lst_values)
            curs.insertRow(tuple(lst_row))


def fillXYZfeatureclassOH(fc_xyz, dct_topo, oh, fld_topo_layer, flds_topo,
                        tbl_topo, fld_topo_x, fld_topo_y, fld_topo_z):
    """Fill the topo XYZ featureclass based on topo table
       filter information on provided OH (soil layer)"""
    flds = ["SHAPE@"]
    flds.extend(flds_topo)
    with arcpy.da.InsertCursor(fc_xyz, flds) as curs:
        for oid, lst_values in dct_topo.items():
            if lst_values[flds_topo.index(fld_topo_layer)] == oh:
                pnt = arcpy.Point(lst_values[flds_topo.index(fld_topo_x)],
                                  lst_values[flds_topo.index(fld_topo_y)],
                                  lst_values[flds_topo.index(fld_topo_z)])
                lst_row = [pnt]
                lst_row.extend(lst_values)
                curs.insertRow(tuple(lst_row))


def addFields2XYZfeatureclass(fc, flds, tbl):
    """Add field to FC based on list and extract settings for field from table"""
    for fld_name in flds:
        fld = arcpy.ListFields(tbl, wild_card=fld_name)[0]
        arcpy.AddField_management(fc, fld.name, fld.type, fld.precision,
                                  fld.scale, fld.length)

def createXYZfeatureclass(fc, sr):
    """Create an empty Z-aware point featureclass"""
    fc_ws, fc_name = os.path.split(fc)
    arcpy.CreateFeatureclass_management(fc_ws, fc_name, "POINT", has_z="ENABLED",
                                        spatial_reference=sr)


def createDictionaryFromTable(tbl, flds_att):
    """Create a dictionary with key is OID and values is a list of attributes"""
    flds = ["OID@"]
    flds.extend(flds_att)
    dct = {}
    with arcpy.da.SearchCursor(tbl, flds) as curs:
        for row in curs:
            oid = row[0]
            lst_vals = []
            for i in range(1, len(flds)):
                lst_vals.append(row[i])
            dct[oid] = lst_vals
    return dct


def unionExtents(ext1, ext2):
    """Union 2 extents"""
    return arcpy.Extent(min(list([ext1.XMin, ext2.XMin])),
                        min(list([ext1.YMin, ext2.YMin])),
                        max(list([ext1.XMax, ext2.XMax])),
                        max(list([ext1.YMax, ext2.YMax])))


def getExtentFromDct(dct, flds, fld_x, fld_y):
    """Determine the XY extent based on a dictionary with values"""
    lst_x = getListValues(dct, flds.index(fld_x))
    lst_y = getListValues(dct, flds.index(fld_y))
    return arcpy.Extent(min(lst_x), min(lst_y),max(lst_x), max(lst_y))


def getListValues(dct, i):
    """extract a list of values from a dictionary
       the value is a list of values and i is the index to the attribute"""
    return [lst_val[i] for lst_val in dct.values()]


def getUniqueListValues(dct, i):
    """Create a list of unique values from a dictionary
       the value is a list of values and i is the index to the attribute"""
    lst = [lst_val[i] for lst_val in dct.values()]
    return list(set(lst))

if __name__ == '__main__':
    main()

Outcomes