I have a raster layer with 2 already determined classes: Low vegetation and High vegetation.

They don't mix, there is a characteristic border between them.

It's a typical situation of altitudinal zonation in the mountains.

I would like to determine height a.s.l. of this border based on a DEM that I have as a different layer. It would be easy just clicking along the line and writing down the DEM pixel value, but the whole point is that I need an automatic process to determine a high amount of borders. Is it possible to do?

I work with ArcGIS mostly.

The result I would like to get is any kind (vector/raster) of a line with height statistics I could interpret, so just a plain vector won't work here.

  • 1
    What is your ArcGIS licence and do you have a licence for spatial analyst extension? – fatih_dur Jun 12 '16 at 7:50
  • Hi, my license is "Advanced, single use" and yes, I have license for this extension. – adamczi Jun 12 '16 at 11:42

Solution assumes your classes are topologically correct polygons, i.e. there are no gaps and overlaps between them. Class defined by value stored in field GRIDCODE

#define boundary between polygons and delete outer ones
arcpy.PolygonToLine_management(in_features="polygons", out_feature_class="D:/Scratch/lines2D.shp", neighbor_option="IDENTIFY_NEIGHBORS")
arcpy.SelectLayerByAttribute_management(in_layer_or_view="lines2D", selection_type="NEW_SELECTION", where_clause=""""LEFT_FID" =-1""")

Create LEFT_CLASS and RIGHT_CLASS fields in lines table and populate them by relevant classes using joins to LEFT_FID and RIGHT_FID and polygons OID:

#transfer class on the left of the line to relevant line field
arcpy.AddJoin_management(in_layer_or_view="lines2D", in_field="LEFT_FID", join_table="polygons", join_field="FID", join_type="KEEP_ALL")
arcpy.CalculateField_management(in_table="lines2D", field="lines2D.LEFT_CODE", expression="[polygons.GRIDCODE]")
arcpy.RemoveJoin_management(in_layer_or_view="lines2D", join_name="")

#transfer class on the RIGHT of the line to relevant line field - repeat 3 processes above

# select lines that are boundaries between between 2 different classes
arcpy.SelectLayerByAttribute_management(in_layer_or_view="lines2D", selection_type="NEW_SELECTION", where_clause=""""LEFT_CODE" <> "RIGHT_CODE""""")

# find mean elevation
arcpy.InterpolateShape_3d(in_surface="dem", in_feature_class="lines2D", out_feature_class="D:/Scratch/lines_3D.shp", sample_distance="", z_factor="1", method="BILINEAR", vertices_only="DENSIFY", pyramid_level_resolution="0")
  • +1 because it will do the job. I know the OP said they need an automatic process, but it's not clear they want python code. – Adam Jun 12 '16 at 22:49
  • python is fine for me. as I wrote I have 2 rasters (not vector polygons), so first I would need to convert them, therefore these polygons should end up topologically correct. I'll try your code and report, thank you for this proposition. – adamczi Jun 13 '16 at 0:51
  • It is not really a Python, these are geoprocessing tools. Once you've got your 3d lines you might consider extracting their vertices to 3d points. It'll give you more freedom in terms of statistics, e.g. median elevation is not in add Z tool – FelixIP Jun 13 '16 at 1:01

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.