terrain evolution simulator
1996 (updated : 2008.07.05)
The software is designed to model the evolution of a square piece of terrain. It uses a fractal pattern to generate the terrain, then uses some basic geographic concepts to model erosion due to rainfall.
If you're a Geography or Geology student, perhaps you'll find this useful.
I wrote this in the early 90's for the DOS platform. It requires the fastgraf code library (gods only know where you'll find a compatible version after all these years).
/* THE TERRAIN EVOLUTION SIMULATOR v 1.21B by Michael Werneburg This program develops a model terrain surface through the use of a fractal (self-similar at all scales) framework. It also then models the pattern of drainage which would develop on that surface, and simulates the effects of the flowing water on the terrain surface. This results in an evolution of the terrain surface. This evolution is enacted through the use of two erosional procedures, one of which accounts for the combined effects of overland flow and rilling, the other which accounts for the effects of channel erosion. The first of these is a function of surface slope and slope length, while the second is a function of slope and runoff. The equations on which the coded implementations are based may be found in the accompanying text. The modelling process involves an artificial terrain 127 x 127 elements in size. For reasons of accuracy with regards to stream modelling (see text), the ground dimensions of these positions has been set to 10 metres by 10 metres, resulting in an array position ground area of 100m˝. If the ground scale is to represent real values, then the elveation values used in the model must also represent measurable real values. The scale of the vertical range was determined by calculating the average slope found by establishing a sound range of values which could be consistently modelled using an integer data structure to hold elevation values. This was accomplished as follows. The integer data type meant that elevation values could range from -32535 to 32535. It was deemed that a greater range in elevation was desirable as a greater vertical range would result in a finer measurement of real elevation values: this finer measurement was imperative as the rates of denudation were in the order of a few millimetres per year. Since the model used to create the terrain elevation values was stochastic, the range of values thus created could not be known until a number of runs offered a good estimate (although this could be roughly controlled by setting the vertical range of possible values for each random generation. The ability to control the results of the vertical scale generation was used with the goal of achieving an average local slope of 4o throughout the terrain model. The value of 4o was judged as a reasonable approximation of global average slope. Through the fiddling of the vertical range generation to result in an average slope of 4o, a real elevation value of 1 mm could be assigned to each vertical unit. This vertical value, described hereafter as the terrain evolution simulator (TES) scale, means that a range of real elevations of 50 metres is possible. The process by which the evolution model was implemented is discussed below, at the head of the declartions of those functions which implement each step. */ ///////////////////////////////////////////////////////////////////////////// // The following inclusions compile standard C library headers. #include//#include #include #include #include #include #include #include "tes.h" #define toupper(c) ((c)+'A'-'a') #define SANITY_TEST ((Y+local_DY)<127)&&((X+local_DX)<127)&&((X+local_DX)>=0)&&((Y+local_DY)>=0) #define SANITY ((Y+(D_xy*dy))<127)&&((X+(D_xy*dx))<127)&&((X+(D_xy*dx))>=0)&&((Y+(D_xy*dy))>=0) #define SANITY_CHECK ((Y+DY)<127)&&((X+DX)<127)&&((X+DX)>=0)&&((Y+DY)>=0) #define XY (X*127)+Y // The following functions are described in the code below where they are defined. void intro(), setup_graphics(), terrain_gen(), terrain_rec(), erode(), legend(int green_or_grey), ldraw(int x, int y, int col), draw(int x, int y, int col); int vari8(unsigned char scale,unsigned char midX,unsigned char midY), initialize(), values_gen(unsigned char isit_first), get_terrain(), smooth_4stats(); unsigned char huge *f_model; // The use of these data structures is outlined below unsigned char huge *river_a; int huge *terrain; int huge *delta_z; int huge *entrain; int huge *lakes; unsigned int X, // The X position within the arrays. Y; // The Y position within the arrays. unsigned char ask=1, // holder for command-line driven request for no DEM output. rivers=1; int highcount = -32000, // This is used to determine the highest terrain elevation. lowcount = 32000, // This is used to determine the lowest terrain elevation. old_vmode, lines, A, // These four variables are used as limits for windowing procedures which are performed during B, // smoothing, slope calculation, aspect generation, and watershed calculation. C, D, tilt=0, // the terrain will be titled.... tiltz=0, maxx, // width of screen maxy; // height of screen long avecount; // This is a holder for average slope and elevation values. char D_type = 'O', // holder for data type recorded to and read from Idrisi DEMs on disk noise, // holder for user key hits, so they're out of the way output_h[30], // string for output of highest elevation, slope values output_a[30], // string for output of average elevation, slope values output_l[30]; // string for output of lowest elevation, slope values float hold_v; // holder for output values. ////////////////////////////////////////MAIN()//FUNCTION/////////////////// /* This function calls the others in order of operation. It first * initializes the arrays f_model[] and terrain[] through the initialize() * function. It then initiates the graphics mode in run_bgi(), and generates * the fractal surface and terrain model with the terrain_gen() function. * These models are simultaneously draWn to screen. Following this, the * terrain model is smoothed within the smooth_4stats() function (the name of * this function comes from its output of the elevation high, low, and * average values). The slope/aspect generation function values_gen() is * then performed to create arrays holding the slope and aspect values at * each point in the model. The initial watershed generation function * river_gen() is then performed; the model is now ready for the erosion * simulation. This simulation is carried out through the use of the erode() * function. Then the graphics mode is exited and the user is alerted to the * completion of the program. */ void main(int argc, char *argv[]) { char yesno; int initmode; erase(); textcolor(12); initmode = fg_getmode(); fg_setmode(initmode); fg_locate(0,26); fg_setcolor(12); fg_text("TERRAIN EVOLUTION SIMULATOR", 27); fg_locate(3,30); fg_text("M. Werneburg 1994-98", 20); fg_locate(6,0); /* These next few lines determine whether the user wishes to view the optional * introduction, export the generated terrain models at various stages * throughout the program) to Idrisi-format DEMs, or display rivers during * the draWing of the terrain model to screen. */ for (yesno = 1; yesno < 5; yesno++) { if (argv[argc-yesno][0] == ':') intro(); if (argv[argc-yesno][0] == '-') { if (argv[argc-yesno][1] == 'o') ask = 0; else if (argv[argc-yesno][1] == 'r') rivers = 0; } } textcolor(15); fg_locate(15,0); fg_setcolor(15); fg_text ("This model is capable of generating tilted terrains. This has the effect of", 76); fg_locate(16,0); fg_text ("adding a trend to the drainage pattern.", 39); fg_locate(17,0); fg_text ("Do you wish to add a tilt? [y/n] ", 33); //textcolor(8); // these textcolor routines draW the user's response in dark grey yesno = getche(); yesno = toupper(yesno); if (yesno == 'Y') { fg_locate(18,0); fg_text ("Please enter a tilt effect in millimetres. The terrain will be titled so", 73); fg_locate(19,0); fg_text ("that the top of the screen is raised while the bottom is lowered.", 65); fg_locate(20,0); fg_text ("Value (mm):", 11); scanf ("%d", &tiltz); // cprintf ("%d", tiltz); tilt=1; } textcolor(15); /*textcolor(15); cprintf("You may load an Idrisi-generated Digital Elevation Model into the simulator.nr"); cprintf("This file must be of dimensions 127 x 127, and it must have been recorded innr"); cprintf("integer format. Do you wish to load a 127 x 127 DEM for simulation?[y/n]:_b"); textcolor(8); // these textcolor routines draW the user's response in dark grey yesno = getche(); yesno = toupper(yesno); textcolor(15);*/ /* The following if statement determines the course of the program as defined * by the user's wish to read in a DEM created by the Idrisi GIS. If the * result is that a DEM is desired, the lone function get_terrain() are * called. Otherwise, the other set (run_bgi(), terrain_gen(), * smooth_4stats(), and values_gen()) are run. These two courses reflect * the two possible paths in establishing the terrain model and the * evolution process: in the first case, values are read from the DEM, and * the slope, aspect, and watershed values are calculated from that; in the * second case, the fractal modelling network (see text) is first built, then * the terrain model is built upon that, and lastly, the slope, aspect, and * watershed values are derived. In either case the initialize() function * is enacted first, and the values_gen() and erode() functions, which * determine the slope, aspect, and watershed arrays in the first case, and * perform the actual erosion simulation in thesecond. The final function * calls (closegraph() and the farfree() functions) do the necessary wrap-up, * first returning to text mode, then closing the dynamically allocated * memory blocks used in the simulator.*/ initialize(); //if (yesno == 'Y') // get_terrain(); //else // { setup_graphics(); terrain_gen(); smooth_4stats(); // } values_gen(1); erode(); fg_setmode(old_vmode); erase(); farfree(f_model); farfree(terrain); farfree(river_a); farfree(delta_z); farfree(entrain); farfree(lakes); textcolor(8); printf("rtt"); cprintf("The TERRAIN EVOLUTION SIMULATOR by Michael Werneburg"); textcolor(7); } /////////////////////////////////////INITIALIZE()//FUNCTION/////////////////////////////////////////////////////////////////////// /* This function initializes the arrays f_model[] and terrain[]. These arrays are used to hold the fractal modelling network and terrain elevation values, respectively. The first array is set to 15 for all points, and the second to 0. The reasoning for this is as follows. The function which follows initialize() - terrain_gen() - hunts the array f_model[] for zero values in its first pass. As there are none but that at the point (63,63) (assigned in the second last line of the function), only that position is selected in the first iteration of the function. The terrain[] array is set to 0 so that all points in the array have no elevation value at the outset of the terrain_gen() function. */ int initialize() { int array_place; // set memory block size to 127 by 127 f_model = (unsigned char huge*) farcalloc(16129, 1); if (f_model == NULL) { printf("error creating f_model array - insufficient memory!n"); return(0); } terrain = (int huge*) farcalloc(16129, 2); if (terrain == NULL) { printf("error creating terrain array - insufficient memory!n"); return(0); } lakes = (int huge*) farcalloc(16129, 2); if (lakes == NULL) { printf("error creating lakes array - insufficient memory!n"); return(0); } for(array_place = 0; array_place < 16129; array_place++) { f_model[array_place] = 15; // flat fractal model framework terrain[array_place] = 0; // no terrain lakes[array_place] = 0; // no lakes } // 63 * 127 = 8001 + 63 = 8064 f_model [8064] = 0; // assigns a zero to (63,63) terrain [8064] = 75; // no elevation at (63,63) return(1); // This line indicates to main() that the function worked. } /////////////////////////////////SETUP_GRAPHICS()//FUNCTION/////////////////////////////////////////////////////////////////////// // This function sets up the video mode, depending on system capabilities. void setup_graphics() { unsigned char colour_NO; fg_svgainit(0); if (fg_testmode(25, 0)) { // tests for 640 x 480 resolution with 256 colours. old_vmode = fg_getmode(); fg_setmode(25); } else if (fg_testmode(18,0)) { // tests for 640 x 480 resolution with 16 colours if 256 colour mode not possible old_vmode = fg_getmode(); fg_setmode(18); } else { printf ("nThe Terrain Evolution Simulator requires VGA graphics capabilities"); printf ("nto perform."); exit(1); } //if (fg_getmode() == 25) // { setrgb(255, 63, 63, 63); // colour 255 is white setrgb(254, 63, 0, 0); // colour 254 is red // setrgb(253, 0, 63, 63); // colour 253 is cyan setrgb(253, 0, 53, 53); // colour 253 is cyan setrgb(252, 63, 63, 0); // colour 252 is yellow setrgb(251, 63, 0, 63); // colour 251 is purple setrgb(151, 0, 0, 0); // black in 236 colour mode setrgb(152, 0, 0, 63); setrgb(153, 63, 0, 0); // colour 153 is red setrgb(154, 0, 63, 0); setrgb(155, 20, 20, 20); setcolor(255); // } //else // { // setrgb(15, 0, 63, 63); // colour 15 is cyan // setcolor(15); // } fg_box (0,200,0,40); fg_justify(0,0); move_to(100, 20); fg_print("Initializing video setup", 24); //if (fg_getmode() == 25) // { for (colour_NO = 1; colour_NO < 55; colour_NO++) // 50 shades of grey setrgb(colour_NO, (1.26*colour_NO), (1.26*colour_NO), (1.26*colour_NO)); for (colour_NO = 50; colour_NO < 100; colour_NO++) // 50 shades of yellow-green setrgb(colour_NO, (0.63*(colour_NO-50)), (1.26*(colour_NO-50)), 0); for (; colour_NO < 150; colour_NO++) // 50 shades of green setrgb(colour_NO, 0, (1.26*(colour_NO-100)), (0.63*(colour_NO-100))); // } //else // { // for (colour_NO = 1; colour_NO < 15; colour_NO++) // setrgb(colour_NO, colour_NO*4, colour_NO*4, colour_NO*4); // } maxx = fg_getmaxx(); // should be 640 maxy = fg_getmaxy(); // should be 480 lines = fg_getlines(); // should be 30 text lines on screen move_to(100, 20); setcolor(0); fg_box (0,200,0,40); fg_print("Initializing video setup", 24); setcolor(155); fg_box(127, 511, 47, 431); setcolor(153); } //////////////////////////////////TERRAIN_GEN()//FUNCTION///////////////////// /* The following block of nested for loops generates the fractal terrain. * This is accomplished through searching the f_model[] array for a value * equal to one less than the count of the iteration. That is, on the first * pass (where the iteration counter "iter8") is equal to one, the array is * searched for positions with the value 0. In this manner the position (63, * 63) is found during the fisrt iteration. In the following iterations - * indicated with an "iter8" value of 2 through 6 - values in the f_model[] * array of 1 through 5 are hunted. With the discovery of the positions * containing the proper iteration value(s), eight surrounding points in the * array are assigned an iteration value equal to iter8. This indicates * that these points are of the next lesser heirarchy (31, 31), (31, 63), (31, 94), (94, 63), * (94, 94), (63, 94), (31, 94), and (31, 63) are found from the current * central point of (63, 63). In the second iteration, the eight points * surrounding each of these first eight positions are assigned with the * D_xy value of 16, assigned in LOOP LINE 4. Since a number of these points * overlap, only 40 positions are determined in iteration two, instead of the * 64 expected (from 8 first-iteration points * 8 second iteration points). * This goes on until the fifth generation, at which point the scale value is * one, and all positions in the array are determined. As each of the new * positions are found, an elevation value is calculated for that position. * This elevation value is derived by adding or subtracting a number randomly * generated for the new position to the elevation value stored in the point * of higher heirarchy from which the new point's coordinates were found. * In this fashion, all elevation values across the entire terrain[] array are * ultimately derived from the value of 0 (which is equal to a real elevation * of 25 m, being 25000 mm above the zero elevation at -25000 in TES units) * seeded at (63, 63). For each random number generation, the range of * possible values is determined by the current iteration, as an argument * passed to the vari8() function listed below. For the first iteration, * the range of values is -22400 to 22400; for the second it is -11200 to * 11200, and so on, until in the final iteration the vertical range is -700 * to 700. This halving of vertical range with every iteration is in keeping * with the halving of the horizontal distance involved in the derivation of * the new points' array positions. Lastly, it was realized that the three * positions with a lesser Y value, as well the two positions with a lesser * X value than the central position in each iteration following the first * may have already been determined by a previous interpolation of the same * iteration. * This is possible due to the overlap of points based on the scaling system. * To demonstrate, the three positions (31,47), (47,47), and (47,63) found in * the second iteration share as a common third-iteration "offspring" * position the point (39, 55). An elevation value for this point will be * first assigned as the value of (31, 47) +/- a random value between -11200 * and 11200. It will then be overwritten when the program re-evaluates it * as the value of (47, 47) +/- Û 11200, and again when the points around * (47, 63) are determined. This problem gets worse in frequency with every * iteration. To counteract this problem without creating as many as four * more memory-swallowing arrays of 127 x 127 two-byte arrays, part of the * code below checks first for a pre-existing value for each new position; * it is for this reason that the terrain[] array was first set to zero for * all points. * If a non-zero value is found, it averages the elevation value it was going * to assign (as found by the process described above) with the pre-existing * value. This means that the value found from the (31, 47) position would * be averaged with that found from (47, 47). This averaged value would then * be averaged with that found from (63, 47), putting an unfortunate but * inescapable weight on the last value. * The procedure of averaging outlined above was the most realistic one found * in creating a model by which each point could be traced in genesis to a * nearby base point. */ void terrain_gen() { int rnd_value, // holder of randomly generated numbers D_xy=64, // the scale of X and Y change with every iteration, defined as 32 divided by 2 to the power of iter8 tg_row, tg_XY, dx, dy, iter8; // counter of iterations in generation of f_model[] array unsigned int aplace; randomize(); // this function initializes the random # generation routine for use. //text_xy(253, maxx/2, maxy-30, "Generating fractal modelling framework"); for (iter8 = 1; iter8 < 7; iter8++) // LOOP LINE 1 { D_xy/=2; for (X = 0; X < 127; X++) // LOOP LINE 9 { tg_row=X*127; for (Y = 0; Y < 127; Y++) // LOOP LINE 11 { tg_XY=tg_row+Y; if (f_model[tg_XY] == (iter8 - 1)) // LOOP LINE 13 { /* These lines assign four of the eight f_model[] heirarchy values found * in each iteration, for those positions diagonally * surrounding the central points found with each iteration. */ /* The following lines assign the f_model[] values to points in the array * which are of a distance set by D_xy and on the same horizontal or * vertical. */ /* These lines create the elevation values at each of the eight positions * for which f_model[] values were found above. */ for (dx=-1;dx<2;dx++) { for (dy=-1;dy<2;dy++) { if ((dx || dy) && SANITY) { aplace=((X+(D_xy*dx))*127)+(Y+(D_xy*dy)); f_model[aplace] = iter8; if (terrain[aplace]) terrain [aplace] = (vari8(iter8,X,Y) + terrain[aplace])/2; else terrain [aplace] = vari8(iter8,X,Y); // if (f_model[tg_XY] == (iter8-1)) // draw(X, Y, (49 - (f_model[tg_XY]*8))); } } } if (tiltz) terrain[aplace] += (int) ((Y - 63) * (tiltz/63.0)); } } // end of y for loop } // end of x for loop } // end of for iter8 loop. /* These lines draw the legend of the f_model[] array to the screen beside the * array. */ /*set_colour(155); draw_rect(2, 118, 49, 178); set_colour(254); fg_box(0, 120, 47, 180); set_colour(255); fg_box(1, 119, 48, 179); draw (-8, 7, 49);draw (-7, 7, 49);draw (-8, 8, 49);draw (-7, 8, 49); text_xy(151, 50, 73,"Iteration 1"); draw (-8, 12, 41);draw (-7, 12, 41);draw (-8, 13, 41);draw (-7, 13, 41); text_xy(151, 50, 88,"Iteration 2"); draw (-8, 17, 33);draw (-7, 17, 33);draw (-8, 18, 33);draw (-7, 18, 33); text_xy(151, 50, 103,"Iteration 3"); draw (-8, 22, 25);draw (-7, 22, 25);draw (-8, 23, 25);draw (-7, 23, 25); text_xy(151, 50, 118,"Iteration 4"); draw (-8, 27, 17);draw (-7, 27, 17);draw (-8, 28, 17);draw (-7, 28, 17); text_xy(151, 50, 133,"Iteration 5"); draw (-8, 32, 9); draw (-7, 32, 9); draw (-8, 33, 9); draw (-7, 33, 9); text_xy(151, 50, 148,"Iteration 6"); text_xy (151, 55, 170, "LEGEND"); text_xy (0, maxx/2, maxy-30, "Generating fractal modelling framework"); text_xy (253, maxx/2, maxy-30, "Fractal modelling framework complete"); text_xy (253, maxx/2, maxy-15, "press any key to continue"); noise = getche(); text_xy (0, maxx/2, maxy-15, "press any key to continue"); set_colour(0); draw_rect(0, 120, 47, 180); */ text_xy (253, maxx/2, maxy-15, "Generating terrain model"); /* The following lines draw the terrain[] model to the screen, overwriting the * previously drawn f_model[] array. */ for (X = 0; X < 127; X++) { tg_row=X*127; for (Y = 0; Y < 127; Y++) { tg_XY=tg_row+Y; if (terrain[tg_XY] > 0) { // draw (X, Y, (int)(24.5 + (terrain[tg_XY]/1000.0))); draw (X, Y, (int)(50 + 3.5 * (terrain[tg_XY]/1000.0))); if (lakes[tg_XY]) ldraw (X, Y, 152); } else draw (X, Y, 253); // Ocean... } } //legend(0); // This call causes the legend to be drawn. text_xy (0, maxx/2, maxy-30, "Fractal modelling framework complete"); text_xy (0, maxx/2, maxy-15, "Generating terrain model"); text_xy (253, maxx/2, maxy-30, "Terrain model generated"); text_xy (253, maxx/2, maxy-15, "Press any key to continue"); noise = getche(); text_xy (0, maxx/2, maxy-15, "Press any key to continue"); text_xy (253, maxx/2, maxy-15, "Smoothening terrain model"); } ///////////////////////////////////VARI8()//FUNCTION///////////////////////// /* This is the function - called in the terrain_gen() function above - which * generates a new elevation value, which it then returns to the calling * function. It takes as arguments the current scale factor (denoted in the * parent function terrain_gen() as iter8), as well as the array position * from which the currently active array position is derived; that is, * (countX, countY) in terrain_gen(). The function begins by assigning the * elevation value of the parent array position (midX, midY) to old_value. * It then creates a new value randomly, defining the range through the * "scale" variable, which is passed to the function as iter8() by * terrain_gen(). If the random number generated was greater than half of * the range, it (rnd_value) is added to old_value. If not, it is subtracted. * The resulting value is then assigned to new_value, which is then * returned to the calling function. */ // modified as of Nov 3, 1995 - I'm trying a new formula for the vertical // difference with each pass. The basic formula is // z(max)=44800/(1+(pow(0.5,(scale-1))) int vari8(unsigned char scale,unsigned char midX,unsigned char midY) { int new_value, old_value, var_factor; unsigned int rnd_value; var_factor= (int) 10000 / (1+(pow(0.75, (scale -1)))); old_value = terrain[(midX*127) +midY]; // the elevation of the central calling position is copied to "old_value". rnd_value = random(var_factor); // the random number (between zero and the value "var_factor") is assigned to rnd_value new_value = old_value + ((var_factor/2) - rnd_value); // The new elevation value is assigned by adding a modified rnd_value to // old_value. Since the modified rnd_value may be below zero, the new value // may be below old_value. return(new_value+500); // This line returns the new elevation value to the calling function. } //////////////////////////////////SMOOTH_4STATS()//FUNCTION////////////////// /* The following code is designed to smoothen the terrain, to make its slopes * more continuous, without jeopardizing the overall shape of the surface. * This is accomplished through the use of a travelling window as outlined in * the blurb above the values_gen() code. The window is used such that the * current active position is set to the average of the four, six, or nine * positions in the active window. This smoothing removes the unrealistic * sharpness and obviously artificial terrain shapes which are a result of the * terrain model generation process. */ int smooth_4stats() { long sum; // the sum of all terrain elevation values for positions in the window. unsigned int divisor; // used to derive an average of all window elevation values as a divisor of "sum". This is assigned a // value by determining the position of the active window in the river_a array. char across, // the horizontal count across the active window. along, // the vertical count across the active window pass; // number of passes the smoothing algorithm makes int s4_row, s4_XY; int huge *tsmooth; // memory block used to hold terrain values during smoothening tsmooth = (int huge*) farcalloc(16129, 2); if (terrain == NULL) { printf("error creating tsmooth array - insufficient memory!n"); return(0); } //setfillstyle(1,0); /* This code determines the position of the active window within the river_a * array by finding the position of the active central point. In this manner, * a window in the corner is assigned a 2 x 2 dimension for the calucation * of sum and divisor values, an edge window is assigned a shape of 3 x 2 * positions, and all other windows 3 x 3. */ for (pass = 0; pass < 2; pass++) { for (X = 0; X < 127; X++) { s4_row=X*127; for (Y = 0; Y < 127; Y++) { s4_XY=s4_row+Y; sum = 0; if (X == 0) { A = 0; B = 1; divisor = 2; } else if (X == 126) { A = -1; B = 0; divisor = 2; } else if ((X > 0) && (X < 126)) { A = -1; B = 1; divisor = 3; } else { A = 0; B = 0; divisor=1; } if (Y == 0) { C = 0; D = 1; divisor *= 2; } else if (Y == 126) { C = -1; D = 0; divisor *= 2; } else if ((Y > 0) && (Y < 126)) { C = -1; D = 1; divisor *= 3; } else { A = 0; B = 0; divisor=1; } /* Using the parameters determined above, the sum of all elevation values * within the active window is found, and then divided by the appropriate * value to come up with a sum. This is then assigned to the entrain model, * and added to the running total avecount for later division by the size of * the array (16129 positions). */ for (across = A; across <= B; across++) for (along = C; along <= D; along++) sum += terrain[((X+across)*127) + (Y+along)]; tsmooth[s4_row + Y] = (int)(sum / divisor); // tsmooth[s4_row+Y] = terrain[s4_row+Y]; } } NOSMOOTH: /* These lines draw the smoothed terrain model[] to the screen, and record * the lowest and highest elevation values in the array. */ for (X = 0; X < 127; X++) { s4_row=X*127; for (Y = 0; Y < 127; Y++) { s4_XY=s4_row+Y; terrain[s4_XY] = tsmooth[s4_XY]; if (pass == 1) { // 24.5 is just a colour.... if (terrain[s4_XY] > 0) { // draw (X, Y, (int)(24.5 + (terrain[s4_XY]/1000.0))); draw (X, Y, (int)(50 + 3.5 * (terrain[s4_XY]/1000.0))); if (lakes[s4_XY]) ldraw (X,Y, 152); } else draw (X,Y, 253); // Ocean if (tsmooth[s4_XY] < lowcount) lowcount = tsmooth[s4_XY]; else if (tsmooth[s4_XY] > highcount) highcount = tsmooth[s4_XY]; avecount += tsmooth[s4_XY]; } } } } // end of pass loop farfree(tsmooth); text_xy (0, maxx/2, maxy-30, "Terrain model generated"); text_xy (0, maxx/2, maxy-15, "Smoothening terrain model"); text_xy (253, maxx/2, maxy-30, "Terrain model smoothened"); text_xy (253, maxx/2, maxy-15, "press any key to continue"); // These lines display the highest, lowest, and average elevation values in the model. hold_v = highcount/1000.0; sprintf (output_h, "Highest elevation %.1f m", hold_v); text_xy (253, maxx/2, 10, output_h); hold_v = avecount / (16129000.0); sprintf (output_a, "Average elevation %.1f m", hold_v); text_xy (253, maxx/2, 22, output_a); hold_v = lowcount/1000.0; sprintf (output_l, "Lowest elevation %.1f m", hold_v); text_xy (253, maxx/2, 35, output_l); noise = getch(); text_xy (0, maxx/2, maxy-15, "press any key to continue"); text_xy (0, maxx/2, 10, output_h); text_xy (0, maxx/2, 22, output_a); text_xy (0, maxx/2, 35, output_l); // These lines offer the choice to record the DEM as it currently is to a DEM. if (!(ask)) { text_xy (253, maxx/2, maxy-15, "Do you wish to save this model to an Idrisi DEM? [y/n]"); noise = getch(); text_xy (0, maxx/2, maxy-15, "Do you wish to save this model to an Idrisi DEM? [y/n]"); noise = toupper(noise); if (noise == 'Y') terrain_rec(); } text_xy (253, maxx/2, maxy-15, "Now generating aspect and slope values"); return(1); } /////////////////////////////////VALUES_GEN()//FUNCTION/////////////////////// /* This function assigns slope and aspect values to each position in the * terrain[] array. These are held in the delta_z[] and f_model[] arrays, * respectively. These values are then used to determine the watershed which * empties into each point in the array, which is equated linearly to the * flow of water through that 10m x 10m position. At the same time, the * amount of material which would be transported from each position - a * function of the already derived slope and watershed values for that * position. Note that f_model[] is reused here out of an unwillingness to * pile up the variables. */ int values_gen(unsigned char isit_first) { int new_value; // variable used to identify locally lowest elevation & highest local slope. int X1, // temporary X position as watersheds are calculated Y1, // temporary Y position as watersheds are calculated DX, // horizontal distance across window DY; // vertical distance across window double base, // holder of slope value for power function "pow()" exponent = 1.4; // holder of exponent of slope for erosional equation int min_dz = 1000, vg_row, tdx, tdy, A1, B1, C1, D1, lakec, toggle, max_dz = 0; // minimum and maximum slope value holders. unsigned int vg_XY, aplace; /* During the first pass of this function only, these lines initiate the * river_a[], delta_z[], and entrain[] arrays. */ if (isit_first) { river_a = (unsigned char huge*) farcalloc(16129, 2); if (river_a == NULL) { printf("error creating river_a array - insufficient memory!n"); noise = getch(); return(0); } delta_z = (int huge*) farcalloc(16129, 2); if (delta_z == NULL) { printf("error creating delta_z array - insufficient memory!n"); noise = getch(); return(0); } entrain = (int huge*) farcalloc(16129, 2); if (entrain == NULL) { printf("error creating entrain array - insufficient memory!n"); noise = getch(); return(0); } } for (X = 0; X < 127; X++) { vg_row=X*127; for (Y = 0; Y < 127; Y++) { vg_XY=vg_row+Y; river_a[vg_XY] = 1; // initializes the points in river_a[] to one. delta_z[vg_XY] = 0; // initializes the points in delta_z[] to zero. entrain[vg_XY] = 0; // initializes the points in entrain[] to zero. /* The following lines determine the slope at each point in the array toward * the aspect of greatest down-hill direction. These are assigned to the * delta_z[] array for all points. The aspect of these recorded slopes are * kept in the array f_model[] */ if (X == 0) { A = 0; B = 1; } else if (X == 126) { A = -1; B = 0; } else { A = -1; B = 1; } if (Y == 0) { C = 0; D = 1; } else if (Y == 126) { C = -1; D = 0; } else { C = -1; D = 1; } f_model[vg_XY]=0; for (DX = A; DX <= B; DX++) { for (DY = C; DY <= D; DY++) { if (DX || DY) { aplace=((X+DX)*127)+Y+DY; switch (DX*DY) { case (-1): // ...or case (1): // We must fake a little trigonometry new_value = (int)((((terrain[vg_XY]) - terrain[aplace])/1.414213562) + 0.5); break; case (0): // Simple math.... new_value = (terrain[vg_XY]) - terrain[aplace]; break; } if (new_value > delta_z[vg_XY]) { delta_z[vg_XY] = new_value; if ((DX == -1) && (DY == -1)) f_model[vg_XY] = 1; // the aspect value for each position is else if ((DX == 0) && (DY == -1)) f_model[vg_XY] = 2; // found with these lines, which determine else if ((DX == 1) && (DY == -1)) f_model[vg_XY] = 3; // the direction based on the X and Y offset, else if ((DX == 1) && (DY == 0)) f_model[vg_XY] = 4; // which may range from -1 to 1 each. The else if ((DX == 1) && (DY == 1)) f_model[vg_XY] = 5; // directions assigned range from 1 through 8 else if ((DX == 0) && (DY == 1)) f_model[vg_XY] = 6; // starting in the upper left. else if ((DX == -1) && (DY == 1)) f_model[vg_XY] = 7; else if ((DX == -1) && (DY == 0)) f_model[vg_XY] = 8; } } } } // end of for DY // if (delta_z[vg_row+Y] == 0) // we're on a plain, so.. // f_model[vg_row+Y] = 0; // there is no aspect if (f_model[vg_XY] == 0) { for (DX = A; DX <= B; DX++) { for (DY = C; DY <= D; DY++) { if ((DX || DY) && ((Y+DY)<127)&&((X+DX)<127)&&((X+DX)>=0)&&((Y+DY)>=0)) { aplace=((X+DX)*127)+Y+DY; if ((terrain[vg_XY] + river_a[vg_XY]) > terrain[aplace] && (terrain[vg_XY] <= terrain[aplace])) { lakes[aplace]=terrain[vg_XY] + river_a[vg_XY]; } else lakes[aplace]=0; } } } } /* The following lines record the minimum and maximum slopes found, a process * which is only enacted in the first call of the function, since doing so * each time would slow the evolution process. */ if (isit_first) { hold_v += delta_z[vg_XY]; if (delta_z[vg_row+Y] > max_dz) max_dz = delta_z[vg_XY]; else if (delta_z[vg_XY] < min_dz) min_dz = delta_z[vg_XY]; } } } /* These torturous lines add the ponding effect. This has a look at the * entire terrain, and finds local lows. It then checks a ring of positions * around the local low, and determines whether those spots are above the * low in elevation, but below the elevation of the water filling the low. * This is how pools form. */ for (X=0;X<127;X++) { for (Y=0;Y<127;Y++) { vg_XY=(X*127)+Y; if (lakes[vg_XY]) { A=0; B=0; C=0; D=0; lakec=1; while (lakec) { toggle=0; // Has at least one value changed? if ((X+A)>0) { A--; toggle++; } if ((X+B)<126) { B++; toggle++; } if ((Y+C)>0) { C--; toggle++; } if ((Y+D)<126) { D++; toggle++; } lakec=0; if (!toggle) break; for (DX=A; DX<=B; DX++) { for (DY=C; DY<=D;DY++) { if ((DX==A || DX==B) || (DY==C || DY==D)) { aplace=(X+DX)*127+(Y+DY); if (terrain[aplace] < lakes[vg_XY] && terrain[aplace] >= terrain[vg_XY] && lakes[aplace] != lakes[vg_XY]) { for (tdx=-1; tdx<=1;tdx++) { for (tdy=-1; tdy<=1; tdy++) { if (lakes[(X+DX+tdx)*127+(Y+DY+tdy)] == lakes[vg_XY]) { lakec++; lakes[aplace] = lakes[vg_XY]; } } } } } } } } } } } // These lines write the minimum, maximum, and mean slope values in the delta_z[] array to the screen. if (isit_first) { base = tan((double)(max_dz/10000.0)) * 57.29577952; sprintf(output_h, "Highest slope %.1f¯", (float)base); text_xy(253, maxx/2, 10, output_h); hold_v /= 16129; base = tan((double)(hold_v/10000.0)) * 57.29577952; sprintf(output_a, "Average slope %.1f¯", (float)base); text_xy(253, maxx/2, 22, output_a); base = tan((double)(min_dz/10000.0)) * 57.29577952; sprintf(output_l, "Lowest slope %.1f¯", (float)base); text_xy(253, maxx/2, 35, output_l); text_xy (0, maxx/2, maxy-30, "Terrain model smoothened"); text_xy (0, maxx/2, maxy-15, "Now generating aspect and slope values"); text_xy (253, maxx/2, maxy-30, "Slope and aspect model generated"); text_xy (253, maxx/2, maxy-15, "press any key to continue"); noise = getch(); text_xy (0, maxx/2, maxy-30, "Slope and aspect model generated"); text_xy (0, maxx/2, maxy-15, "press any key to continue"); text_xy (253, maxx/2, maxy-30, "Now generating watersheds & sediment values"); // text_xy (253, maxx/2, maxy-15, "This takes time - the dot traces progress"); } /* This code determines the watershed values at each point by searching * through the entire array for non-zero slopes. For each such value found, * a subroutine follows the path downhill until a zero slope value is found. * With every position found, the watershed value is incremented by one, to * indicate that the first non-zero slope value position (at which the * subroutine began) is contained within that position's watershed. This is * accomplished through the determination of DX and DY from the aspect * direction given in f_model[] for each position. These values, in turn, * are added to the active X and Y coordinates so that the new, downhill * (X,Y) are found, and then checked for slope value. In the case of * non-zero slope, the value in f_model[] for that position is incremented by * one, and the next downslope position is sought. For zero slope positions * (local lows), no downslope position is sought.*/ for (X = 0; X < 127; X++) { for (Y = 0; Y < 127; Y++) { X1 = X; Y1 = Y; if (terrain[X1*127+Y1] > 0) { do { vg_row=X1*127; switch ((f_model[vg_row + Y1] % 10)) { case 1: DX = -1; DY = -1; break; // These lines recreate the X and Y offset values used in the case 2: DX = 0; DY = -1; break; // previous function for use below. case 3: DX = 1; DY = -1; break; case 4: DX = 1; DY = 0; break; case 5: DX = 1; DY = 1; break; case 6: DX = 0; DY = 1; break; case 7: DX = -1; DY = 1; break; case 8: DX = -1; DY = 0; break; case 0: DX = 0; DY = 0; break; } // The following line adds one to the watershed value stored in river_a[] to the position downhill of the active position. if (terrain[((X1+DX)*127)+(Y1+DY)] > 0) river_a[((X1+DX)*127)+(Y1+DY)]++; /* The following lines calculate the erosion due to wash and rilling which take * place at each position. If several positions drain to the same position, * that position is passed through by the algorithm several times, and the * wash/rilling effect would be calculated each time, if the position were * not flagged somehow after the first calculation. The flag used is the * addition of ten to the f_model[] array for that position. While this has * no effect on the aspect values held in that array, it allows for the * determination of a previous calculation of that position's wash/rill * erosion effects. The equation used is Î = (tan1.4‡ * L0.6)/10, where * ‡ is slope angle, and L is ground length of slope. */ if (f_model[(X1*127)+Y1] < 10) { base = (double) (delta_z[vg_row+Y1]/10000.0); new_value = (int)(0.5 + (25.11886431 * pow(base, exponent))); entrain[vg_row+Y1] += new_value; entrain[vg_row+Y1+DY] -= new_value; f_model[vg_row+Y1] += 10; } /* if (isit_first) { if (terrain[vg_row+Y] > 0) { //draw(X1,Y1,(int)(24.5 + (terrain[vg_row+Y1]/1000.0))); draw (X, Y, (int)(50 + 3.5 * (terrain[vg_XY]/1000.0))); if (lakes[vg_row+Y1]) ldraw(X1,Y1,152); } else draw(X1,Y1,253); // Ocean }*/ X1 += DX; Y1 += DY; } while ((f_model[vg_row+Y1]) && (X1>=0) && (Y1>=0) && (X1<=126) && (Y1<=126) && (delta_z[vg_row+Y1])); } } } /* At the end of every running of the values_gen() function, the following * code draws the new terrain[] model to the screen, including the streams * which might form upon the surface. No streams are assumed to form below * a watershed area of 1000m˝, which yields a flow of 400 cubic metres per * year if annual rainfall is assumed to approximate that of southern Ontario * (800mm/yr), and that half of all precipitation is lost to evaporation