1

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.

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

1 Answer 1

Reset to default
2

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""")
arcpy.DeleteFeatures_management(in_features="lines2D")

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")
Improve this answer
3
  • +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, 2016 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, 2016 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, 2016 at 1:01

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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