diff --git a/history.asc b/history.asc index 8f6eedd..178daa0 100644 --- a/history.asc +++ b/history.asc @@ -35,19 +35,27 @@ == 0.7 (2022-09-13) -* This includes a number of improvements to the valley width extraction algorithms including those used in the paper https://doi.org/10.5194/esurf-10-437-2022 -* This release has a new dependency on the OpenCV library. The programs will compile without this library but the valley metrics component will only compile with OpenCV. -* This version removed dependency on the Point Cloud Library. -* It also includes updated functionality for the cosmogenic tools, so the analyses from the CAIRN (Mudd et al 2016) method are streamlined. +* This includes a number of improvements to the valley width extraction algorithms including those used in the paper https://doi.org/10.5194/esurf-10-437-2022 (SMM and FJC) +* This release has a new dependency on the OpenCV library. The programs will compile without this library but the valley metrics component will only compile with OpenCV. (FJC) +* This version removed dependency on the Point Cloud Library. (FJC) +* It also includes updated functionality for the cosmogenic tools, so the analyses from the CAIRN (Mudd et al 2016) method are streamlined. (SMM) == 0.8 (2023-05-03) -* New contributions from SMM and FJC -* Additions to channel extraction include getting longest channel, including more information in the output files, and an option to print channels as shapefiles. -* A number of changes to the cosmo tool including better accounting for a column, an interface to bring in transient erosion rates, better support for CRN, smoother operation of the nesting functions, and fixes to the shielding calculations that allow you to compute erosion rates in one step. -* Added a function for creating steady state fluvial landscapes -* Included some new metrics, most notably the normalised concavity index or NCI that was used in Chen et al Nature paper. -* Added a routine that facilitates tagging of major drainage divides. -* Added perona-malik filter to hillslope analysis -* added a small function to analyse channel tips for a channel extraction routine -* Other minor bug and typo fixes. \ No newline at end of file +* Additions to channel extraction include getting longest channel, including more information in the output files, and an option to print channels as shapefiles. (SMM) +* A number of changes to the cosmo tool including better accounting for a column, an interface to bring in transient erosion rates, better support for CRN, smoother operation of the nesting functions, and fixes to the shielding calculations that allow you to compute erosion rates in one step. (SMM) +* Added a function for creating steady state fluvial landscapes (SMM) +* Included some new metrics, most notably the normalised concavity index or NCI that was used in Chen et al Nature paper. (SMM) +* Added a routine that facilitates tagging of major drainage divides. (SMM) +* Added perona-malik filter to hillslope analysis (FJC) +* added a small function to analyse channel tips for a channel extraction routine (SMM) +* Other minor bug and typo fixes. (SMM and FJC) + +== 0.9 (2023-06-23) + +* Added function in lsdtt-basic-metrics to extract basin outlines (SMM) +* Added function in lsdtt-basic-metrics to get the bearing of channel segments (SMM) +* Added function in lsdtt-basic-metrics to print flow direction codes in ArcMap format (SMM) +* More testing of lsdtt-valley-metrics to reduce bugs and give more informative error messages. (SMM) +* Added swath output to the terrace routine (SMM) +* Some testing of ridge extraction and minor bug fixes (FJC) \ No newline at end of file diff --git a/src/LSDBasin.cpp b/src/LSDBasin.cpp index b459a8d..d016284 100644 --- a/src/LSDBasin.cpp +++ b/src/LSDBasin.cpp @@ -873,7 +873,7 @@ void LSDBasin::set_Perimeter(LSDFlowInfo& FlowInfo) { FlowInfo.retrieve_current_row_and_col(BasinNodes[q], i, j); - BasinData[i][j] = BasinNodes[q]; + BasinData[i][j] = BasinNodes[q]; } @@ -884,31 +884,31 @@ void LSDBasin::set_Perimeter(LSDFlowInfo& FlowInfo) NDVCount = 0; // cout << i << "||" << j << endl; - if (i != 0 && j != 0 && i= 1 && NDVCount < 8) - { //increase the first value to get a simpler polygon (changed to 1 by FJC 23/03/15 to get only internal hilltops. - //edge pixel // Otherwise not all external ridges were being excluded from the analysis). - I.push_back(i); - J.push_back(j); - B.push_back(BasinNodes[q]); - } - } - else - { - ++i; + if (i != 0 && j != 0 && i= 1 && NDVCount < 8) + { //increase the first value to get a simpler polygon (changed to 1 by FJC 23/03/15 to get only internal hilltops. + //edge pixel // Otherwise not all external ridges were being excluded from the analysis). + I.push_back(i); + J.push_back(j); + B.push_back(BasinNodes[q]); } + } + else + { + ++i; + } } @@ -916,8 +916,6 @@ void LSDBasin::set_Perimeter(LSDFlowInfo& FlowInfo) Perimeter_i = I; Perimeter_j = J; Perimeter_nodes = B; - - } //=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- @@ -967,10 +965,6 @@ void LSDBasin::print_perimeter_hypsometry_to_csv(LSDFlowInfo& FlowInfo, string p set_Perimeter(FlowInfo); clean_perimeter(FlowInfo); - // TESTING - string perim_test = "/home/clubb/Data_for_papers/drainage_capture/Santa_Cruz/FUCK_THIS.csv"; - print_perimeter_to_csv(FlowInfo, perim_test); - // open the file ofstream perim_out; perim_out.open(perimeter_fname.c_str()); @@ -1100,10 +1094,10 @@ vector LSDBasin::order_perimeter_nodes(LSDFlowInfo& FlowInfo) next_j = Outlet_j+1; cout << "The closest node is to the SE" << endl; } - else if (PerimeterNodes[Outlet_i-1][Outlet_j+1] == 1) // southwest + else if (PerimeterNodes[Outlet_i+1][Outlet_j-1] == 1) // southwest { - next_i = Outlet_i-1; - next_j = Outlet_j+1; + next_i = Outlet_i+1; + next_j = Outlet_j-1; cout << "The closest node is to the SW" << endl; } else @@ -1357,6 +1351,209 @@ vector LSDBasin::order_perimeter_nodes(LSDFlowInfo& FlowInfo) return sorted_nodes; } + + + + + + + +//=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- +// Order perimeter nodes from the outlet +// Must CLEAN the perimeter first using clean_perimeter function +// FJC 16/01/18 +//=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- +vector LSDBasin::order_perimeter_nodes_SMM(LSDFlowInfo& FlowInfo) +{ + // first clean the perimeter nodes + //clean_perimeter(FlowInfo); + cout << "Ordering the perimeter nodes..." << endl; + vector sorted_nodes; + Array2D PerimeterNodes(NRows,NCols,0); + + // first get the outlet node + int outlet_node = get_Outlet_node(); + int outlet_row, outlet_col; + FlowInfo.retrieve_current_row_and_col(outlet_node, outlet_row, outlet_col); + Outlet_i = outlet_row; + Outlet_j = outlet_col; + cout << "The outlet node is: " << outlet_node << endl; + + // get an array of perimeter nodes + for (int i = 0; i < int(Perimeter_nodes.size()); i++) + { + int this_row, this_col; + FlowInfo.retrieve_current_row_and_col(Perimeter_nodes[i], this_row, this_col); + PerimeterNodes[this_row][this_col] = 1; + } + // The outlet always needs to be in the perimeter + // The outlet is already added so it gets an increment to 2 + PerimeterNodes[Outlet_i][Outlet_j] = 2; + cout << "Got the array of perimeter nodes" << endl; + cout << "The outlet row and column are: " << Outlet_i << "," << Outlet_j << endl; + + // push back the outlet node, node 0 + sorted_nodes.push_back(outlet_node); + int next_i, next_j; + bool first_node = true; // bool to check if you're at the first node. This makes a difference because if you're at the first node we don't want to go back to the outlet. + + + // N, S, E, and W will always be the shortest distances, so do these first + if (PerimeterNodes[Outlet_i][Outlet_j-1] == 1) // West + { + next_i = Outlet_i; + next_j = Outlet_j-1; + //cout << "The closest node is to the west" << endl; + } + else if (PerimeterNodes[Outlet_i-1][Outlet_j] == 1) // North + { + next_i = Outlet_i-1; + next_j = Outlet_j; + //cout << "The closest node is to the north" << endl; + } + else if (PerimeterNodes[Outlet_i][Outlet_j+1] == 1) // east + { + next_i = Outlet_i; + next_j = Outlet_j+1; + //cout << "The closest node is to the east" << endl; + } + else if (PerimeterNodes[Outlet_i+1][Outlet_j] == 1) // south + { + next_i = Outlet_i+1; + next_j = Outlet_j; + //cout << "The closest node is to the south" << endl; + } + else if (PerimeterNodes[Outlet_i-1][Outlet_j-1] == 1) // northwest + { + next_i = Outlet_i-1; + next_j = Outlet_j-1; + //cout << "The closest node is to the NW" << endl; + } + else if (PerimeterNodes[Outlet_i-1][Outlet_j+1] == 1) // northeast + { + next_i = Outlet_i-1; + next_j = Outlet_j+1; + //cout << "The closest node is to the NE" << endl; + } + else if (PerimeterNodes[Outlet_i+1][Outlet_j+1] == 1) // southeast + { + next_i = Outlet_i+1; + next_j = Outlet_j+1; + //cout << "The closest node is to the SE" << endl; + } + else if (PerimeterNodes[Outlet_i+1][Outlet_j-1] == 1) // southwest + { + next_i = Outlet_i+1; + next_j = Outlet_j-1; + //cout << "The closest node is to the SW" << endl; + } + else + { + //cout << "None of these were perimeter nodes, oops" << endl; + } + + cout << "The first node is at: " << next_i << "," << next_j << endl; + + // push back the next node to the sorted node vector + int next_node = FlowInfo.retrieve_node_from_row_and_column(next_i, next_j); + sorted_nodes.push_back(next_node); + + bool reached_outlet = false; + int this_i, this_j; + // now start at the outlet node and find the nearest perimeter node. + while (reached_outlet == false) + { + // start at the next node and find the one with the closest distance that + // hasn't already been visited + this_i = next_i; + this_j = next_j; + PerimeterNodes[this_i][this_j] = 2; + + // N, S, E, and W will always be the shortest distances, so do these first + if (PerimeterNodes[this_i][this_j-1] == 1) // West + { + next_i = this_i; + next_j = this_j-1; + //cout << "The closest node is to the west" << endl; + } + else if (PerimeterNodes[this_i-1][this_j] == 1) // North + { + next_i = this_i-1; + next_j = this_j; + //cout << "The closest node is to the north" << endl; + } + else if (PerimeterNodes[this_i][this_j+1] == 1) // east + { + next_i = this_i; + next_j = this_j+1; + //cout << "The closest node is to the east" << endl; + } + else if (PerimeterNodes[this_i+1][this_j] == 1) // south + { + next_i = this_i+1; + next_j = this_j; + //cout << "The closest node is to the south" << endl; + } + else if (PerimeterNodes[this_i-1][this_j-1] == 1) // northwest + { + next_i = this_i-1; + next_j = this_j-1; + //cout << "The closest node is to the NW" << endl; + } + else if (PerimeterNodes[this_i-1][this_j+1] == 1) // northeast + { + next_i = this_i-1; + next_j = this_j+1; + //cout << "The closest node is to the NE" << endl; + } + else if (PerimeterNodes[this_i+1][this_j+1] == 1) // southeast + { + next_i = this_i+1; + next_j = this_j+1; + //cout << "The closest node is to the SE" << endl; + } + else if (PerimeterNodes[this_i+1][this_j-1] == 1) // southwest + { + next_i = this_i+1; + next_j = this_j-1; + //cout << "The closest node is to the SW" << endl; + } + else + { + reached_outlet = true; + } + + //cout << "The next node is at: " << next_i << "," << next_j << endl; + + if (!reached_outlet) + { + next_node = FlowInfo.retrieve_node_from_row_and_column(next_i, next_j); + sorted_nodes.push_back(next_node); + } + } + + return sorted_nodes; +} + + + + + + + + + + + + + + + + + + + + //=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- // Set the four different hillslope length measurements for the basin. @@ -2527,7 +2724,7 @@ void LSDBasin::organise_perimeter(LSDFlowInfo& flowpy) // TESTING FUNCTION, DELETE IT AFTERWARDS BORIS!!!! Perimeter_nodes = Perimeter_nodes_sorted; - print_perimeter_to_csv(flowpy, "/home/boris/Desktop/LSD/capture/sorbas/peritest_AFTER.csv"); + print_perimeter_to_csv(flowpy, "/home/boris/Desktop/LSD/capture/sorbas/peritest_AFTER.csv"); diff --git a/src/LSDBasin.hpp b/src/LSDBasin.hpp index bcc49e8..ff10428 100644 --- a/src/LSDBasin.hpp +++ b/src/LSDBasin.hpp @@ -464,6 +464,12 @@ class LSDBasin /// @date 16/01/18 vector order_perimeter_nodes(LSDFlowInfo& FlowInfo); + /// @brief Orders perimeter nodes from the outlet. Uses a different approach + /// @param FlowInfo the LSDFlowInfo object + /// @author SMM + /// @date 29/05/23 + vector order_perimeter_nodes_SMM(LSDFlowInfo& FlowInfo); + /// @brief Set the four different hillslope length measurements for the basin. /// @param FlowInfo Flowinfo object. /// @param HillslopeLengths LSDRaster of hillslope lengths from the hilltop flow routing method. diff --git a/src/LSDJunctionNetwork.cpp b/src/LSDJunctionNetwork.cpp index 46317ae..ed11637 100644 --- a/src/LSDJunctionNetwork.cpp +++ b/src/LSDJunctionNetwork.cpp @@ -2002,7 +2002,400 @@ void LSDJunctionNetwork::calculate_junction_angles_complete(vector Junction +//=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- +// This the bearings and gradients of each channel segment +// +// The vector of floats actually contains a number of elements +// The int map has: +// [0] Stream order of junction +// The float map has +// [0] drainage area junction +// [1] bearing between junction and penultimate node +// [2] bearing between junction and vertical interval node +// [3] x junction +// [4] y junction +// [5] z junction +// [6] x penultimate +// [7] y penultimate +// [8] z penultimate +// [9] x vertical interval +// [10] y vertical interval +// [11] z vertical interval +// [12] FD junction +// [13] FD penultimate +// [14] FD vertical interval +//=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- +void LSDJunctionNetwork::calculate_segment_bearings_and_gradients_complete(vector JunctionList, + LSDFlowInfo& FlowInfo, LSDRaster& Elevation, + LSDRaster& FlowDistance, float vertical_interval, + map >& JA_int_info, + map >& JA_float_info ) +{ + cout << "I'm calculating complete bearings of all channel segments associated with your chosen junctions." << endl; + + map > JA_float_map; + map > JA_int_map; + vector temp_float_junctioninfo; + vector temp_int_junctioninfo; + vector nine_element_float(17,0); + vector four_element_int(1,0); + + int NJuncs = int(JunctionList.size()); + if (NJuncs == 0) + { + NJuncs = get_NJunctions(); + cout << "You gave me an empty junction list, I am getting bearings and gradients for all segments." << endl; + vector JL; + for(int i = 0; i= NJunctions) + { + cout << "FATAL ERROR LSDJunctionNetwork::calculate_segment_bearings_and_gradients_complete" << endl; + cout << "You have called a junction that doesn't exist." << endl; + exit(EXIT_FAILURE); + } + + + // get the receiver junction and the junction node + receiver_junc = get_Receiver_of_Junction(this_junc); + junction_node = get_Node_of_Junction(this_junc); + + //cout << "Junction " << junc << " of " << NJuncs; + //cout << ", this junc: " << this_junc << " and receiver is: " << receiver_junc << endl; + + + + // if not a baselevel see if it has donors + if (this_junc == receiver_junc) + { + //cout << ", this is a baselevel junction." << endl; + //JunctionAngles.push_back(NoDataValue); + } + else + { + // Get the junction order. + junction_order = get_StreamOrder_of_Junction(this_junc); + //cout << "Got junction order, is: " << junction_order << endl; + + // Now the pu node: + end_node_penultimate = get_penultimate_node_from_stream_link(this_junc,FlowInfo); + //cout << "The penultamite end node is: " << end_node_penultimate << endl; + + // we need to deal with an edge case where the next junction is one pixel from the previous junction + if (junction_node!=end_node_penultimate) + { + // reset the temp vectors + temp_float_junctioninfo = nine_element_float; + temp_int_junctioninfo = four_element_int; + //cout << "reset vectors" << endl; + + // get all the elevations and flow distances + // This is for the junction itself + FlowInfo.retrieve_current_row_and_col(junction_node,row,col); + FD_J = FlowDistance.get_data_element(row,col); + z_J = Elevation.get_data_element(row,col); + FlowInfo.get_x_and_y_locations(row, col, x_j, y_j); + //cout << "Got some locations and flow distances for the junction " << endl; + + // this is for the penultimate node + FlowInfo.retrieve_current_row_and_col(end_node_penultimate,row,col); + FD_PU = FlowDistance.get_data_element(row,col); + z_PU = Elevation.get_data_element(row,col); + FlowInfo.get_x_and_y_locations(row, col, x_pu, y_pu); + //cout << "Got some locations and flow distances for the PU node " << endl; + + //cout << ", FD_J: " << FD_J << ", FD_PU: " << FD_PU; + + // get the gradient + grad_PU = (z_J-z_PU)/(FD_J-FD_PU); + + // We now go through each donor and the receiver to get the fixed drop + float z_search; + + // Initial logic: if the vertical interval is greater than the vertical drop to the penulimate node, + // we use that node to calculate vertical interval, etc. + if ( (z_J-z_PU) <= vertical_interval) + { + z_VI = z_PU; + FD_VI = FD_PU; + grad_VI = grad_PU; + x_vi = x_pu; + y_vi = y_pu; + } + else + { + // This is repeated in case the if statement for the vertical interval is + // not triggered. Ensures this component returns some value. + z_VI = z_PU; + FD_VI = FD_PU; + grad_VI = grad_PU; + + // Logic if the fixed vertical interval is somewhere within the channel segment + int this_node = junction_node; + int next_node; + + while(this_node != end_node_penultimate) + { + FlowInfo.retrieve_receiver_information(this_node,next_node, row, col); + + // Now check the elevation of the search node + z_search = Elevation.get_data_element(row,col); + + if ((z_J - z_search) <= vertical_interval) + { + z_VI = z_search; + FD_VI = FlowDistance.get_data_element(row,col); + grad_VI = (z_J-z_search)/(FD_J-FD_VI); + FlowInfo.get_x_and_y_locations(row, col, x_vi, y_vi); + this_node = next_node; + } + else + { + this_node = next_node; + } + } + } + //cout << ", done with this junction" << endl; + + + // now we get the bearings + JPU_bearing = clockwise_angle_between_vector_and_north(x_j, y_j, x_pu, y_pu ); + JVI_bearing = clockwise_angle_between_vector_and_north(x_j, y_j, x_vi, y_vi ); + + DA_J = FlowInfo.get_DrainageArea_square_m(junction_node); + + temp_int_junctioninfo[0] = junction_order; + + temp_float_junctioninfo[0] = DA_J; + temp_float_junctioninfo[1] = JPU_bearing; + temp_float_junctioninfo[2] = JVI_bearing; + temp_float_junctioninfo[3] = grad_PU; + temp_float_junctioninfo[4] = grad_VI; + temp_float_junctioninfo[5] = x_j; + temp_float_junctioninfo[6] = y_j; + temp_float_junctioninfo[7] = z_J; + temp_float_junctioninfo[8] = x_pu; + temp_float_junctioninfo[9] = y_pu; + temp_float_junctioninfo[10] = z_PU; + temp_float_junctioninfo[11] = x_vi; + temp_float_junctioninfo[12] = x_vi; + temp_float_junctioninfo[13] = z_VI; + temp_float_junctioninfo[14] = FD_J; + temp_float_junctioninfo[15] = FD_PU; + temp_float_junctioninfo[16] = FD_VI; + + JA_int_map[junc] = temp_int_junctioninfo; + JA_float_map[junc] = temp_float_junctioninfo; + + } + } // end if statement that make sure this junction is not a baselevel + } // end loop through junctions + + JA_int_info = JA_int_map; + JA_float_info = JA_float_map; +} + + + +//=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- +// +// This prints junction angles to a csv file +// Uses the complete junction angle code so much more extensive statistics +// +//=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- +void LSDJunctionNetwork::print_complete_segment_bearings_and_gradients_to_csv(vector JunctionList, + LSDFlowInfo& FlowInfo, LSDRaster& Elevations, + LSDRaster& FlowDistance, float vertical_interval, + string csv_name) +{ + + ofstream csv_out; + csv_out.open(csv_name.c_str()); + csv_out.precision(9); + + + // get the junction information + map > JuncInfo_int; + map > JuncInfo_float; + + cout << "Now I will calculate the bearings and segment gradients." << endl; + calculate_segment_bearings_and_gradients_complete(JunctionList, FlowInfo, Elevations, + FlowDistance, vertical_interval, + JuncInfo_int,JuncInfo_float); + map >::iterator iter; + vector this_JI_float; + vector this_JI_int; + int this_junc; + int this_node,curr_row,curr_col; + double latitude, longitude; + + csv_out << "latitude,longitude,junction_number,"; + csv_out << "junction_stream_order,"; + csv_out << "junction_drainage_area,segment_bearing_whole_segment, segment_bearing_vertical_interval,"; + csv_out << "gradient_whole_segment,gradient_vertical_interval,"; + csv_out << "x_j,y_j,z_j,x_pu,y_pu,z_pu,x_vi,y_pi,z_vi,"; + csv_out << "FlowDistance_junction,FlowDistance_pu,FlowDistance_vi" << endl; + + cout << "csv Header printed, moving on to looping through the junctions" << endl; + for(iter = JuncInfo_float.begin(); iter != JuncInfo_float.end(); ++iter) + { + this_junc = iter->first; + this_JI_float = iter->second; + this_JI_int = JuncInfo_int[this_junc]; + + // get the row and column of the junction from the junction node + this_node = JunctionVector[this_junc]; + LSDCoordinateConverterLLandUTM Converter; + FlowInfo.retrieve_current_row_and_col(this_node,curr_row,curr_col); + get_lat_and_long_locations(curr_row, curr_col, latitude,longitude, Converter); + // print to the csv file + csv_out << latitude <<"," << longitude <<"," << this_junc <<"," + << this_JI_int[0] << "," + << this_JI_float[0] << "," << deg(this_JI_float[1]) << "," << deg(this_JI_float[2]) << "," + << this_JI_float[3] << "," << this_JI_float[4] << "," << this_JI_float[5] << "," + << this_JI_float[6] << "," << this_JI_float[7] << "," << this_JI_float[8] << "," + << this_JI_float[9] << "," << this_JI_float[10] << "," << this_JI_float[11] << "," + << this_JI_float[12] << "," << this_JI_float[13] << "," << this_JI_float[14] << "," + << this_JI_float[15] << "," << this_JI_float[16] << endl; + } + + csv_out.close(); +} + + + + +void LSDJunctionNetwork::print_complete_segment_bearings_and_gradients_to_geojson(vector JunctionList, + LSDFlowInfo& FlowInfo, LSDRaster& Elevations, + LSDRaster& FlowDistance, float vertical_interval, + string json_name) +{ + + // open the geojson for orthogonal lines + ofstream bearing_lines; + bearing_lines.precision(9); + bearing_lines.open((json_name).c_str()); + + // get the junction information + map > JuncInfo_int; + map > JuncInfo_float; + + cout << "Now I will calculate the complete junction angle information." << endl; + calculate_segment_bearings_and_gradients_complete(JunctionList, FlowInfo, Elevations, + FlowDistance, vertical_interval, + JuncInfo_int,JuncInfo_float); + map >::iterator iter; + vector this_JI_float; + vector this_JI_int; + int this_junc; + int this_node,curr_row,curr_col; + double latitude, longitude; + + + bearing_lines << "{" << endl; + bearing_lines << "\"type\": \"FeatureCollection\"," << endl; + bearing_lines << "\"crs\": { \"type\": \"name\", \"properties\": { \"name\": \"urn:ogc:def:crs:OGC:1.3:CRS84\" } }," << endl; + bearing_lines << "\"features\": [" << endl; + + LSDCoordinateConverterLLandUTM Converter; + + int count = 0; + cout << "csv Header printed, moving on to looping through the junctions" << endl; + for(iter = JuncInfo_float.begin(); iter != JuncInfo_float.end(); ++iter) + { + count++; + + this_junc = iter->first; + this_JI_float = iter->second; + this_JI_int = JuncInfo_int[this_junc]; + + // this is for latitude and longitude + double lat_j, long_j; + double lat_pu, long_pu; + FlowInfo.get_lat_and_long_locations(this_JI_float[5], this_JI_float[6], lat_j, long_j, Converter); + FlowInfo.get_lat_and_long_locations(this_JI_float[8], this_JI_float[9], lat_pu, long_pu, Converter); + + // now write the geojson with orthogonal lines + string first_bit = "{ \"type\": \"Feature\", \"properties\": { \"full_segment_bearing\": "; + string second_bit = ", \"vi_segment_bearing\": "; + string third_bit = ", \"full_segment_gradient\": "; + string fourth_bit = ", \"vi_gradient\": "; + string fifth_bit = ", \"stream_order\": "; + string sixth_bit = ", \"drainage_area\": "; + + string twelvth_bit = " }, \"geometry\": { \"type\": \"Linestring\", \"coordinates\": [ ["; + + string last_bit = " ] ] } }"; + if (count == 1) + { + bearing_lines << first_bit << deg(this_JI_float[1]) << second_bit << deg(this_JI_float[2]) + << third_bit << this_JI_float[3] << fourth_bit << this_JI_float[4] << fifth_bit << this_JI_int[0] << sixth_bit << this_JI_float[0] + << twelvth_bit << long_j << "," << lat_j << "], [" << long_pu << "," << lat_pu << last_bit + << endl; + } + else + { + bearing_lines << "," << first_bit << deg(this_JI_float[1]) << second_bit << deg(this_JI_float[2]) + << third_bit << this_JI_float[3] << fourth_bit << this_JI_float[4] << fifth_bit << this_JI_int[0] << sixth_bit << this_JI_float[0] + << twelvth_bit << long_j << "," << lat_j << "], [" << long_pu << "," << lat_pu << last_bit + << endl; + } + + } + + bearing_lines << "]" << endl; + bearing_lines << "}" << endl; + + bearing_lines.close(); + +} @@ -2152,6 +2545,8 @@ void LSDJunctionNetwork::print_junction_angles_from_basin_list(vector Junct csv_out.close(); } + + //=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- // // This gets the junction angle stats for all basins of a given size @@ -2545,7 +2940,13 @@ LSDIndexChannel LSDJunctionNetwork::generate_longest_index_channel_from_junction } //=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= - +vector LSDJunctionNetwork::generate_longest_nodeindex_channel_from_junction(int outlet_junction, LSDFlowInfo& FInfo, + LSDRaster& dist_from_outlet) +{ + LSDIndexChannel this_channel = generate_longest_index_channel_from_junction(outlet_junction, FInfo,dist_from_outlet); + vector nodeindex_chan = this_channel.get_NodeSequence(); + return nodeindex_chan; +} //=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= //=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= @@ -2600,6 +3001,15 @@ LSDIndexChannel LSDJunctionNetwork::generate_longest_index_channel_in_basin(int } //=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=--=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= + +vector LSDJunctionNetwork::generate_longest_nodeindex_channel_in_basin(int basin_junction, LSDFlowInfo& FInfo, + LSDRaster& dist_from_outlet) +{ + LSDIndexChannel this_channel = generate_longest_index_channel_in_basin(basin_junction, FInfo,dist_from_outlet); + vector nodeindex_chan = this_channel.get_NodeSequence(); + return nodeindex_chan; +} + //=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= //=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= // this function takes a vector of basin junctions and returns a vector of the farthest @@ -5527,6 +5937,22 @@ map LSDJunctionNetwork::ExtractAllRidges(LSDFlowInfo& FlowInfo) // we use the stack to get the junction ordering vector pen_nodes = get_node_list_of_penultimate_node_from_junction_list(SVector,FlowInfo); + // checking on these vectors + if (upslope_junc_nodes.size() != pen_nodes.size()) + { + cout << "Warning LSDJunctionNetwork::ExtractAllRidges you have different sizes of junction nodes and penultamite nodes" << endl; + } + else + { + for (int jn = 0; jn< int(upslope_junc_nodes.size()); jn++) + { + if (upslope_junc_nodes[jn] == pen_nodes[jn]) + { + cout << "Warning, junction node " << jn << " is at the same place as the penultimate node and this will likeley cause an error" << endl; + } + } + } + // get the stream orders of these vectors vector so_stack; int row,col; @@ -5543,10 +5969,11 @@ map LSDJunctionNetwork::ExtractAllRidges(LSDFlowInfo& FlowInfo) // for debugging //cout << "O: " << StreamOrderVector[ SVector[i] ] << endl; - if(pen_nodes[i] != NoDataValue && pen_nodes[i] > 0 && pen_nodes[i] < NDataNodes-1) + //if(pen_nodes[i] != NoDataValue && pen_nodes[i] >= 0 && pen_nodes[i] <= NDataNodes-1) + if(pen_nodes[i] != NoDataValue) { upslope_junc_nodes = FlowInfo.get_upslope_nodes(pen_nodes[i]); - //cout << "Number of upslope nodes: " << int(upslope_junc_nodes.size()) << endl; + cout << "Number of upslope nodes: " << int(upslope_junc_nodes.size()) << endl; // now get the nodes that are internal vector is_internal = FlowInfo.internal_nodes(upslope_junc_nodes); @@ -5563,6 +5990,7 @@ map LSDJunctionNetwork::ExtractAllRidges(LSDFlowInfo& FlowInfo) //RidgeNetwork[row][col] = StreamOrderVector[ SVector[i] ]; node_map[upslope_junc_nodes[usn]] = StreamOrderVector[ SVector[i] ]; + // TO DO - check if these keys exist. Key = node index of the ridgeline. Value = stream order of the downslope junction } } } @@ -5587,6 +6015,23 @@ map LSDJunctionNetwork::ExtractAllRidges(LSDFlowInfo& FlowInfo, string // we use the stack to get the junction ordering vector pen_nodes = get_node_list_of_penultimate_node_from_junction_list(SVector,FlowInfo); + // checking on these vectors + if (upslope_junc_nodes.size() != pen_nodes.size()) + { + cout << "Warning LSDJunctionNetwork::ExtractAllRidges you have different sizes of junction nodes and penultamite nodes" << endl; + } + else + { + for (int jn = 0; jn< int(upslope_junc_nodes.size()); jn++) + { + if (upslope_junc_nodes[jn] == pen_nodes[jn]) + { + cout << "Warning, junction node " << jn << " is at the same place as the penultimate node and this will likeley cause an error" << endl; + } + } + } + + // get the stream orders of these vectors vector so_stack; int row,col; @@ -11600,9 +12045,9 @@ map LSDJunctionNetwork::filter_nodes_by_distance_to_channel(mapfirst; this_SO = it->second; int chan_node = find_nearest_downslope_channel(this_node, threshold_SO, FlowInfo); - cout << "this node: " << this_node << " chan node: " << chan_node < threshold_dist) { // this node is farther away than the threshold distance, so we keep it diff --git a/src/LSDJunctionNetwork.hpp b/src/LSDJunctionNetwork.hpp index 4816fcf..d193068 100644 --- a/src/LSDJunctionNetwork.hpp +++ b/src/LSDJunctionNetwork.hpp @@ -273,11 +273,16 @@ class LSDJunctionNetwork /// @brief This is a more complete function for junction angles - /// It overwirtes two maps, containing all sorts of information about + /// It overwrites two maps, containing all sorts of information about /// the junction angles. /// This is an overloaded version that gives some slope and flow distance information /// @param JunctionList a list of junctions /// @param FlowInfo an LSDFlowInfo object + /// @param Elevation The elevation of the studied landscape + /// @param FlowDistance The flow distance raster + /// @param vertical_interval A vertical interval over which gradient is measured + /// note that both this interval and the vertical offset between the junction pixel and the last pixel before the + /// next junction are used to generate two gradients. /// @param JA_int_info A map where the key is the junction number /// and the vector is a series of integer data about that juction /// @param JA_float_info A map where the key is the junction number @@ -377,6 +382,60 @@ class LSDJunctionNetwork LSDRaster& FlowDistance, float vertical_interval, string csv_name); + + /// @brief This calculates the channel bearings and gradients so you can make + /// a stereonet using the channel segment information. + /// @param JunctionList a list of junctions + /// @param FlowInfo an LSDFlowInfo object + /// @param Elevation The elevation of the studied landscape + /// @param FlowDistance The flow distance raster + /// @param vertical_interval A vertical interval over which gradient is measured + /// note that both this interval and the vertical offset between the junction pixel and the last pixel before the + /// next junction are used to generate two gradients. + /// @param JA_int_info A map where the key is the junction number + /// and the vector is a series of integer data about that juction + /// @param JA_float_info A map where the key is the junction number + /// and the vector is a series of float data about that juction + /// @author SMM + /// @date 01/06/2023 + void calculate_segment_bearings_and_gradients_complete(vector JunctionList, + LSDFlowInfo& FlowInfo, LSDRaster& Elevation, + LSDRaster& FlowDistance, float vertical_interval, + map >& JA_int_info, + map >& JA_float_info ); + + + /// @brief This prints the bearings and gradients to a csv file + /// @param JunctionList The list of junctions to analyze. If this is an empty vector, + /// the code analyses all junctions in the DEM + /// @param FlowInfo The LSDFlowInfo object + /// @param Elevations An elevation raster + /// @param FlowDistance A flow distance raster + /// @param vertical_interval: the vertical interval around which you want the slope measured + /// @param csv_name The name of the file. Needs full path and csv extension + /// @author SMM + /// @date 01/06/2023 + void print_complete_segment_bearings_and_gradients_to_csv(vector JunctionList, + LSDFlowInfo& FlowInfo, LSDRaster& Elevations, + LSDRaster& FlowDistance, float vertical_interval, + string csv_name); + + + /// @brief This prints the bearings and gradients to a geojson file + /// @param JunctionList The list of junctions to analyze. If this is an empty vector, + /// the code analyses all junctions in the DEM + /// @param FlowInfo The LSDFlowInfo object + /// @param Elevations An elevation raster + /// @param FlowDistance A flow distance raster + /// @param vertical_interval: the vertical interval around which you want the slope measured + /// @param json_name The name of the file. Needs full path and geojson extension + /// @author SMM + /// @date 01/06/2023 + void print_complete_segment_bearings_and_gradients_to_geojson(vector JunctionList, + LSDFlowInfo& FlowInfo, LSDRaster& Elevations, + LSDRaster& FlowDistance, float vertical_interval, + string json_name); + /// @brief This gets the junction number of a given node. /// @param Node /// @param FlowInfo Flow Info object @@ -754,6 +813,17 @@ class LSDJunctionNetwork /// @date 01/09/12 LSDIndexChannel generate_longest_index_channel_from_junction(int outlet_junction,LSDFlowInfo& FInfo,LSDRaster& dist_from_outlet); + /// @brief This function extracts the longest channel originating from a junction number + /// outlet_junction. + /// @param outlet_junction Outlet of junction. + /// @param FInfo LSDFlowInfo object. + /// @param dist_from_outlet Distance from outlet junction. + /// @return nodeindex sequence of the channel + /// @author SMM + /// @date 11/06/23 + vector generate_longest_nodeindex_channel_from_junction(int outlet_junction, LSDFlowInfo& FInfo, + LSDRaster& dist_from_outlet); + // this generates the longest channel in a basin. The basin starts where a channel of // some order intersects with a channel of higher order. So the bain includes the // basin junction, but also the channel flowing downstream from this basin @@ -777,6 +847,30 @@ class LSDJunctionNetwork LSDIndexChannel generate_longest_index_channel_in_basin(int basin_junction, LSDFlowInfo& FInfo, LSDRaster& dist_from_outlet); + + // this generates the longest channel in a basin. The basin starts where a channel of + // some order intersects with a channel of higher order. So the bain includes the + // basin junction, but also the channel flowing downstream from this basin + // junction + // It starts from the node of the receiver junction, so if one were to extract + // the basin from this node one would get a basin that starts one node upstream from + // the lowest node in this + /// @brief This generates the longest channel in a basin. + /// + /// @details The basin starts where a channel of some order intersects with a + /// channel of higher order. So the bain includes the basin junction, but also + /// the channel flowing downstream from this basin junction. It starts from the + /// node of the receiver junction, so if one were to extract the basin from + /// this node one would get a basin that starts one node upstream from the lowest node in this. + /// @param basin_junction + /// @param FInfo LSDFlowInfo object. + /// @param dist_from_outlet + /// @return nodeindex of the longest channel. + /// @author SMM + /// @date 11/06/23 + vector generate_longest_nodeindex_channel_in_basin(int basin_junction, LSDFlowInfo& FInfo, + LSDRaster& dist_from_outlet); + /// @brief This generates the upstream source nodes from a vector of basin junctions /// /// @details The basin starts where a channel of some order intersects with a diff --git a/src/LSDLithoCube.cpp b/src/LSDLithoCube.cpp index 34f3d92..678ba30 100644 --- a/src/LSDLithoCube.cpp +++ b/src/LSDLithoCube.cpp @@ -441,9 +441,7 @@ void LSDLithoCube::ingest_litho_data(string filename) // get the header. Print and discard. string str; getline(data_in,str); - //getline(data_in,str); cout << "The header is: " << str << endl; - getline(data_in,str); cout << "Let me get the data for you." << endl; cout << "If this is a big file loading it will take a little time." << endl; @@ -480,6 +478,42 @@ void LSDLithoCube::ingest_litho_data(string filename) } +// Get a count of all the strat keys in the lithocube. +map LSDLithoCube::get_strat_key_counts() +{ + + map strat_map; + int this_strat; + + // get the counts and put them in a map + for (int row = 0; row::iterator it = strat_map.begin(); it != strat_map.end(); ++it) + { + cout << "Strat_code: " << it->first << " with count: " << it->second << endl; + } + return strat_map; +} + + //=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- // // This prints a layter as a raster @@ -1163,7 +1197,8 @@ void LSDLithoCube::convert_EC_to_K_file(vector K_values, string EC_to_K_o EC_to_K_out << "13," << K_values[3] << endl; EC_to_K_out << "14," << K_values[2] << endl; EC_to_K_out << "15," << K_values[3] << endl; - EC_to_K_out << "-9999," << K_values[2] << endl; + EC_to_K_out << "-9999," << K_values[2] << endl; + EC_to_K_out << "-99999," << K_values[2] << endl; EC_to_K_out.close(); } diff --git a/src/LSDLithoCube.hpp b/src/LSDLithoCube.hpp index 0c1c23c..e2a78c8 100644 --- a/src/LSDLithoCube.hpp +++ b/src/LSDLithoCube.hpp @@ -161,6 +161,11 @@ class LSDLithoCube /// @date 25/01/2020 void ingest_stratigraphic_info(string filename); + /// @brief This counts and returns the number of pixels in the lithocube with each stratigraphic code + /// @author SMM + /// @param return a map holding the counts of each strategraphic code + /// @date 02/06/2023 + map get_strat_key_counts(); /// @brief This gets the row and column of a particular point /// @param x_location the x location of a particular point in local coordinates diff --git a/src/lsdtt-drivers/lsdtt-basic-metrics.cpp b/src/lsdtt-drivers/lsdtt-basic-metrics.cpp index 03856af..8c0dd95 100644 --- a/src/lsdtt-drivers/lsdtt-basic-metrics.cpp +++ b/src/lsdtt-drivers/lsdtt-basic-metrics.cpp @@ -68,7 +68,7 @@ int main (int nNumberofArgs,char *argv[]) { - string version_number = "0.8"; + string version_number = "0.9"; string citation = "http://doi.org/10.5281/zenodo.4577879"; cout << "=========================================================" << endl; @@ -241,7 +241,7 @@ int main (int nNumberofArgs,char *argv[]) help_map["swath_points_csv"] = { "string","swath.csv","The name of the csv file for swath mapping. It can be the start and (optional) end point of the swath channel or a series of points that form a polyline that serves as the baseline for the swath.","Filename needs to include the csv extension and have latitude and longitude in the column headers."}; float_default_map["swath_point_spacing"] = 500; - help_map["swath_point_spacing"] = { "float","500","Only used by calcualte_swath_along_line. How closely spaced the points are (in metres) along the baseline.","If this number is large (many times the DEM resolution) the swath will have an irregular edge so this should be close to the DEM resolution."}; + help_map["swath_point_spacing"] = { "float","500","Only used by calculate_swath_along_line. How closely spaced the points are (in metres) along the baseline.","If this number is large (many times the DEM resolution) the swath will have an irregular edge so this should be close to the DEM resolution."}; float_default_map["swath_bin_spacing"] = 1000; help_map["swath_bin_spacing"] = { "float","1000","The spacing of the bins in metres that are used to calculate statistics along the swath.","Should be at least a few times bigger than swath_point_spacing or the pixel size of the DEM."}; @@ -394,7 +394,10 @@ int main (int nNumberofArgs,char *argv[]) help_map["print_directional_gradients"] = { "bool","false","Prints two rasters: the x and y direction gradients after polynomial fitting.","Part of surface fitting metrics."}; bool_default_map["calculate_basin_statistics"] = false; - help_map["calculate_basin_statistics"] = { "bool","false","Prints some statistics for each basin. SMM needs to check what this does.","Part of surface fitting metrics."}; + help_map["calculate_basin_statistics"] = { "bool","false","Prints some statistics for each basin. does mean, median, MAD and other stuff.","You should supply basin locations."}; + + bool_default_map["print_basin_perimeters"] = false; + help_map["print_basin_perimeters"] = {"bool","false","Prints some basin perimeters with elevations.","You will get out an odered csv file."}; //======================================================== // Window size estimation @@ -514,6 +517,9 @@ int main (int nNumberofArgs,char *argv[]) bool_default_map["print_stream_order_raster"] = false; help_map["print_stream_order_raster"] = { "bool","false","Prints a raster with _SO in filename with stream orders of channel in the appropriate pixel.","Generates a big file so we suggest printing the network to csv."}; + bool_default_map["print_flow_direction_raster"] = false; + help_map["print_flow_direction_raster"] = { "bool","false","Prints a raster with the arcmap flow direction codes.","Codes are 32 64 128 for NW N NE, 16 0 1 for W none W, 8 4 2 for SW S SE."}; + bool_default_map["print_channels_to_csv"] = false; help_map["print_channels_to_csv"] = { "bool","false","Prints the channel network to a csv file.","This version produces smaller files than the raster version."}; @@ -674,7 +680,14 @@ int main (int nNumberofArgs,char *argv[]) help_map["print_junction_angles_to_csv_in_basins"] = { "bool","false","Prints a csv with the locations of the junctions within selected basins and associated statistics.","csv file contains junction angles; bending angles; areas; gradients; and other information."}; float_default_map["SA_vertical_interval"] = 10; - help_map["SA_vertical_interval"] = { "float","10","This is used in both slope-area routines and also junction angle routines. It sets the vertical drop over which the gradient is measured and the junction angle is measured on points between a tributary within the height and the junction.","For S-A analysis the value should be greater than the vertical uncertainty of your DEM. For junction angles if you set this to a large number it will measure angles between the junction the donor junctions and the receiver junction."}; + help_map["SA_vertical_interval"] = { "float","10","This is used in both slope-area routines and also junction angle and junction bearing routines. It sets the vertical drop over which the gradient is measured and the junction angle is measured on points between a tributary within the height and the junction.","For S-A analysis the value should be greater than the vertical uncertainty of your DEM. For junction angles if you set this to a large number it will measure angles between the junction the donor junctions and the receiver junction."}; + + // This is for the angles of each + bool_default_map["print_segment_bearings_and_gradients_to_csv"] = false; + help_map["print_segment_bearings_and_gradients_to_csv"] = { "bool","false","Prints a csv with ALL the locations of channel segments (denoted by their upstream junction) giving their bearing and their gradient.","csv file contains bearings, stream orders, gradients, and other information."}; + + bool_default_map["print_segment_bearings_and_gradients_to_csv_in_basins"] = false; + help_map["print_segment_bearings_and_gradients_to_csv_in_basins"] = { "bool","false","Prints a csv with the locations of channel segments (denoted by their upstream junction) giving their bearing and their gradient within selected basins and associated statistics.","csv file contains bearings, stream orders, gradients, and other information."}; // Now for connectivity @@ -788,6 +801,12 @@ int main (int nNumberofArgs,char *argv[]) cout << "Read filename is: " << DATA_DIR+DEM_ID << endl; cout << "Write filename is: " << OUT_DIR+OUT_ID << endl; + if(this_string_map["basin_outlet_csv"]!= "NULL") + { + cout << "You gave me a basin outlet file, I am going to assume you want to select basins from outlets." << endl; + this_bool_map["get_basins_from_outlets"] = true; + } + // A little switch to find basins if you supply an outlet file if( this_bool_map["get_basins_from_outlets"] ) { @@ -2091,6 +2110,7 @@ int main (int nNumberofArgs,char *argv[]) || this_bool_map["print_MD_drainage_area_raster"] || this_bool_map["print_fill_raster"] || this_bool_map["print_stream_order_raster"] + || this_bool_map["print_flow_direction_raster"] || this_bool_map["print_channels_to_csv"] || this_bool_map["print_junction_index_raster"] || this_bool_map["print_junctions_to_csv"] @@ -2098,6 +2118,7 @@ int main (int nNumberofArgs,char *argv[]) || this_bool_map["print_channel_tips_raster"] || this_bool_map["print_chi_data_maps"] || this_bool_map["print_junction_angles_to_csv"] + || this_bool_map["print_segment_bearings_and_gradients_to_csv"] || this_bool_map["extract_single_channel"] || this_bool_map["remove_nodes_influenced_by_edge"] || this_bool_map["isolate_pixels_draining_to_fixed_channel"] @@ -2175,6 +2196,14 @@ int main (int nNumberofArgs,char *argv[]) LSDFlowInfo FlowInfo(boundary_conditions,filled_topography); cout << "Finished flow routing." << endl; + if (this_bool_map["print_flow_direction_raster"]) + { + cout << "Let me print the flow direction raster. It will use the ArcMap conventions for flow directions." << endl; + string direction_raster_name = OUT_DIR+OUT_ID+"_FlowDirection"; + LSDIndexRaster FlowDir = FlowInfo.write_FlowDirection_to_LSDIndexRaster_Arcformat(); + FlowDir.write_raster(direction_raster_name,raster_ext); + } + if (this_bool_map["remove_nodes_influenced_by_edge"]) { cout << "I am going to print a raster that has nodata for all nodes influenced by the edge or nodata" << endl; @@ -2491,6 +2520,7 @@ int main (int nNumberofArgs,char *argv[]) || this_bool_map["print_channel_tips_raster"] || this_bool_map["print_chi_data_maps"] || this_bool_map["print_junction_angles_to_csv"] + || this_bool_map["print_segment_bearings_and_gradients_to_csv"] || this_bool_map["extract_ridges"] || this_bool_map["calculate_hypsometric_integral"] || this_bool_map["clip_raster_to_basins"]) @@ -2620,7 +2650,9 @@ int main (int nNumberofArgs,char *argv[]) cout << "I am testing the junction angle code with elevations." << endl; string JAngles_csv_name = OUT_DIR+OUT_ID+"_FULL_JAngles.csv"; vector JunctionList; - JunctionNetwork.print_complete_junction_angles_to_csv(JunctionList, FlowInfo, filled_topography, FlowDistance, this_float_map["SA_vertical_interval"], JAngles_csv_name); + JunctionNetwork.print_complete_junction_angles_to_csv(JunctionList, FlowInfo, + filled_topography, FlowDistance, + this_float_map["SA_vertical_interval"], JAngles_csv_name); if ( this_bool_map["convert_csv_to_geojson"]) { @@ -2631,6 +2663,35 @@ int main (int nNumberofArgs,char *argv[]) } + if( this_bool_map["print_segment_bearings_and_gradients_to_csv"] ) + { + // Calculate flow distance + LSDRaster FlowDistance = FlowInfo.distance_from_outlet(); + + cout << endl << endl << "=======================================" << endl; + cout << "I am running the segment bearing routine. Use this to make a fun stereonet of channels and see if it matches your rocks!" << endl; + string SBearings_csv_name = OUT_DIR+OUT_ID+"_segment_bearings.csv"; + vector JunctionList; + JunctionNetwork.print_complete_segment_bearings_and_gradients_to_csv(JunctionList, FlowInfo, + filled_topography, FlowDistance, + this_float_map["SA_vertical_interval"], SBearings_csv_name); + + + if ( this_bool_map["convert_csv_to_geojson"]) + { + string gjson_name = OUT_DIR+OUT_ID+"_segment_bearings.geojson"; + LSDSpatialCSVReader thiscsv(SBearings_csv_name); + thiscsv.print_data_to_geojson(gjson_name); + + string gjson_vector_name = OUT_DIR+OUT_ID+"_segment_bearings_vector.geojson"; + JunctionNetwork.print_complete_segment_bearings_and_gradients_to_geojson(JunctionList, FlowInfo, + filled_topography, FlowDistance, + this_float_map["SA_vertical_interval"], gjson_vector_name); + } + cout << "Done with segment bearing routine." << endl << endl << endl; + } + + // Print sources if( this_bool_map["print_sources_to_csv"]) { @@ -2662,9 +2723,21 @@ int main (int nNumberofArgs,char *argv[]) if(this_bool_map["find_basins"] || this_bool_map["print_chi_data_maps"] || this_bool_map["calculate_basin_statistics"] || - this_bool_map["clip_raster_to_basins"]) + this_bool_map["clip_raster_to_basins"] || + this_bool_map["print_basin_perimeters"] || + this_bool_map["print_segment_bearings_and_gradients_to_csv_in_basins"] || + this_bool_map["print_junction_angles_to_csv_in_basins"]) { + cout << "I am now going to extract some basins for you." << endl; + cout << "The trigger for this is: " << endl; + if (this_bool_map["find_basins"]) {cout << "find basins" << endl;} + if (this_bool_map["print_chi_data_maps"]) {cout << "print_chi_data_maps" << endl;} + if (this_bool_map["calculate_basin_statistics"]) {cout << "calculate_basin_statistics" << endl;} + if (this_bool_map["clip_raster_to_basins"]) {cout << "clip_raster_to_basins" << endl;} + if (this_bool_map["print_basin_perimeters"]) {cout << "print_basin_perimeters" << endl;} + if (this_bool_map["print_segment_bearings_and_gradients_to_csv_in_basins"]) {cout << "print_segment_bearings_and_gradients_to_csv_in_basins" << endl;} + if (this_bool_map["print_junction_angles_to_csv_in_basins"]) {cout << "print_junction_angles_to_csv_in_basins" << endl;} vector BaseLevelJunctions; vector BaseLevelJunctions_Initial; @@ -2866,6 +2939,49 @@ int main (int nNumberofArgs,char *argv[]) + // Now for junction angles only within the selected basins + if(this_bool_map["print_segment_bearings_and_gradients_to_csv_in_basins"]) + { + vector basin_junctions; + for(int i = 0; i< N_BaseLevelJuncs; i++) + { + // get the upslope junctions + vector these_junctions = JunctionNetwork.get_upslope_junctions(BaseLevelJunctions[i]); + + // now append these to the master list + basin_junctions.insert( basin_junctions.end(), these_junctions.begin(), these_junctions.end() ); + } + + int N_basin_junctions = int(basin_junctions.size()); + if (N_basin_junctions == 0) + { + cout << "I am stopping here since I don't have any junctions over which to measure the segment gradients and bearings." << endl; + exit(EXIT_FAILURE); + } + + // Calculate flow distance + LSDRaster FlowDistance = FlowInfo.distance_from_outlet(); + + cout << "I am running the segment bearing routine. Use this to ." << endl; + string SBearings_csv_name = OUT_DIR+OUT_ID+"_segment_bearings.csv"; + JunctionNetwork.print_complete_segment_bearings_and_gradients_to_csv(basin_junctions, FlowInfo, + filled_topography, FlowDistance, + this_float_map["SA_vertical_interval"], SBearings_csv_name); + + + if ( this_bool_map["convert_csv_to_geojson"]) + { + string gjson_name = OUT_DIR+OUT_ID+"_segment_bearings.geojson"; + LSDSpatialCSVReader thiscsv(SBearings_csv_name); + thiscsv.print_data_to_geojson(gjson_name); + + string gjson_vector_name = OUT_DIR+OUT_ID+"_segment_bearings_vector.geojson"; + JunctionNetwork.print_complete_segment_bearings_and_gradients_to_geojson(basin_junctions, FlowInfo, + filled_topography, FlowDistance, + this_float_map["SA_vertical_interval"], SBearings_csv_name); + } + } + // Now to clip basins if(this_bool_map["clip_raster_to_basins"]) { @@ -2941,6 +3057,39 @@ int main (int nNumberofArgs,char *argv[]) } } + // Now we get the basin perimeters + if (this_bool_map["print_basin_perimeters"]) + { + cout << "I am now going to print some basin perimeter information." << endl; + cout << "The will be a different file for each basin, use the basin info csv to tell you which basin is which." << endl; + cout << "The basin perimeter files will have the basin number in the filename" << endl; + + int N_BL = int(BaseLevelJunctions.size()); + for(int bsn = 0; bsn < N_BL; bsn++) + { + cout << "Processing basin " << bsn << " of " << N_BL << endl; + LSDBasin ThisBasin(BaseLevelJunctions[bsn], FlowInfo, JunctionNetwork); + + // dont need this stuff it is in the printing routine + ThisBasin.set_Perimeter(FlowInfo); + string basin_perimeter_fname = OUT_DIR+OUT_ID+"_basin_"+itoa(bsn)+"_perimeter.csv"; + ThisBasin.print_perimeter_to_csv(FlowInfo, basin_perimeter_fname); + + + ThisBasin.clean_perimeter(FlowInfo); + basin_perimeter_fname = OUT_DIR+OUT_ID+"_basin_"+itoa(bsn)+"_CLEANperimeter.csv"; + ThisBasin.print_perimeter_to_csv(FlowInfo, basin_perimeter_fname); + + + vector Reordered_nodes = ThisBasin.order_perimeter_nodes_SMM(FlowInfo); + basin_perimeter_fname = OUT_DIR+OUT_ID+"_basin_"+itoa(bsn)+"_CLEAN_ORDEREDperimeter.csv"; + string add_column_name = "elevation(m)"; + FlowInfo.print_vector_of_nodeindices_to_csv_file_with_latlong(Reordered_nodes, basin_perimeter_fname, filled_topography, add_column_name); + + //ThisBasin.print_perimeter_hypsometry_to_csv(FlowInfo, basin_perimeter_fname, filled_topography); + } + } + //==================================================================================== // //..####...##..##..######..........######..##..##..######..#####....####....####...######..######...####...##..##. diff --git a/src/lsdtt-drivers/lsdtt-channel-extraction.cpp b/src/lsdtt-drivers/lsdtt-channel-extraction.cpp index e5b2a98..02ed64b 100644 --- a/src/lsdtt-drivers/lsdtt-channel-extraction.cpp +++ b/src/lsdtt-drivers/lsdtt-channel-extraction.cpp @@ -85,7 +85,7 @@ int main (int nNumberofArgs,char *argv[]) //start the clock clock_t begin = clock(); - string version_number = "0.8"; + string version_number = "0.9"; string citation = "http://doi.org/10.5281/zenodo.4577879"; cout << "=========================================================" << endl; diff --git a/src/lsdtt-drivers/lsdtt-chi-mapping.cpp b/src/lsdtt-drivers/lsdtt-chi-mapping.cpp index 85a6ce2..f576889 100644 --- a/src/lsdtt-drivers/lsdtt-chi-mapping.cpp +++ b/src/lsdtt-drivers/lsdtt-chi-mapping.cpp @@ -80,7 +80,7 @@ int main (int nNumberofArgs,char *argv[]) { - string version_number = "0.8"; + string version_number = "0.9"; string citation = "http://doi.org/10.5281/zenodo.4577879"; cout << "=========================================================" << endl; diff --git a/src/lsdtt-drivers/lsdtt-cosmo-tool.cpp b/src/lsdtt-drivers/lsdtt-cosmo-tool.cpp index adc02bb..94e2121 100644 --- a/src/lsdtt-drivers/lsdtt-cosmo-tool.cpp +++ b/src/lsdtt-drivers/lsdtt-cosmo-tool.cpp @@ -88,7 +88,7 @@ int main (int nNumberofArgs,char *argv[]) { - string version_number = "0.8"; + string version_number = "0.9"; string citation = "http://doi.org/10.5281/zenodo.4577879"; cout << "=========================================================" << endl; @@ -822,7 +822,7 @@ int main (int nNumberofArgs,char *argv[]) cout << "I'm afraid I can't do any cosmogenic computations without this information." << endl; cout << "Make sure you designate the cosmo_parameter_prefix in the parameter file." << endl; - cout << "I can still print a production raste or a pressure raster from here though." << endl; + cout << "I can still print a production raster or a pressure raster from here though." << endl; cout << "I am switching off the other analyses." << endl; this_bool_map["check_cosmo_basins"] = false; this_bool_map["spawn_cosmo_basins"] = false; @@ -1375,6 +1375,13 @@ int main (int nNumberofArgs,char *argv[]) cout << "I am implementing a transient column. " << endl; cout << "This requires a transient erosion file." << endl; cout << "See the help file for the format of this file." << endl; + cout << "It is a csv with three columns with headers: " << endl; + cout << "time,rate,removal" << endl; + cout << "time is the time before present of the change. So must be increasing in time" << endl; + cout << "rate erosion rate in m per year." << endl; + cout << "removal is a depth of removal in the past." << endl; + cout << "removal will be triggered if the subsequent time is the same as the past time, otherwise ignored" << endl; + cout << "The last row is the erosion rate as the background and must have a time of -9999" << endl; LSDParticleColumn PartCol; diff --git a/src/lsdtt-drivers/lsdtt-hillslope-channel-coupling.cpp b/src/lsdtt-drivers/lsdtt-hillslope-channel-coupling.cpp index 9f4a0d3..c01dcd1 100644 --- a/src/lsdtt-drivers/lsdtt-hillslope-channel-coupling.cpp +++ b/src/lsdtt-drivers/lsdtt-hillslope-channel-coupling.cpp @@ -81,7 +81,7 @@ int main (int nNumberofArgs,char *argv[]) { - string version_number = "0.8"; + string version_number = "0.9"; string citation = "http://doi.org/10.5281/zenodo.4577879"; cout << "=========================================================" << endl; @@ -218,7 +218,6 @@ int main (int nNumberofArgs,char *argv[]) int_default_map["connected_components_threshold"] = 100; help_map["connected_components_threshold"] = { "int","100","Number of connected pixels to classify as a channel.","Used in Pelletier and Wiener methods."}; - //Defining hilltops int_default_map["StreamNetworkPadding"] = 0; help_map["StreamNetworkPadding"] = { "int","0","Distance in pixels from channel network that can't be classed as a hilltop","Use if you want to remove hilltops near channels."}; @@ -335,6 +334,9 @@ int main (int nNumberofArgs,char *argv[]) bool_default_map["print_channels_to_csv"] = false; help_map["print_channels_to_csv"] = { "bool","false","Prints the channel network to a csv file.","This version produces smaller files than the raster version."}; + bool_default_map["use_extended_channel_data"] = false; + help_map["use_extended_channel_data"] = { "bool","false","If this is true you get more data columns in your channel network csv.","I will tell you what these columns are one day."}; + bool_default_map["print_junction_index_raster"] = false; help_map["print_junction_index_raster"] = { "bool","false","Prints a raster with junctions and their number.","Makes big files. It is better to use the csv version."}; @@ -897,6 +899,7 @@ int main (int nNumberofArgs,char *argv[]) //Now we have the final channel heads, so we can generate a channel network from them //LSDJunctionNetwork ChanNetwork(FinalSources, FlowInfo); + } else { @@ -953,7 +956,15 @@ int main (int nNumberofArgs,char *argv[]) { cout << "\t\tprinting channels to csv..." << endl; string channel_csv_name = OUT_DIR+OUT_ID+"_CN"; - JunctionNetwork.PrintChannelNetworkToCSV(FlowInfo, channel_csv_name); + if ( this_bool_map["use_extended_channel_data"]) + { + cout << "I am going to use the extended channel network data outputs." << endl; + JunctionNetwork.PrintChannelNetworkToCSV_WithElevation_WithDonorJunction(FlowInfo, channel_csv_name, topography_raster); + } + else + { + JunctionNetwork.PrintChannelNetworkToCSV(FlowInfo, channel_csv_name); + } // convert to geojson if that is what the user wants @@ -997,7 +1008,7 @@ int main (int nNumberofArgs,char *argv[]) if ( this_bool_map["convert_csv_to_geojson"]) { - string gjson_name = OUT_DIR+OUT_ID+"_ATsources.geojson"; + string gjson_name = OUT_DIR+OUT_ID+"_sources.geojson"; LSDSpatialCSVReader thiscsv(sources_csv_name); thiscsv.print_data_to_geojson(gjson_name); } diff --git a/src/lsdtt-drivers/lsdtt-valley-metrics.cpp b/src/lsdtt-drivers/lsdtt-valley-metrics.cpp index ac4b4be..2bbc662 100644 --- a/src/lsdtt-drivers/lsdtt-valley-metrics.cpp +++ b/src/lsdtt-drivers/lsdtt-valley-metrics.cpp @@ -74,7 +74,7 @@ int main (int nNumberofArgs,char *argv[]) //start the clock clock_t begin = clock(); - string version_number = "0.8"; + string version_number = "0.9"; string citation = "http://doi.org/10.5281/zenodo.4577879"; cout << "=========================================================" << endl; @@ -288,6 +288,9 @@ int main (int nNumberofArgs,char *argv[]) float_default_map["swath_width"] = 1000; help_map["swath_width"] = { "float","1000","Width in metres of the swath along the channel that is used to find terraces","Make this wider if you want terraces far from the channel."}; + float_default_map["swath_bin_spacing"] = 100; + help_map["swath_bin_spacing"] = { "float","1000","The spacing of the bins in metres that are used to calculate statistics along the swath.","Should be at least a few times bigger than swath_point_spacing or the pixel size of the DEM."}; + bool_default_map["print_swath_raster"] = false; help_map["print_swath_raster"] = { "bool","false","Prints the raster of the swath used to find terraces","Good for bug checking to make sure the correct river is being sampled"}; @@ -391,6 +394,7 @@ int main (int nNumberofArgs,char *argv[]) // //========================================================================= // First we need to test for filtering + cout << endl << endl << endl; cout << "====================================================" << endl; cout << " I need to check for filtering. " << endl; string filtered_prefix = DATA_DIR+DEM_ID+"_filtered"; @@ -580,56 +584,77 @@ int main (int nNumberofArgs,char *argv[]) //========================================================================= if ( this_bool_map["extract_single_channel"]) { + cout << endl << endl << endl; + cout << "I am going to extract a channel before I get to the floodplain analysis." << endl; + cout << "If you chose this option you can use it to search along terraces. " << endl; + cout << endl << endl << endl; + string filename = OUT_DIR+OUT_ID+"_swath_channel_nodes.csv"; ifstream test_single_chan(filename.c_str()); if (test_single_chan) { - cout << "I found a channel file already. I'll just load that" << endl; + cout << "I found a channel file already. I'll just load that later when you need it." << endl; + cout << "If you updated information about the single channel" << endl; + cout << "YOU NEED TO DELETE THE FILE and run this again. " << endl; } else { - cout << "I am going to extract a channel. " << endl; - cout << "To do this I need to get the flow info object, which will take a little while." << endl; - cout << "Also I am assuming your first point is the upstream point" << endl; - cout << "and the second point is the downstream point." << endl; - cout << "If you don't have a downstream point" << endl; - cout << "I just use the first point as the source and follow it to" << endl; - cout << "the edge of the DEM." << endl; + int start_NI,downstream_NI; + bool end_node_in_csv = false; + vector NI_of_points; - cout << "I am reading points from the file: "+ this_string_map["channel_source_fname"] << endl; - LSDSpatialCSVReader CSVFile( RI, (DATA_DIR+this_string_map["channel_source_fname"]) ); + if (this_string_map["channel_source_fname"] != "NULL") + { + LSDRasterInfo this_RI(filled_topography); + cout << "I am reading points from the file: "+ this_string_map["channel_source_fname"] << endl; + LSDSpatialCSVReader CSVFile( this_RI, (DATA_DIR+this_string_map["channel_source_fname"]) ); - // get the x and y locations of the points - vector UTME; - vector UTMN; - vector latitudes = CSVFile.get_latitude(); - vector longitudes = CSVFile.get_longitude(); - CSVFile.get_x_and_y_from_latlong(UTME,UTMN); + // get the x and y locations of the points + vector UTME; + vector UTMN; + vector latitudes = CSVFile.get_latitude(); + vector longitudes = CSVFile.get_longitude(); - int n_nodes = int(UTME.size()); - bool end_node_in_csv = false; - if (n_nodes > 1) - { - end_node_in_csv = true; - cout << "You have given me a channel source file with a second node." << endl; - cout << "I am assuming this is an end node. " << endl; - cout << "I will search for this node." << endl; + CSVFile.get_x_and_y_from_latlong(UTME,UTMN); + + NI_of_points = CSVFile.get_nodeindices_from_lat_long(FlowInfo); + + int n_nodes = NI_of_points.size(); + if (n_nodes > 1) + { + end_node_in_csv = true; + cout << "You have given me a channel source file with a second node." << endl; + cout << "I am assuming this is an end node. " << endl; + cout << "I will search for this node." << endl; + } + // find the nearest channel to the starting node. This is necessary so that the start node is definitely in the centreline + // this is a bit of a pain because of the different flow infos... + int snapping_SO = 1; + int snapped_NI = ChanNetwork.find_nearest_downslope_channel(UTME[0], UTMN[0], snapping_SO, FlowInfo); + float snapped_UTME, snapped_UTMN; + FlowInfo.get_x_and_y_from_current_node(snapped_NI, snapped_UTME, snapped_UTMN); + start_NI = FlowInfo.get_node_index_of_coordinate_point(snapped_UTME, snapped_UTMN); + cout << "start ni " << start_NI << endl; + cout << snapped_UTME << " " << snapped_UTMN << endl; } + else + { + cout << "I didn't find a channel source, so I am going to route from the longest channel. " << endl; + // get the maximum elevation + int max_row,max_col; + float max_dist = DistanceFromOutlet.max_elevation(max_row,max_col); - vector NI_of_points = CSVFile.get_nodeindices_from_lat_long(FlowInfo); + start_NI = FlowInfo.get_NodeIndex_from_row_col(max_row,max_col); + } - // flow path in NI - vector final_channel_list; - vector flow_path = FlowInfo.get_flow_path(UTME[0],UTMN[0]); - int downstream_NI; + vector flow_path = FlowInfo.get_flow_path(start_NI); + cout << "n nodes " << int(flow_path.size()) << endl; + // loop through the flow path and find the node nearest to the end node (if one exists) if(end_node_in_csv) { // snap the lower point to a channel cout << "The search radius is: " << this_int_map["search_radius_nodes"] << endl; - downstream_NI = ChanNetwork.get_nodeindex_of_nearest_channel_for_specified_coordinates(UTME[1],UTMN[1], - this_int_map["search_radius_nodes"], - this_int_map["threshold_stream_order_for_snapping"], - FlowInfo); + downstream_NI = FlowInfo.find_nearest_node_in_vector(NI_of_points[1], flow_path); } else { @@ -640,6 +665,7 @@ int main (int nNumberofArgs,char *argv[]) // travel down the flow path until you either reach the outlet node or // get to the end + vector final_channel_list; for(int node = 0; node < int(flow_path.size()); node ++) { final_channel_list.push_back(flow_path[node]); @@ -676,20 +702,21 @@ int main (int nNumberofArgs,char *argv[]) } } -//========================================================================= -// .__ __ .__ .__ -// _____ __ __| |_/ |_|__|_____ | | ____ -// / \| | \ |\ __\ \____ \| | _/ __ \ -// | Y Y \ | / |_| | | | |_> > |_\ ___/ -// |__|_| /____/|____/__| |__| __/|____/\___ > -// \/ |__| \/ -// .__ .__ -// ____ | |__ _____ ____ ____ ____ | | ______ -// _/ ___\| | \\__ \ / \ / \_/ __ \| | / ___/ -// \ \___| Y \/ __ \| | \ | \ ___/| |__\___ \ -// \___ >___| (____ /___| /___| /\___ >____/____ > -// \/ \/ \/ \/ \/ \/ \/ -//========================================================================= + //========================================================================= + // + //.##...##..##..##..##......######..######..#####...##......######. + //.###.###..##..##..##........##......##....##..##..##......##..... + //.##.#.##..##..##..##........##......##....#####...##......####... + //.##...##..##..##..##........##......##....##......##......##..... + //.##...##...####...######....##....######..##......######..######. + //................................................................. + //..####...##..##...####...##..##..##..##..######..##.......####.. + //.##..##..##..##..##..##..###.##..###.##..##......##......##..... + //.##......######..######..##.###..##.###..####....##.......####.. + //.##..##..##..##..##..##..##..##..##..##..##......##..........##. + //..####...##..##..##..##..##..##..##..##..######..######...####.. + // + //========================================================================= if(this_bool_map["extract_multiple_channels"]) { cout << "I am going to extract multiple channels. " << endl; @@ -758,21 +785,21 @@ int main (int nNumberofArgs,char *argv[]) } } -//-------------------------------------------------------------------------// -// __ ______ .__ __. _______ _______ _______.___________. -// | | / __ \ | \ | | / _____|| ____| / | | -// | | | | | | | \| | | | __ | |__ | (----`---| |----` -// | | | | | | | . ` | | | |_ | | __| \ \ | | -// | `----.| `--' | | |\ | | |__| | | |____.----) | | | -// |_______| \______/ |__| \__| \______| |_______|_______/ |__| - -// ______ __ __ ___ .__ __. .__ __. _______ __ -// / || | | | / \ | \ | | | \ | | | ____|| | -// | ,----'| |__| | / ^ \ | \| | | \| | | |__ | | -// | | | __ | / /_\ \ | . ` | | . ` | | __| | | -// | `----.| | | | / _____ \ | |\ | | |\ | | |____ | `----. -// \______||__| |__| /__/ \__\ |__| \__| |__| \__| |_______||_______| -//-------------------------------------------------------------------------// + //-------------------------------------------------------------------------// + // __ ______ .__ __. _______ _______ _______.___________. + // | | / __ \ | \ | | / _____|| ____| / | | + // | | | | | | | \| | | | __ | |__ | (----`---| |----` + // | | | | | | | . ` | | | |_ | | __| \ \ | | + // | `----.| `--' | | |\ | | |__| | | |____.----) | | | + // |_______| \______/ |__| \__| \______| |_______|_______/ |__| + + // ______ __ __ ___ .__ __. .__ __. _______ __ + // / || | | | / \ | \ | | | \ | | | ____|| | + // | ,----'| |__| | / ^ \ | \| | | \| | | |__ | | + // | | | __ | / /_\ \ | . ` | | . ` | | __| | | + // | `----.| | | | / _____ \ | |\ | | |\ | | |____ | `----. + // \______||__| |__| /__/ \__\ |__| \__| |__| \__| |_______||_______| + //-------------------------------------------------------------------------// if(this_bool_map["extract_longest_channel"]) { cout << "I am going to extract the longest channel in the DEM." << endl; @@ -786,50 +813,62 @@ int main (int nNumberofArgs,char *argv[]) cout << " BL Junction: " << baselevel_junction << endl; // get the source of this baselevel node - LSDIndexChannel LongestChan = ChanNetwork.generate_longest_index_channel_in_basin(baselevel_junction, FlowInfo, - DistanceFromOutlet); - - // write the whole channel to a CSV - string out_filename = OUT_ID+"_longest_channel"; - LongestChan.write_channel_to_csv(OUT_DIR, out_filename, FlowInfo, DistanceFromOutlet, filled_topography); + vector this_long_channel = ChanNetwork.generate_longest_nodeindex_channel_in_basin(baselevel_junction, FlowInfo, + DistanceFromOutlet); + + cout << "writing to a csv file" << endl; + + string out_filename = OUT_ID+"_swath_channel_longest"; + FlowInfo.print_vector_of_nodeindices_to_csv_file_with_latlong(this_long_channel, + OUT_DIR, out_filename,filled_topography, DistanceFromOutlet, + DA_d8); + + + if ( this_bool_map["convert_csv_to_geojson"]) + { + string in_name = OUT_DIR+OUT_ID+"_swath_channel_longest_nodes.csv"; + cout << "Printing the geojson of the swath channel." << endl; + string gjson_name = OUT_DIR+OUT_ID+"_swath_channel_longest.geojson"; + LSDSpatialCSVReader thiscsv(in_name); + thiscsv.print_data_to_geojson(gjson_name); + } - // just write the starting and ending nodes to a csv - ofstream source_out; - string source_out_fname = OUT_DIR+OUT_ID+"_longest_channel_source_points.csv"; - source_out.open(source_out_fname.c_str()); - source_out << "id,latitude,longitude" << endl; + // just write the starting and ending nodes to a csv + //ofstream source_out; + //string source_out_fname = OUT_DIR+OUT_ID+"_longest_channel_source_points.csv"; + //source_out.open(source_out_fname.c_str()); + //source_out << "id,latitude,longitude" << endl; // get lat and longitude of source and outlet nodes - LSDCoordinateConverterLLandUTM Converter; - vector X_coords, Y_coords; - double lat, lon; - LongestChan.get_coordinates_of_channel_nodes(X_coords, Y_coords, FlowInfo); + //LSDCoordinateConverterLLandUTM Converter; + //vector X_coords, Y_coords; + //double lat, lon; + //LongestChan.get_coordinates_of_channel_nodes(X_coords, Y_coords, FlowInfo); // source node - FlowInfo.get_lat_and_long_locations(X_coords.front(), Y_coords.front(), lat, lon, Converter); - source_out << "1," << lat << "," << lon << endl; + //FlowInfo.get_lat_and_long_locations(X_coords.front(), Y_coords.front(), lat, lon, Converter); + //source_out << "1," << lat << "," << lon << endl; // outlet node - FlowInfo.get_lat_and_long_locations(X_coords.back(), Y_coords.back(), lat, lon, Converter); - source_out << "2," << lat << "," << lon << endl; - - source_out.close(); + //FlowInfo.get_lat_and_long_locations(X_coords.back(), Y_coords.back(), lat, lon, Converter); + //source_out << "2," << lat << "," << lon << endl; + //source_out.close(); } //========================================================================= // - // o 'O o o Oo o - // O o O O o O O o - // o O o o O o o - // o o O O oOooOoOo O - // O O' .oOoO' o o .oOo. O o o O 'OoOo. .oOoO' o O o .oOo O .oOo - // `o o O o O O OooO' o O O o o O O o O o O `Ooo. o `Ooo. - // `o O o O o o O O o o O O o o O o O o O O O - // `o' `OoO'o Oo Oo `OoO' `OoOO O. O o O `OoO'o Oo `OoOO `OoO' o' `OoO' - // o o - // OoO' OoO' + //.######..##..##..######..#####....####....####...######. + //.##.......####.....##....##..##..##..##..##..##....##... + //.####......##......##....#####...######..##........##... + //.##.......####.....##....##..##..##..##..##..##....##... + //.######..##..##....##....##..##..##..##...####.....##... + //........................................................ + //.######..##.......####....####...#####...#####...##.......####...######..##..##. + //.##......##......##..##..##..##..##..##..##..##..##......##..##....##....###.##. + //.####....##......##..##..##..##..##..##..#####...##......######....##....##.###. + //.##......##......##..##..##..##..##..##..##......##......##..##....##....##..##. + //.##......######...####....####...#####...##......######..##..##..######..##..##. // //========================================================================= - cout << endl << endl << "=========================================================" << endl; cout << "In this analysis (lsdtt-valley-metrics) I automatically extract floodplains" << endl; cout << "I am going to extract several rasters for you." << endl; @@ -891,8 +930,12 @@ int main (int nNumberofArgs,char *argv[]) if(this_bool_map["use_absolute_thresholds"]) { + cout << endl << endl << endl << "--------------------------------------------" << endl; + cout << "I am going to calculate the floodplain pixels using absolute thresholds for slope and relief." << endl; relief_threshold = this_int_map["relief_threshold"]; slope_threshold = this_float_map["slope_threshold"]; + cout << "Thresholds are R: " << relief_threshold << " S: " << slope_threshold << endl; + cout << "--------------------------------------------" << endl << endl << endl << endl; } else { @@ -962,12 +1005,12 @@ int main (int nNumberofArgs,char *argv[]) cout << "=========================================================" << endl << endl; //============================================================================================= - // _ _ _ _ _ _ _ - // __ __ __ _ | | | | ___ | || | o O O __ ___ _ _ | |_ _ _ ___ | | (_) _ _ ___ o O O - // \ V / / _` | | | | | / -_) \_, | o / _| / -_) | ' \ | _| | '_| / -_) | | | | | ' \ / -_) o - // _\_/_ \__,_| _|_|_ _|_|_ \___| _|__/ TS__[O] \__|_ \___| |_||_| _\__| _|_|_ \___| _|_|_ _|_|_ |_||_| \___| TS__[O] - // _|"""""|_|"""""|_|"""""|_|"""""|_|"""""|_| """"| {======|_|"""""|_|"""""|_|"""""|_|"""""|_|"""""|_|"""""|_|"""""|_|"""""|_|"""""|_|"""""| {======| - // "`-0-0-'"`-0-0-'"`-0-0-'"`-0-0-'"`-0-0-'"`-0-0-'./o--000'"`-0-0-'"`-0-0-'"`-0-0-'"`-0-0-'"`-0-0-'"`-0-0-'"`-0-0-'"`-0-0-'"`-0-0-'"`-0-0-'./o--000' + // + //.##..##...####...##......##......######..##..##...........####...######..##..##..######..#####...######..##......######..##..##..######. + //.##..##..##..##..##......##......##.......####...........##..##..##......###.##....##....##..##..##......##........##....###.##..##..... + //.##..##..######..##......##......####......##............##......####....##.###....##....#####...####....##........##....##.###..####... + //..####...##..##..##......##......##........##............##..##..##......##..##....##....##..##..##......##........##....##..##..##..... + //...##....##..##..######..######..######....##.............####...######..##..##....##....##..##..######..######..######..##..##..######. // //============================================================================================= if(this_bool_map["get_valley_centreline"]) @@ -984,6 +1027,9 @@ int main (int nNumberofArgs,char *argv[]) { cout << "I have not found a centreline file so I need to generate one." << endl; + cout << "This differs from a single channel in that it goes down the centre of the valley" << endl; + cout << "and not down the channel pathway. " << endl; + cout << "See Clubb et al 2022 ESURF https://doi.org/10.5194/esurf-10-437-2022 for details." << endl; cout << "First what I will do is take the floodplain and punch a nodata raster" << endl; LSDRaster Punched = Floodplain.get_NoData_FloodplainRaster(filled_topography); string punched_name = OUT_DIR+OUT_ID+"_Punched"; @@ -1030,24 +1076,8 @@ int main (int nNumberofArgs,char *argv[]) river_path = river_path.MapAlgebra_subtract(trough); LSDRaster carved_topography = river_path.Breaching_Lindsay2016(); river_path = carved_topography.fill(this_float_map["min_slope_for_fill"]); - - // We keep track of the trough additions so at the end we add these back on except for the last one - // this means that the depth of the middle of the trough should be set by the trough_scaling_factor - // and not related to the number of loops - // if (i < n_loops-2) - // { - // add_trough = add_trough.MapAlgebra_add(trough); - // } } - //river_path = river_path.MapAlgebra_add(add_trough); - //river_path.AdjustElevation(this_float_map["river_depth"]); - - - // Now merge the rasters - //cout << "Combining rasters." << endl; - //LSDRaster topo_copy = topography_raster; - //topo_copy.OverwriteRaster(river_path); string imposedriver_R_name = OUT_DIR+OUT_ID+"_RiverTrough"; //topo_copy.write_raster(imposedriver_R_name,DEM_extension); @@ -1065,78 +1095,118 @@ int main (int nNumberofArgs,char *argv[]) hs_raster.write_raster(hs_fname,DEM_extension); } - // Now we need the flow path + // Now we need the flow path this is a new flow info object that includes + // the trough topography. cout << "\t Flow routing for centreline. Note this is memory intensive. If your DEM is very large you may get a segmentation fault here..." << endl; // get a flow info object LSDFlowInfo FI(boundary_conditions,river_path); cout << "Finished flow routing." << endl; // get the junction network - int trough_threshold = 1; - LSDIndexRaster FA = FlowInfo.write_NContributingNodes_to_LSDIndexRaster(); + //int trough_threshold = 1; + //LSDIndexRaster FA = FlowInfo.write_NContributingNodes_to_LSDIndexRaster(); + + cout << "I need to calculate the flow distance now." << endl; + LSDRaster FD = FI.distance_from_outlet(); + LSDRaster DA = FI.write_DrainageArea_to_LSDRaster(); + cout << "Got the dist from outlet and DA" << endl; + // now get the starting node for the flow path from your csv file - // this is assuming you want a single channel int start_NI, downstream_NI; bool end_node_in_csv = false; vector NI_of_points; - if (this_bool_map["extract_single_channel"] || this_bool_map["extract_longest_channel"]) + vector flow_path; + if (this_bool_map["extract_single_channel"]) { - cout << "I am reading points from the file: "+ this_string_map["channel_source_fname"] << endl; - LSDSpatialCSVReader CSVFile( RI, (DATA_DIR+this_string_map["channel_source_fname"]) ); + if (this_string_map["channel_source_fname"] != "NULL") + { + cout << "I am reading points from the file: "+ this_string_map["channel_source_fname"] << endl; + LSDSpatialCSVReader CSVFile( RI, (DATA_DIR+this_string_map["channel_source_fname"]) ); + + cout << "Because I am reading a channel source file, I am switching off the longest channel." << endl; + this_bool_map["extract_longest_channel"] = false; + + // get the x and y locations of the points + vector UTME; + vector UTMN; + vector latitudes = CSVFile.get_latitude(); + vector longitudes = CSVFile.get_longitude(); + CSVFile.get_x_and_y_from_latlong(UTME,UTMN); + + NI_of_points = CSVFile.get_nodeindices_from_lat_long(FI); + + int n_nodes = NI_of_points.size(); + if (n_nodes > 1) + { + end_node_in_csv = true; + cout << "You have given me a channel source file with a second node." << endl; + cout << "I am assuming this is an end node. " << endl; + cout << "I will search for this node." << endl; + } + // find the nearest channel to the starting node. This is necessary so that the start node is definitely in the centreline + // this is a bit of a pain because of the different flow infos... + int snapping_SO = 1; + int snapped_NI = ChanNetwork.find_nearest_downslope_channel(UTME[0], UTMN[0], snapping_SO, FlowInfo); + float snapped_UTME, snapped_UTMN; + FlowInfo.get_x_and_y_from_current_node(snapped_NI, snapped_UTME, snapped_UTMN); + start_NI = FI.get_node_index_of_coordinate_point(snapped_UTME, snapped_UTMN); + cout << "start ni " << start_NI << endl; + cout << "easting: "<< snapped_UTME << " northing: " << snapped_UTMN << endl; + } + else if (this_bool_map["extract_longest_channel"]) + { + cout << "You have chosen to extract the longest channel in the dataset. " << endl; + int long_row, long_col; + float longest_dist = FD.max_elevation(long_row,long_col); - // get the x and y locations of the points - vector UTME; - vector UTMN; - vector latitudes = CSVFile.get_latitude(); - vector longitudes = CSVFile.get_longitude(); - CSVFile.get_x_and_y_from_latlong(UTME,UTMN); + start_NI = FI.get_NodeIndex_from_row_col(long_row,long_col); + downstream_NI = -1; + } + else + { + cout << "I didn't find a channel source, so I am going to route from the longest channel in the river path. " << endl; + // get the maximum elevation + int long_row, long_col; + float longest_dist = FD.max_elevation(long_row,long_col); - NI_of_points = CSVFile.get_nodeindices_from_lat_long(FI); + start_NI = FI.get_NodeIndex_from_row_col(long_row,long_col); + downstream_NI = -1; + } - int n_nodes = NI_of_points.size(); - if (n_nodes > 1) + flow_path = FI.get_flow_path(start_NI); + if( int(flow_path.size()) < 10) { - end_node_in_csv = true; - cout << "You have given me a channel source file with a second node." << endl; - cout << "I am assuming this is an end node. " << endl; - cout << "I will search for this node." << endl; + cout << endl << endl << "======================================" << endl; + cout << "n nodes in the intitial flow path for the centreline is: " << int(flow_path.size()) << endl; + cout << "This is small and the routing through hte valley centre may have not worked" << endl; + cout << "You can try increasing the trough_scaling_factor to get a deeper channel trough" << endl; + cout << "Default is 0.1 but you can raise this to 0.25 or even 0.5 in extreme cases" << endl; + cout << "======================================" << endl << endl; } - // find the nearest channel to the starting node. This is necessary so that the start node is definitely in the centreline - // this is a bit of a pain because of the different flow infos... - int snapping_SO = 1; - int snapped_NI = ChanNetwork.find_nearest_downslope_channel(UTME[0], UTMN[0], snapping_SO, FlowInfo); - float snapped_UTME, snapped_UTMN; - FlowInfo.get_x_and_y_from_current_node(snapped_NI, snapped_UTME, snapped_UTMN); - start_NI = FI.get_node_index_of_coordinate_point(snapped_UTME, snapped_UTMN); - cout << "start ni " << start_NI << endl; - cout << snapped_UTME << " " << snapped_UTMN << endl; - } - else - { - cout << "I didn't find a channel source, so I am going to route from the highest elevation in the river path. " << endl; - // get the maximum elevation - int max_row,max_col; - float max_elev = river_path.max_elevation(max_row,max_col); + + // loop through the flow path and find the node nearest to the end node (if one exists) + if(end_node_in_csv) + { + // snap the lower point to a channel + cout << "The search radius is: " << this_int_map["search_radius_nodes"] << endl; + downstream_NI = FI.find_nearest_node_in_vector(NI_of_points[1], flow_path); + cout << "Found nearest downstream node. it is: " << downstream_NI << endl; + if (downstream_NI = -9999) + { + downstream_NI = -1; + } + } + else + { + // There is no end node in the csv so we set the end node to a pixel that does not exist + cout << "You didn't specify a downstream node, I'll follow down to the outlet" << endl; + downstream_NI = -1; + } + } // end extract channel loop + - start_NI = FI.get_NodeIndex_from_row_col(max_row,max_col); - } - vector flow_path = FI.get_flow_path(start_NI); - cout << "n nodes " << int(flow_path.size()) << endl; - // loop through the flow path and find the node nearest to the end node (if one exists) - if(end_node_in_csv) - { - // snap the lower point to a channel - cout << "The search radius is: " << this_int_map["search_radius_nodes"] << endl; - downstream_NI = FI.find_nearest_node_in_vector(NI_of_points[1], flow_path); - } - else - { - // There is no end node in the csv so we set the end node to a pixel that does not exist - cout << "You didn't specify a downstream node, I'll follow down to the outlet" << endl; - downstream_NI = -1; - } // travel down the flow path until you either reach the outlet node or // get to the end @@ -1159,13 +1229,9 @@ int main (int nNumberofArgs,char *argv[]) } } - cout << "I need to calculate the flow distance now." << endl; - LSDRaster FD = FI.distance_from_outlet(); - LSDRaster DA = FI.write_DrainageArea_to_LSDRaster(); - cout << "Got the dist from outlet and DA" << endl; string fname = OUT_ID+"_valley_centreline"; - FI.print_vector_of_nodeindices_to_csv_file_with_latlong(flow_path,OUT_DIR, fname, + FI.print_vector_of_nodeindices_to_csv_file_with_latlong(final_channel_list,OUT_DIR, fname, river_path, FD, DA); if ( this_bool_map["convert_csv_to_geojson"]) @@ -1191,7 +1257,7 @@ int main (int nNumberofArgs,char *argv[]) if (this_bool_map["extract_multiple_channels"]) { cout << "You chose to extract multiple channels earlier in this routine." << endl; - cout << "I am assuming you want to use those channel so I am overriding" << endl; + cout << "I am assuming you want to use those channels so I am overriding" << endl; cout << "Your valley_points_csv filename." << endl; LSDSpatialCSVReader CSVFile( RI, (DATA_DIR+this_string_map["channel_source_fname"]) ); @@ -1226,7 +1292,7 @@ int main (int nNumberofArgs,char *argv[]) } } - if (this_bool_map["extract_all_channels"]) + else if (this_bool_map["extract_all_channels"]) { cout << endl << endl << "=========================================================" << endl; cout << "I am going to calculate width for all channels in the DEM. This might be slow. " << endl; @@ -1264,7 +1330,7 @@ int main (int nNumberofArgs,char *argv[]) cout << "You chose to extract the longest channel earlier in this routine." << endl; cout << "I am assuming you want to use that channel so I am overriding" << endl; cout << "Your valley_points_csv filename." << endl; - this_string_map["valley_points_csv"] = OUT_DIR+OUT_ID+"_longest_channel_index_chan.csv"; + this_string_map["valley_points_csv"] = OUT_DIR+OUT_ID+"_swath_channel_longest_nodes.csv"; } if (this_bool_map["get_valley_centreline"]) { @@ -1306,19 +1372,44 @@ int main (int nNumberofArgs,char *argv[]) //========================================================================= if (this_bool_map["get_terraces"]) { + cout << endl << endl << endl << "===============================================" << endl; + cout << "Let me compute some terraces for you" << endl; + // Read in the channel if (this_bool_map["extract_single_channel"]) { cout << "You chose to extract a single channel earier in this routine." << endl; cout << "I am assuming you want to use that channel so I am overriding" << endl; cout << "Your valley_points_csv filename." << endl; + cout << "I am also overriding instructions to extract the longest channel," << endl; + cout << "valley centreline, and a uder defined valley points csv" << endl; + cout << "That is, extract single channel takes precedence over other methods." << endl; + this_bool_map["extract_longest_channel"] = false; + this_bool_map["get_valley_centreline"] = false; + this_bool_map["use_valley_csv_for_terraces"] = false; this_string_map["valley_points_csv"] = OUT_DIR+OUT_ID+"_swath_channel_nodes.csv"; } + else if (this_bool_map["extract_longest_channel"]) + { + cout << "You chose to extract the longest channel earier in this routine." << endl; + cout << "I am assuming you want to use that channel so I am overriding" << endl; + cout << "options for the valley centreline and a valley csv." << endl; + this_bool_map["get_valley_centreline"] = false; + this_bool_map["use_valley_csv_for_terraces"] = false; + this_string_map["valley_points_csv"] = OUT_DIR+OUT_ID+"_swath_channel_longest_nodes.csv"; + } + else if (this_bool_map["get_valley_centreline"]) + { + cout << "I am going to use the file designated by the valley centreline" << endl; + cout << "that you calculated previously to calculate relief. " << endl; + cout << "I am overriding instructions to use a valley csv" << endl; + this_bool_map["use_valley_csv_for_terraces"] = false; + this_string_map["valley_points_csv"] = OUT_DIR+OUT_ID+"__valley_centreline_nodes.csv"; + } else if (this_bool_map["use_valley_csv_for_terraces"]) { - cout << "I am going to use the file designated by valley_points_csv" << endl; - cout << "to calculate relief. " << endl; - cout << "This could just be a channel network." << endl; + cout << "I am going to use the file designated by the user suppied" << endl; + cout << "valley_points_csv, which is: " << this_string_map["valley_points_csv"] << endl; } else { @@ -1401,11 +1492,12 @@ int main (int nNumberofArgs,char *argv[]) } else { + cout << endl << endl << "I am calculating the swath raster. This is computationally expensive and might take a while." << endl; vector< LSDRaster > swath_rasters = filled_topography.find_nearest_point_from_list_of_points_with_relief(spaced_eastings, spaced_northings, spaced_distances, this_float_map["swath_width"]); vector channel_node_list = CSVFile.data_column_to_int("id"); // now get the channel relief along the swath - LSDRaster swath_channel_relief = swath_rasters[3]; + swath_channel_relief = swath_rasters[3]; string relief_fname = "_swath_channel_relief"; swath_channel_relief.write_raster(OUT_DIR+OUT_ID+relief_fname, DEM_extension); @@ -1416,7 +1508,7 @@ int main (int nNumberofArgs,char *argv[]) // remove any stupid slope values LSDRaster Slope_new = Slope.mask_to_nodata_using_threshold(mask_threshold, below); - bool debug_terrace = true; + bool debug_terrace = false; if(debug_terrace) { LSDRasterInfo SC_RI(swath_channel_relief); @@ -1427,8 +1519,13 @@ int main (int nNumberofArgs,char *argv[]) } // get the terrace pixels + cout << endl << "-----------------" << endl; cout << "Printing the terraces tagged by their terrace number" << endl; + cout << "The minimum terrace height is the relief threshold plus user defined minimum_terrace_height, the latter of which is the height above the floodplain." << endl; + cout << "The floodplains are under: " << relief_threshold << " and the height_above_this is: " << this_float_map["minimum_terrace_height"] << endl; + float terrace_relief_threshold = relief_threshold+this_float_map["minimum_terrace_height"]; + cout << "So top of terraces will be at: " << terrace_relief_threshold << " m and bottom at: " << this_float_map["minimum_terrace_height"] << endl; LSDTerrace Terraces(swath_channel_relief, Slope_new, ChanNetwork, FlowInfo, terrace_relief_threshold, slope_threshold, this_int_map["minimum_patch_size"], this_int_map["threshold_SO"], this_float_map["minimum_terrace_height"]); LSDIndexRaster ConnectedComponents = Terraces.print_ConnectedComponents_to_Raster(); string CC_ext = "_terraces"; @@ -1441,14 +1538,14 @@ int main (int nNumberofArgs,char *argv[]) string R_ext = "_terrace_relief"; TerraceRelief.write_raster((OUT_DIR+OUT_ID+R_ext), DEM_extension); - //cout << "\t Testing connected components" << endl; - //vector > CC_vector = TestSwath.get_connected_components_along_swath(ConnectedComponents, RasterTemplate, this_int_map["NormaliseToBaseline"]); + string swath_data_prefix = OUT_DIR+OUT_ID+"_categorised"; + bool we_dont_need_to_do_this_again = false; + cout << "Making the final terrace swath" << endl; + TerraceRelief.make_swaths_from_categorised(spaced_eastings, spaced_northings,spaced_distances, this_float_map["swath_width"], + this_float_map["swath_bin_spacing"], swath_data_prefix, we_dont_need_to_do_this_again, + ConnectedComponents); + - // print the terrace information to a csv - //string csv_fname = "_terrace_info.csv"; - //string full_csv_name = DATA_DIR+DEM_ID+csv_fname; - //cout << "The full csv filename is: " << full_csv_name << endl; - //Terraces.print_TerraceInfo_to_csv(full_csv_name, filled_topography, swath_channel_relief, FlowInfo, TestSwath); } // Done, check how long it took