How would you do this in Pro?

[ThomasColson]

I have a new requirement that, in some 200+ ArcGIS Feature Services on Server 10.6.1 (They are not hosted via a PTL), they all get an "ELEVATION" field and users will somehow magically populate the elevation using a local, Lidar-derived DEM (so it can't be the AGOL elevation service). 

Here's what's not going to work:

Use Extract Values to Points (which creates another FC), join back to the original feature service, calculate field using RASTERVALUE: To many steps, users mostly stop understanding the workflow after the word "join". GIS is not these folks job, they GPS some things, add the points to a feature service, need to know the elevation...

The requirement is the Lidar-derived elevation, not the GPS one

Z-enabling everything. Not possible in our environment, nor feasible given the number of SDE feature classes we'd have to modify. 

Can't use SQL on the back end like 

ELEVATION = (SELECT pdata.getValueByLoc(1,p.SHAPE.Long,p.SHAPE.Lat)  FROM [dbo].[DEM10MP]),‍

which uses ST_RASTER, an unsafe assembly which we're not allowed to use. 

Surely someone (many hundreds in fact) have had this requirement....and solved it? Some python toolbox is out there somewhere......

[KoryKramer]

I haven't checked but would Add Surface Information—Help | ArcGIS Desktop work with the services?

[ThomasColson]

yes and no, the problem is I then have to script, or otherwise automate, the calculation of Z for each new point added. 

[ThomasColson]

Based on arcpy - Extracting value of raster for specific points (with X and Y coord) and write to column in A...  I baked the code sample into a toolbox which seems to be working well based on a morning of comparing tool results to spot checks: 

import arcpy
import numpy as np
update_feature_class = arcpy.GetParameterAsText(0)
source_dem= arcpy.GetParameterAsText(1)
elev_attrib= arcpy.GetParameterAsText(2)

rast = arcpy.Raster(source_dem)
desc = arcpy.Describe(rast)

ulx = desc.Extent.XMin
uly = desc.Extent.YMax
lly = desc.Extent.YMin
#Spatialrefernce
sr = desc.spatialReference


rstArray  = arcpy.RasterToNumPyArray(rast)

with arcpy.da.UpdateCursor(update_feature_class,["SHAPE@",elev_attrib]) as uc:
    for row in uc:
        pnt = row[0].projectAs(sr)
        #assuming every point falls to the left and below uperleft corner
        deltaX = pnt.centroid.X - ulx
        deltaY = lly- pnt.centroid.Y
        arow = int(deltaY/rast.meanCellHeight)
        acol = int(deltaX/rast.meanCellWidth)
        row[1] = rstArray[arow,acol]
        uc.updateRow(row)
‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍