topic Re: How to find Elevation difference along the river line? in ArcGIS Spatial Analyst Questions

How to find Elevation difference along the river line?

PawanThapa — Fri, 06 Oct 2017 11:12:07 GMT

Dear GIS experts,

I would be very thankful to receive your help/support/suggestion for solving my following problem. I'm using ArcGIS 10.1 with Spatial Analyst and 3D Analyst(no ArcInfo license) and Arcview 3x (still for many plusses).

I have a river network - both a shape file and raster grid - generated from a DEM through Flow Direction and Flow Accumulation grid for a very big region (a country).

What I want to know/calculate is - for every river cell, what is the difference in elevation value between the elevation of that river cell and the elevation of the river cell at a given distance (say, 1 km, 2km, 5km) downstream from that river cell along that river i.e. flow direction?

I tried to do this by using the focal statistic function running over the River-DEM (using minimum as statistic) and then subtracting the resulting grid from the original River-DEM; but this dos not work completely correctly because at many places the minimum in the neighbourhood is found in the adjacent river (tributary or another parallel river at another side of the ridge) and it gives me wrong elevation difference because the minimum river cell is not in the same river downstream but in another river although within the defined search/neighbhourhood distance. How can I get this correctly- i.e. the elevation difference in the same river downstream in a given distance?

I hope that I could make this my problem sufficiently clear. What would be your suggestion to solve this?

Thanking you a lot!!

Re: How to find Elevation difference along the river line?

DanPatterson_Retired — Fri, 06 Oct 2017 12:39:58 GMT

You would have to do this for each river channel, either moving upstream or down stream, branching where tributaries meet. (densify a channel polyline at a spacing of 1 cell width, convert to points, extract values to points)

Once the data set(s) are constructed, you need to perform a rolling difference, since you specify that you want the distance difference between two fixed locations, but for every location as you move upstream. (rolling anything isn't readily implemented in arcmap, but is easy in numpy and python, so you may wish to examine those options)

So basically what is you purpose in doing this since it isn't clear what the information is required for and how an elevation difference at 500m spacing at 10 m increments is going to be more useful that extracting the dem values along each stream reach and processing from there. What terrain derivatives do you need? slope? aspect? curvature? etc?

Any more insights on you final goal may spawn more suggestions.

Re: How to find Elevation difference along the river line?

XanderBakker — Fri, 06 Oct 2017 12:45:16 GMT

Interesting problem. I guess doing this in raster would be too complicated or inefficient. Maybe that Arc Hydro‌ has something to solve this particular problem, but if not I would probably go for a vector approach. For each inlet trace downstream to the end of the stream. Make sure the lines will be z-enabled. For each interval over each line (will create duplicate points further down the streams) get the elevation at that points and at the position 1, 2 or 5 km further down the stream. Read out the elevation of that point and return the values to the point. This will require some programming to get the job done. Hopefully there is something in ArcHydro that can help to get the job done.

Re: How to find Elevation difference along the river line?

PawanThapa — Fri, 06 Oct 2017 13:34:02 GMT

Sure Dan, let me give more insight on my goal to obtain exact solution so far as possible.

I'm doing this work in the field of "mapping the hydropower potential" in the river network of a country. Aim is to show the hydropower potential in each river cell of the river network (which is 30m X 30m from SRTM-DEM) , when the waterway (or powerhouse) is kept at, say 5km below/downstream from that river cell in the same river. So, the calculation should successively move down cell by cell, each time/cell calculating the elevation difference between that cell and the cell 5km downstream, through all river network.

As obvious the power potential at a river cell is function of the product of flow and elevation head at that cell (Head is the elevation difference between that considered cell and the cell in the same river 5km downstream). Through hydrological analysis and GIS, I have river raster with flow values at each river cell; and now remains to be calculated is a river raster where each river cell is showing the elevation head for that 5km river stretch downstream starting from that river cell. This is where I'm getting stucked.

I hope I could make my problem and aim clearer to you!

Re: How to find Elevation difference along the river line?

PawanThapa — Fri, 06 Oct 2017 13:38:19 GMT

Thanks Xander for the reply/ suggestion ..... So far I could not find solution through ArcHydro as well .. I was thinking somehow I could use flow direction map as well in combination to find appropriate solution method doing in raster but could not find a way out.

I think I have to learn Python for this (time-permitting for this quick job) -- I know only FORTRAN programming- quite old..

Re: How to find Elevation difference along the river line?

DanPatterson_Retired — Fri, 10 Dec 2021 20:13:35 GMT

A small illustration on calculating moving sequential differences. Assume that you have your river sampled at a finite increment.

a = np.array([1, 2, 4, 6, 7, 8, 8.5, 9, 9.2, 9.7, 9.9, 10])  # ---- river elevation ----

# samples

a[2:] - a[:-2]  # difference between every 2nd
array([ 3. ,  4. ,  3. ,  2. ,  1.5,  1. ,  0.7,  0.7,  0.7,  0.3])

a[3:] - a[:-3]  # difference between every 3rd
array([ 5. ,  5. ,  4. ,  2.5,  2. ,  1.2,  1.2,  0.9,  0.8])

# ---- etcetera ----‍‍‍‍‍‍‍‍‍‍‍

Just converted this to a numpy array since it is the easiest to deal with these types of questions.

Re: How to find Elevation difference along the river line?

PawanThapa — Mon, 09 Oct 2017 10:06:37 GMT

Dear Dan,

Thank you for the suggestion. But main problem with the array calculation as you have suggested is, as far as I have understood- it does not follow the flow direction from the calculating cell to the end cell (2nd, 3rd etc.. or 166th cell for 5km elevation difference --my raster is 30m x 30m..), which will become problem when there are several branching rivers/tributaries within that 2-5 km distance conflencing towards a main river-- -- or am I missing or misunderstanding something in your suggestion?

Re: How to find Elevation difference along the river line?

XanderBakker — Mon, 09 Oct 2017 11:28:42 GMT

Could you share a sample of your data (DEM and the flows as lines)? I wonder what can be done with some scripting. It won't be avenue, since I left that behind me a long time ago (just like AML and ArcObjects)

Re: How to find Elevation difference along the river line?

PawanThapa — Mon, 09 Oct 2017 12:55:20 GMT

Dear Xander,

Thanks a lot.

I have extracted a small area from my working-region as a sample data and attached herewith. But these are two raster grids-- one-RiverDEM showing Elevation value at each cell (30m X30m) and another-Flowdir showing at each cell its flow direction to next cell (standard D-8 algorithm of flow direction in GIS). At this location of sample data as well, I had that problem while trying to find the elevation difference (Head) for each river cell in a 5km river stretch from that each river cell.

I hope this gives you now further clearer info about my need/problem.

Thanks a lot!

Re: How to find Elevation difference along the river line?

PawanThapa — Mon, 09 Oct 2017 13:02:31 GMT

Forgot to mention, but as already said--how I was doing--- I use focal statistics to calculate minimum with 5km search radius over the river DEM and i subtract the original RivDEM with the result of the focal statistics. The problem here is, in that search radius of 5km, for several cells of RiverDEM it found the minimum at adjacent rivers instead of in the same river downstream, which in reality is at another side of the mountain ridge, hence giving me unwanted results. When the adjacent rivers are more than 5km away, then the results are ok but that will not be the case everywhere, mainly toward confluencing or branching of the rivers....

Re: How to find Elevation difference along the river line?

XanderBakker — Fri, 10 Dec 2021 20:13:38 GMT

Hi Pawan Thapa , thanks for sharing the data.

This is what I have come up with. Part of the process was done manually, but a large part of that can be easily automated.

Since the raster data did not have any spatial reference defined I simply assigned one, although it will probably not be the correct one.

I created 2 featureclasses; one with the outlet of the network (red point) and one with some of the inlets (green points):

Then I took your flowdir raster and created a new raster called cost that simply has value 1 where the rivers are defined and NoData for all the other pixels. Python snippet:

arcpy.gp.RasterCalculator_sa('Con("flowdir", 1)', "C:/GeoNet/DEMstreams/test.gdb/cost")‍

Next I used this raster in a Cost Distance calculation to generate a cost distance and backlink raster. Python snippet:

arcpy.gp.CostDistance_sa("outlet", "cost", "C:/GeoNet/DEMstreams/test.gdb/CostDist", "", "C:/GeoNet/DEMstreams/test.gdb/Backlink", "", "", "", "", "")‍

The following step was a Cost Path on a single inlet to obtain the route from the inlet to the outlet. Python snippet:

arcpy.gp.CostPath_sa("inlets", "CostDist", "Backlink", "C:/GeoNet/DEMstreams/test.gdb/CostPath_oid01", "EACH_CELL", "OBJECTID")‍

This raster I converted to vector to get the line (not applying simplification):

arcpy.RasterToPolyline_conversion(in_raster="CostPath_oid01", out_polyline_features="C:/GeoNet/DEMstreams/test.gdb/path_oid01a", background_value="ZERO", minimum_dangle_length="0", simplify="NO_SIMPLIFY", raster_field="Value")‍

The resulting featureclass has two features and the one with the longest length (OID=2) is the one we need. However, the data does not have any z values, so I added them with the Interpolate shape 3D tool. Python snippet:

arcpy.InterpolateShape_3d(in_surface="rivdem", in_feature_class="path_oid01a", out_feature_class="C:/GeoNet/DEMstreams/test.gdb/path_oid01a_3d", sample_distance="", z_factor="1", method="BILINEAR", vertices_only="VERTICES_ONLY", pyramid_level_resolution="0")‍

Now we have a line that has z values and we are ready to start the process. The next part is to evaluate each of the vertices of the line and find the corresponding point 1, 2 or 5 km further along the river and compare the heights. For this I created the following code:

def main():
    import arcpy
    import os
    arcpy.env.overwriteOutput = True

    fc = r'C:\GeoNet\DEMstreams\test.gdb\path_oid01a_3d'
    oid_line = 2
    fc_out = r'C:\GeoNet\DEMstreams\test.gdb\result_v05'

    # output fields
    fld_ed1km = "ElevationDrop1km"
    fld_ed2km = "ElevationDrop2km"
    fld_ed5km = "ElevationDrop5km"

    # create output fc
    sr = arcpy.Describe(fc).spatialReference
    ws, fc_name = os.path.split(fc_out)
    arcpy.CreateFeatureclass_management(ws, fc_name, "POINT", None, "ENABLED", "ENABLED", sr)

    # add fields
    arcpy.AddField_management(fc_out, fld_ed1km, "DOUBLE")
    arcpy.AddField_management(fc_out, fld_ed2km, "DOUBLE")
    arcpy.AddField_management(fc_out, fld_ed5km, "DOUBLE")


    # create list M,Z for easy data access
    lst_mz = []
    fld_oid = arcpy.Describe(fc).OIDFieldName
    where = "{0} = {1}".format(fld_oid, oid_line)
    polyline = arcpy.da.SearchCursor(fc, ('SHAPE@'), where).next()[0]
    for part in polyline:
        for pnt in part:
            m = polyline.measureOnLine(pnt, False)
            z = pnt.Z
            lst_mz.append([m, z, pnt]) #

    print "Records to process:", len(lst_mz)
##    print lst_mz
##    exit("bailing out...")

    # start insert cursor
    flds = ('SHAPE@', fld_ed1km, fld_ed2km, fld_ed5km)
    with arcpy.da.InsertCursor(fc_out, flds) as curs:

        # evaluate each point
        cnt = 0
        for lst in lst_mz:
            cnt += 1
            if cnt % 10 == 0:
                print "Processing vertex: ", cnt, drop_1km, drop_2km, drop_5km
            m = lst[0]
            z = lst[1]
            pnt = lst[2]
            drop_1km = GetDropatM(lst_mz, z, m + 1000, cnt - 1)
            drop_2km = GetDropatM(lst_mz, z, m + 2000, cnt - 1)
            drop_5km = GetDropatM(lst_mz, z, m + 5000, cnt - 1)

            pntg = arcpy.PointGeometry(arcpy.Point(pnt.X, pnt.Y, z, m), sr)
            curs.insertRow((pntg, drop_1km, drop_2km, drop_5km, ))


def GetDropatM(lst_mz, z, m_search, i_init):
    # first find the points nearest to given m value (m_search)
    index = i_init
    m = 0
    ok = True
    while m < m_search:
        if index < len(lst_mz):
            lst = lst_mz[index]
            m = lst[0]
            index += 1
        else:
            # no value available
            ok = False
            break
    if ok:
        m_start = lst_mz[index - 2][0]
        m_end = lst_mz[index - 1][0]
        # print "M:", m_search, m_start, m_end
        z_start = lst_mz[index - 2][1]
        z_end = lst_mz[index - 1][1]
        # print "Z:", z, z_start, z_end
        z_search = (m_search - m_start) / (m_end - m_start) * z_start + (m_end - m_search) / (m_end - m_start) * z_end
        return z - z_search

    else:
        return None


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

The result is a featureclass with all the vertices and for each vertice the drop at 1, 2 and 5 km is calculated. The result displayed using the drop in 5 km field:

Drop at 1km:

Re: How to find Elevation difference along the river line?

SepheFox1 — Thu, 30 Jul 2020 14:33:11 GMT

Wow Xander, this is awesome. What a shame the OP never responded to indicate if it was the solution they needed.

Re: How to find Elevation difference along the river line?

XanderBakker — Thu, 30 Jul 2020 16:42:29 GMT

Hi sfox_cpud ,

It sometimes happens that an OP does not respond to a possible solution (you may call it a "hit and run"), but I like to think that a possible solution might be useful not only to the OP but to others too. And in this case, thinking about this solution has helped me to create ideas that have been implemented at clients. So I don't feel I have wasted my time in any way. Thank for the shout out!