How would you do this in Pro?

Question

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

ThomasColson · Answer

Based on arcpy - Extracting value of raster for specific points (with X and Y coord) and write to column in Attribute Table? - Ge… 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) ‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

ThomasColson · Answer

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

GISMallku · Answer

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