Getting elevation for x,y when having multiple rasters

[BKuiper]

Hi,

I made a script to get the Z value for a list of terrains files at a specific X,Y, but the script can not be converted into a GPK because of the methods i'm using.

Does anybody have another good simple approach to do this, which can be packaged into a GPK ?

I will need to replace MosaicToNewRaster_management and GetCellValue_management, but i'm not sure how! I was thinking on using the Raster class and/or some NumPy approach, but can't figure out how to do it. Anybody have some leads on how to achieve this ?

import arcpy
from arcpy import env
from arcpy.sa import *

# Check out the ArcGIS Spatial Analyst extension license
arcpy.CheckOutExtension("Spatial")

arcpy.gp.overwriteOutput = 1

print "Start"

# original terrain data
terrainOne = r"C:\MyFolder\TerrainFiles\ned_25581390.flt";
terrainTwo = r"C:\MyFolder\TerrainFiles\ned_57214367.flt";

tempRasterFolder = r"C:\Myfolder\scratch.gdb"
tempRaster = "tempRaster"

arcpy.MosaicToNewRaster_management(terrainOne + ";" + terrainTwo, tempRasterFolder, tempRaster, number_of_bands=1, mosaic_method="LAST")

# location to get Elevation from
source = arcpy.Point( )
source.X = -87.6815094515746
source.Y =  41.811765962018300

# define string
sourceString = str(source.X) +  str(" ") + str(source.Y)
sourceElevationString = arcpy.GetCellValue_management(tempRasterFolder +"\\"+ tempRaster, sourceString)
#sourceElevationString = arcpy.GetCellValue_management(terrainOne, sourceString)

# get elevation
result = sourceElevationString.getOutput(0)

if result != 'NoData':
    source.Z = float(result)

    # print (converted meters to feet)
    print str(source.Z)
    print str(source.Z / 0.3048)
else:
    print "value (x,y) not found on raster"

[KevinHibma]

If you're still working on this, heres some code which you can use to get the cell value out of a raster.
I'm not really sure about the mosaic conversion part. Do you need to do this? You have rasters originally right? If thats the case, you should be able to run this code against them..

import arcpy

inPt = arcpy.GetParameterAsText(0)

ds = arcpy.Describe(inPt)

p = arcpy.Point(ds.extent.XMin, ds.extent.YMin)
cellValue = arcpy.RasterToNumPyArray("raster.tif", p, 1,1)

arcpy.AddMessage(cellValue)

This case takes a FeatureSet point from a user (expects a single point), gets the extent from that point, creates a Point object off it, then uses that as the start point to the RasterToNumPyArray function adding just 1 cell value in size from the extent.

Hopefully this moves you forward.

[BKuiper]

thanks Kevin. I'm still trying to do this and will try out your sample later this week.

Thank you!

[BKuiper]

Hi Kevin,

this works. Thanks! How would this work for a folder with .FLT files or a Mosaic Dataset within a GDB?