Extracting highest point within a polygon using Terrain 5 dataset

[THomasDanks]

Hello,

This maybe be a stupid question/simple answer scenario. I have a Terrain 5 dataset (OS Terrain 5 – height DTM and contour data ), existentially an ASCII grid. Its loaded in my ArcMap as point, line and area layer. I am looking to extract the highest point within a polygon. For example I create a polygon shapefile around a section of town and then want to find the highest point  (location X,Y) within this polygon I've drawn. Is this possible to do this just using this shapefile and Terrain 5 dataset? Can I use a standard toolbox tool for this, or does it require a few more steps?

Any help would be greatly appreciated.

Regards,

[XanderBakker]

You can use the Zonal Statistics as Table tool (Zonal Statistics as Table—Help | ArcGIS Desktop), use the polygon as zone and the raster as raster value data to get statistics from a raster at that polygon (like the maximum value). In case you have multiple rasters, you might want to mosaic them together (Mosaic To New Raster—Help | ArcGIS Desktop )

[THomasDanks]

Thanks for your reply, I'll try this out and get back to you. Think I will have to mosaic a lot of rasters, as its an entire Country (Scotland). We have a corporate .mxd we use and this links to the Terrain 5 dataset for whole of Scotland which is stitched together. I can't get access to this file directly however. Is there anyway I can reference the Terrain 5 layer in my Arcmap view within the .mxd and run this tool? The "input raster value" box seems to want to select a file from your hard drive. Apologies if I'm missing the obvious here.

[XanderBakker]

Zonal Statistics as Table allows selection of a raster that you have in your TOC:

If it doesn't show up in the list or the option is not provided, it is not a valid format. 

Can you double click on the layer and check the datasource:

[THomasDanks]

Thanks again. It looks like the Terrain 5 layer I'm trying to use is not compatible with this tool as I am not getting that pulldown arrow next to "Input value raster". (I should of said I'm on Arcmap 10.2.1). I've had a look Terrain 5 layer source and the layer is in SDE Terrain format. I haven't got much experience with this format, will it be possible to use this dataset to achieve the above? Will I need to use another tool / combination of tools to accomplish it?

Thanks again

[XanderBakker]

You will first have to create a raster from it using the tool: Terrain To Raster—Help | ArcGIS Desktop  A Terrain dataset is not a raster format and the Zonal Statistics as Table requires a raster. 

If you just want to find the maximum value you could also adjust the stretching to the range of value (range closely around the maximum value) and simply visually identify it.

[THomasDanks]

Thanks, I'll look into doing this and post back if I fall foul of anyting

[THomasDanks]

Just coming back to say thanks for the help. All sorted now and have built a model to automate the process of obtaining the highest elevation point within polygons. Converts the results to a point layer containing elevation and coordinates.

Cheers again

[XanderBakker]

Hi THomas Danks ,

I have some additional questions: how often do you need to determine the maximum height and location within a given polygon from a terrain dataset? Is it always 1 polygon at the time? Are you the only end user or is this a model (or tool) you will share with other users? Is the polygon always available as featureclass or will the end user draw the polygon on the map?

Does the model convert the terrain to raster? This would probably mean long processing times, when it might be better to have a raster version pre-processed ready so you could skip that step.

I just created a terrain and tried something (manually) and came up with the following steps that seems to get the result:

• Set the polygon as extent
• Convert the terrain to point (will only create points for the extent)
• Select the points by the polygon (to remove those points outside the polygon)
• Add the Z value as attribute to the points
• Sort by Z value (descending) 
• Select the first record (highest Z value) and export to new featureclass.

Processing time is relatively short and the result is correct.

The whole process can be scripted, the last steps would require some lines of Python code, but nothing to worry about.

[XanderBakker]

To show you what this script could look like:

def main():
    import arcpy

    # inputs (could be made parameters to use script as tool)
    terrain = r'C:\GeoNet\TerrainDataset\test.gdb\FDS\Terrain'
    fc_tmp_pnts = r'C:\GeoNet\TerrainDataset\test.gdb\terrain2point_area4'
    fc_pol = r'C:\GeoNet\TerrainDataset\test.gdb\area4'
    fc_pnt = r'C:\GeoNet\TerrainDataset\test.gdb\highest_point4'

    # get first polygon from polygon featureclass
    polygon = arcpy.da.SearchCursor(fc_pol, ('SHAPE@')).next()[0]

    # set extent
    arcpy.env.extent = polygon.extent

    # convert terrain to points
    arcpy.CheckOutExtension("3D")
    arcpy.TerrainToPoints_3d(in_terrain=terrain,
                             out_feature_class=fc_tmp_pnts,
                             pyramid_level_resolution="0",
                             source_embedded_feature_class="<NONE>",
                             out_geometry_type="POINT")

    # get highest point (first one with highest value)
    z_max = -9999
    pntg_found = None
    cnt = 0
    with arcpy.da.SearchCursor(fc_tmp_pnts, 'SHAPE@') as curs:
        for row in curs:
            cnt += 1
            if cnt % 1000 == 0:
                print("Processing point: {0}".format(cnt))
            pntg = row[0]
            if pntg.firstPoint.Z > z_max:
                if polygon.contains(pntg):
                    pntg_found = pntg
                    z_max = pntg.firstPoint.Z

    # export output punt
    if pntg_found:
        arcpy.CopyFeatures_management([pntg_found], fc_pnt)
        print("Coords: X={0}, Y={1}, Z={2}".format(pntg_found.firstPoint.X,
                                                   pntg_found.firstPoint.Y,
                                                   pntg_found.firstPoint.Z))

if __name__ == '__main__':
    main()
‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍