/**************************************************************************** * * MODULE: r.path * * AUTHOR(S): based on r.drain * Markus Metz * * PURPOSE: Tracing paths from starting points following * input directions. * COPYRIGHT: (C) 2017 by the GRASS Development Team * * This program is free software under the GNU General Public * License (>=v2). Read the file COPYING that comes with GRASS * for details. * *****************************************************************************/ #include #include #include #include #include #include /* for using the "open" statement */ #include #include #include /* for using the close statement */ #include #include #include #include #include #include "local.h" #include "pavlrc.h" /* should probably be updated to a pointer array & malloc/realloc as needed */ #define POINTS_INCREMENT 1024 /* start points */ struct point { int row; int col; double value; struct point *next; }; /* stack points for bitmask directions */ struct spoint { int row; int col; int dir; double value; struct spoint *next; }; /* output path points */ struct ppoint { int row; int col; double value; }; /* managed list of output path points */ struct point_list { struct ppoint *p; int n; int nalloc; }; int dir_degree(int dir_fd, int val_fd, struct point *startp, struct Cell_head *window, struct Map_info *Out, struct point_list *pl, int out_mode); int dir_bitmask(int dir_fd, int val_fd, struct point *startp, struct Cell_head *window, struct Map_info *Out, struct point_list *pl, int out_mode); void pl_add(struct point_list *, struct ppoint *); /* comparing path points with qsort */ int cmp_pp(const void *a, const void *b) { const struct ppoint *ap = (const struct ppoint *)a; const struct ppoint *bp = (const struct ppoint *)b; /* 1. ascending by row */ if (ap->row != bp->row) return (ap->row - bp->row); /* 2. ascending by col */ if (ap->col != bp->col) return (ap->col - bp->col); /* 3. descending by value */ if (ap->value > bp->value) return -1; return (ap->value < bp->value); } int main(int argc, char **argv) { int fd, dir_fd, val_fd; int i, j, have_points = 0; int nrows, ncols; int npoints; int out_id, dir_id, dir_format; struct FPRange drange; DCELL dmin, dmax; char map_name[GNAME_MAX], out_name[GNAME_MAX], dir_name[GNAME_MAX]; char *tempfile1, *tempfile2; struct History history; struct Cell_head window; struct { struct Option *dir; struct Option *format; struct Option *val; struct Option *coord; struct Option *vpoint; struct Option *rast; struct Option *vect; } opt; struct { struct Flag *copy; struct Flag *accum; struct Flag *count; } flag; struct GModule *module; void *dir_buf; struct point *head_start_pt = NULL; struct point *next_start_pt; struct point_list pl, *ppl; int start_row, start_col, out_mode; double east, north; struct line_pnts *Points; struct line_cats *Cats; struct Map_info vout, *pvout; char *desc = NULL; G_gisinit(argv[0]); module = G_define_module(); G_add_keyword(_("raster")); G_add_keyword(_("hydrology")); G_add_keyword(_("cost surface")); module->description = _("Traces paths from starting points following input directions."); opt.dir = G_define_standard_option(G_OPT_R_INPUT); opt.dir->label = _("Name of input direction"); opt.dir->description = _("Direction in degrees CCW from east, or bitmask encoded"); opt.format = G_define_option(); opt.format->type = TYPE_STRING; opt.format->key = "format"; opt.format->label = _("Format of the input direction map"); opt.format->required = YES; opt.format->options = "auto,degree,45degree,bitmask"; opt.format->answer = "auto"; G_asprintf(&desc, "auto;%s;degree;%s;45degree;%s;bitmask;%s", _("auto-detect direction format"), _("degrees CCW from East"), _("degrees CCW from East divided by 45 (e.g. r.watershed directions)"), _("bitmask encoded directions (e.g. r.cost -b)")); opt.format->descriptions = desc; opt.val = G_define_standard_option(G_OPT_R_INPUT); opt.val->key = "values"; opt.val->label = _("Name of input raster values to be used for output"); opt.val->required = NO; opt.rast = G_define_standard_option(G_OPT_R_OUTPUT); opt.rast->key = "raster_path"; opt.rast->required = NO; opt.rast->label = _("Name for output raster path map"); opt.vect = G_define_standard_option(G_OPT_V_OUTPUT); opt.vect->key = "vector_path"; opt.vect->required = NO; opt.vect->label = _("Name for output vector path map"); opt.coord = G_define_standard_option(G_OPT_M_COORDS); opt.coord->key = "start_coordinates"; opt.coord->multiple = YES; opt.coord->description = _("Coordinates of starting point(s) (E,N)"); opt.coord->guisection = _("Start"); opt.vpoint = G_define_standard_option(G_OPT_V_INPUTS); opt.vpoint->key = "start_points"; opt.vpoint->required = NO; opt.vpoint->label = _("Name of starting vector points map(s)"); opt.vpoint->guisection = _("Start"); flag.copy = G_define_flag(); flag.copy->key = 'c'; flag.copy->description = _("Copy input cell values on output"); flag.copy->guisection = _("Path settings"); flag.accum = G_define_flag(); flag.accum->key = 'a'; flag.accum->description = _("Accumulate input values along the path"); flag.accum->guisection = _("Path settings"); flag.count = G_define_flag(); flag.count->key = 'n'; flag.count->description = _("Count cell numbers along the path"); flag.count->guisection = _("Path settings"); G_option_required(opt.rast, opt.vect, NULL); G_option_exclusive(flag.copy, flag.accum, flag.count, NULL); G_option_requires_all(flag.copy, opt.rast, opt.val, NULL); G_option_requires_all(flag.accum, opt.rast, opt.val, NULL); G_option_requires_all(flag.count, opt.rast, NULL); if (G_parser(argc, argv)) exit(EXIT_FAILURE); strcpy(dir_name, opt.dir->answer); *map_name = '\0'; *out_name = '\0'; if (opt.rast->answer) { strcpy(out_name, opt.rast->answer); if (opt.val->answer) strcpy(map_name, opt.val->answer); } pvout = NULL; if (opt.vect->answer) { if (0 > Vect_open_new(&vout, opt.vect->answer, 0)) { G_fatal_error(_("Unable to create vector map <%s>"), opt.vect->answer); } Vect_hist_command(&vout); pvout = &vout; } if (flag.copy->answer) out_mode = OUT_CPY; else if (flag.accum->answer) out_mode = OUT_ACC; else if (flag.count->answer) out_mode = OUT_CNT; else out_mode = OUT_PID; /* get the window information */ G_get_window(&window); nrows = Rast_window_rows(); ncols = Rast_window_cols(); npoints = 0; /* TODO: use r.cost method to create a list of start points */ if (opt.coord->answer) { for (i = 0; opt.coord->answers[i] != NULL; i += 2) { G_scan_easting(opt.coord->answers[i], &east, G_projection()); G_scan_northing(opt.coord->answers[i + 1], &north, G_projection()); start_col = (int)Rast_easting_to_col(east, &window); start_row = (int)Rast_northing_to_row(north, &window); if (start_row < 0 || start_row > nrows || start_col < 0 || start_col > ncols) { G_warning(_("Starting point %d is outside the current region"), i + 1); continue; } npoints++; have_points = 1; next_start_pt = (struct point *)(G_malloc(sizeof(struct point))); next_start_pt->row = start_row; next_start_pt->col = start_col; next_start_pt->value = npoints; next_start_pt->next = head_start_pt; head_start_pt = next_start_pt; } } if (opt.vpoint->answers) { for (i = 0; opt.vpoint->answers[i] != NULL; i++) { struct Map_info In; struct bound_box box; int cat, type; Points = Vect_new_line_struct(); Cats = Vect_new_cats_struct(); Vect_set_open_level(1); /* topology not required */ if (1 > Vect_open_old(&In, opt.vpoint->answers[i], "")) G_fatal_error(_("Unable to open vector map <%s>"), opt.vpoint->answers[i]); G_verbose_message(_("Reading vector map <%s> with start points..."), Vect_get_full_name(&In)); Vect_rewind(&In); Vect_region_box(&window, &box); box.T = box.B = 0; while (1) { /* register line */ type = Vect_read_next_line(&In, Points, Cats); /* Note: check for dead lines is not needed, because they are skipped by V1_read_next_line_nat() */ if (type == -1) { G_warning(_("Unable to read vector map")); continue; } else if (type == -2) { break; } if (!Vect_point_in_box(Points->x[0], Points->y[0], 0, &box)) continue; start_col = (int)Rast_easting_to_col(Points->x[0], &window); start_row = (int)Rast_northing_to_row(Points->y[0], &window); /* effectively just a duplicate check to G_site_in_region() ??? */ if (start_row < 0 || start_row > nrows || start_col < 0 || start_col > ncols) continue; npoints++; have_points = 1; next_start_pt = (struct point *)(G_malloc(sizeof(struct point))); next_start_pt->row = start_row; next_start_pt->col = start_col; Vect_cat_get(Cats, 1, &cat); next_start_pt->value = cat; next_start_pt->next = head_start_pt; head_start_pt = next_start_pt; } Vect_close(&In); /* only catches maps out of range until something is found, not after */ if (!have_points) { G_warning(_("Starting vector map <%s> contains no points in the current region"), opt.vpoint->answers[i]); } Vect_destroy_line_struct(Points); Vect_destroy_cats_struct(Cats); } } if (have_points == 0) G_fatal_error(_("No start point(s) specified")); /* determine the drainage paths */ /* get some temp files */ val_fd = -1; tempfile1 = NULL; if (opt.val->answer && opt.rast->answer) { DCELL *map_buf; G_verbose_message(_("Reading raster values map <%s> ..."), map_name); tempfile1 = G_tempfile(); val_fd = open(tempfile1, O_RDWR | O_CREAT, 0666); map_buf = Rast_allocate_d_buf(); fd = Rast_open_old(map_name, ""); /* transfer the values map to a temp file */ for (i = 0; i < nrows; i++) { Rast_get_d_row(fd, map_buf, i); if (write(val_fd, map_buf, ncols * sizeof(DCELL)) != ncols * sizeof(DCELL)) { G_fatal_error(_("Unable to write to tempfile")); } } Rast_close(fd); G_free(map_buf); } if (Rast_read_fp_range(dir_name, "", &drange) < 0) G_fatal_error(_("Unable to read range file")); Rast_get_fp_range_min_max(&drange, &dmin, &dmax); if (dmax <= 0) G_fatal_error(_("Invalid directions map <%s>"), dir_name); dir_format = -1; if (strcmp(opt.format->answer, "degree") == 0) { if (dmax > 360) G_fatal_error(_("Directional degrees can not be > 360")); dir_format = DIR_DEG; } else if (strcmp(opt.format->answer, "45degree") == 0) { if (dmax > 8) G_fatal_error(_("Directional degrees divided by 45 can not be > 8")); dir_format = DIR_DEG45; } else if (strcmp(opt.format->answer, "bitmask") == 0) { if (dmax > (1 << 16) - 1) G_fatal_error(_("Bitmask encoded directions can not be > %d"), (1 << 16) - 1); dir_format = DIR_BIT; } else if (strcmp(opt.format->answer, "auto") == 0) { if (dmax <= 8) { dir_format = DIR_DEG45; G_important_message(_("Input direction format assumed to be degrees CCW from East divided by 45")); } else if (dmax <= (1 << 8) - 1) { dir_format = DIR_BIT; G_important_message(_("Input direction format assumed to be bitmask encoded without Knight's move")); } else if (dmax <= 360) { dir_format = DIR_DEG; G_important_message(_("Input direction format assumed to be degrees CCW from East")); } else if (dmax <= (1 << 16) - 1) { dir_format = DIR_BIT; G_important_message(_("Input direction format assumed to be bitmask encoded with Knight's move")); } else G_fatal_error(_("Unable to detect format of input direction map <%s>"), dir_name); } if (dir_format <= 0) G_fatal_error(_("Invalid directions format '%s'"), opt.format->answer); G_verbose_message(_("Reading direction map <%s> ..."), dir_name); dir_id = Rast_open_old(dir_name, ""); tempfile2 = G_tempfile(); dir_fd = open(tempfile2, O_RDWR | O_CREAT, 0666); if (dir_format == DIR_BIT) { dir_buf = Rast_allocate_c_buf(); for (i = 0; i < nrows; i++) { Rast_get_c_row(dir_id, dir_buf, i); if (write(dir_fd, dir_buf, ncols * sizeof(CELL)) != ncols * sizeof(CELL)) { G_fatal_error(_("Unable to write to tempfile")); } } } else { dir_buf = Rast_allocate_d_buf(); for (i = 0; i < nrows; i++) { Rast_get_d_row(dir_id, dir_buf, i); if (dir_format == DIR_DEG45) { DCELL *dp; dp = (DCELL *)dir_buf; for (j = 0; j < ncols; j++, dp++) *dp *= 45; } if (write(dir_fd, dir_buf, ncols * sizeof(DCELL)) != ncols * sizeof(DCELL)) { G_fatal_error(_("Unable to write to tempfile")); } } } Rast_close(dir_id); G_free(dir_buf); ppl = NULL; if (opt.rast->answer) { /* raster output */ pl.p = NULL; pl.n = 0; pl.nalloc = 0; ppl = &pl; } /* determine the drainage paths */ /* repeat for each starting point */ G_verbose_message(_("Processing start points...")); next_start_pt = head_start_pt; while (next_start_pt) { /* follow directions from start points to determine paths */ /* path tracing algorithm selection */ if (dir_format == DIR_BIT) { struct Map_info Tmp; if (pvout) { if (Vect_open_tmp_new(&Tmp, NULL, 0) < 0) G_fatal_error(_("Unable to create temporary vector map")); pvout = &Tmp; } if (!dir_bitmask(dir_fd, val_fd, next_start_pt, &window, pvout, ppl, out_mode)) { G_warning(_("No path at row %d, col %d"), next_start_pt->row, next_start_pt->col); } if (pvout) { Vect_build_partial(&Tmp, GV_BUILD_BASE); G_message(_("Breaking lines...")); Vect_break_lines(&Tmp, GV_LINE, NULL); Vect_copy_map_lines(&Tmp, &vout); Vect_set_release_support(&Tmp); Vect_close(&Tmp); /* temporary map is deleted automatically */ pvout = &vout; } } else { if (!dir_degree(dir_fd, val_fd, next_start_pt, &window, pvout, ppl, out_mode)) { G_warning(_("No path at row %d, col %d"), next_start_pt->row, next_start_pt->col); } } next_start_pt = next_start_pt->next; } /* raster output */ if (opt.rast->answer) { int row; if (pl.n > 1) { /* sort points */ qsort(pl.p, pl.n, sizeof(struct ppoint), cmp_pp); /* remove duplicates */ j = 1; for (i = 1; i < pl.n; i++) { if (pl.p[j - 1].row != pl.p[i].row || pl.p[j - 1].col != pl.p[i].col) { pl.p[j].row = pl.p[i].row; pl.p[j].col = pl.p[i].col; pl.p[j].value = pl.p[i].value; j++; } } pl.n = j; } if (out_mode == OUT_PID || out_mode == OUT_CNT) { CELL *out_buf; /* Output will be a cell map */ /* open a new file and allocate an output buffer */ out_id = Rast_open_c_new(out_name); out_buf = Rast_allocate_c_buf(); Rast_set_c_null_value(out_buf, window.cols); row = 0; /* build the output map */ G_message(_("Writing output raster map...")); for (i = 0; i < pl.n; i++) { while (row < pl.p[i].row) { G_percent(row, nrows, 2); Rast_put_c_row(out_id, out_buf); Rast_set_c_null_value(out_buf, window.cols); row++; } out_buf[pl.p[i].col] = pl.p[i].value; } while (row < window.rows) { G_percent(row, nrows, 2); Rast_put_c_row(out_id, out_buf); Rast_set_c_null_value(out_buf, window.cols); row++; } G_percent(1, 1, 1); G_free(out_buf); } else { /* mode = OUT_CPY or OUT_ACC */ DCELL *out_buf; /* Output could be of the same type as input */ /* open a new file and allocate an output buffer */ out_id = Rast_open_new(out_name, DCELL_TYPE); out_buf = Rast_allocate_d_buf(); Rast_set_d_null_value(out_buf, window.cols); row = 0; /* build the output map */ G_message(_("Writing output raster map...")); for (i = 0; i < pl.n; i++) { while (row < pl.p[i].row) { G_percent(row, nrows, 2); Rast_put_d_row(out_id, out_buf); Rast_set_d_null_value(out_buf, window.cols); row++; } out_buf[pl.p[i].col] = pl.p[i].value; } while (row < window.rows) { G_percent(row, nrows, 2); Rast_put_d_row(out_id, out_buf); Rast_set_d_null_value(out_buf, window.cols); row++; } G_percent(1, 1, 1); G_free(out_buf); } Rast_close(out_id); Rast_put_cell_title(out_name, "Path trace"); Rast_short_history(out_name, "raster", &history); Rast_command_history(&history); Rast_write_history(out_name, &history); } /* vector output */ if (opt.vect->answer) { Vect_build(&vout); Vect_close(&vout); } /* close files and free buffers */ close(dir_fd); unlink(tempfile2); if (opt.val->answer && opt.rast->answer) { close(val_fd); unlink(tempfile1); } G_done_msg(" "); exit(EXIT_SUCCESS); } void pl_add(struct point_list *pl, struct ppoint *p) { if (pl->n == pl->nalloc) { pl->nalloc += POINTS_INCREMENT; pl->p = (struct ppoint *)G_realloc(pl->p, pl->nalloc * sizeof(struct ppoint)); if (!pl->p) G_fatal_error(_("Unable to increase point list")); } pl->p[pl->n] = *p; pl->n++; } /* Jenson Domingue 1988, r.fill.dir: * bitmask encoded directions CW from North * value log2(value) direction * 1 0 NE * 2 1 E * 4 2 SE * 8 3 S * 16 4 SW * 32 5 W * 64 6 NW * 128 7 N * 256 8 NEN * 512 9 ENE * 1024 10 ESE * 2048 11 SES * 4096 12 SWS * 8192 13 WSW * 16384 14 WNW * 32768 15 NWN * * * 8 neighbours * 64 128 1 * 32 x 2 * 16 8 4 * * log2(value) CW from North starting at NE * expanded for Knight's move * * 15 8 * 14 6 7 0 9 * 5 X 1 * 13 4 3 2 10 * 12 11 */ int dir_bitmask(int dir_fd, int val_fd, struct point *startp, struct Cell_head *window, struct Map_info *Out, struct point_list *pl, int out_mode) { int go = 1, next_row, next_col, next_dir; int npoints; int dir_row = -1, val_row = -1; CELL direction; CELL *dir_buf; DCELL *val_buf; struct spoint *stackp, *newp; struct pavlrc_table *visited; struct pavlrc ngbr_rc; struct ppoint pp; int is_stack; int cur_dir, i, npaths; struct line_pnts *Points; struct line_cats *Cats; double x, y; double value; int col_offset[16] = { 1, 1, 1, 0, -1, -1, -1, 0, 1, 2, 2, 1, -1, -2, -2, -1 }; int row_offset[16] = { -1, 0, 1, 1, 1, 0, -1, -1, -2, -1, 1, 2, 2, 1, -1, -2 }; dir_buf = Rast_allocate_c_buf(); val_buf = NULL; stackp = (struct spoint *)G_malloc(sizeof(struct spoint)); stackp->row = startp->row; stackp->col = startp->col; stackp->value = startp->value; value = startp->value; stackp->dir = -1; stackp->next = NULL; visited = pavlrc_create(NULL); ngbr_rc.row = stackp->row; ngbr_rc.col = stackp->col; pavlrc_insert(visited, &ngbr_rc); if (Out) { Points = Vect_new_line_struct(); Cats = Vect_new_cats_struct(); Vect_cat_set(Cats, 1, (int)stackp->value); } if (pl) { if (out_mode == OUT_CNT) value = 1; else if (out_mode == OUT_CPY || out_mode == OUT_ACC) { /* read input raster */ val_buf = Rast_allocate_d_buf(); if (val_row != stackp->row) { lseek(val_fd, (off_t) stackp->row * window->cols * sizeof(DCELL), SEEK_SET); if (read(val_fd, val_buf, window->cols * sizeof(DCELL)) != window->cols * sizeof(DCELL)) { G_fatal_error(_("Unable to read from temp file")); } val_row = stackp->row; } value = val_buf[stackp->col]; } stackp->value = value; } /* multi-path tracing: * multiple paths from the same starting point */ npoints = 0; while (stackp) { is_stack = 1; stackp->dir++; next_row = stackp->row; next_col = stackp->col; value = stackp->value; /* begin loop */ go = 1; while (go) { go = 0; /* find direction from this point */ if (dir_row != next_row) { lseek(dir_fd, (off_t) next_row * window->cols * sizeof(CELL), SEEK_SET); if (read(dir_fd, dir_buf, window->cols * sizeof(CELL)) != window->cols * sizeof(CELL)) { G_fatal_error(_("Unable to read from temp file")); } dir_row = next_row; } direction = *(dir_buf + next_col); if (direction <= 0 || Rast_is_c_null_value(&direction)) { /* end of current path: write out the result */ if (Out && Points->n_points > 1) { Vect_write_line(Out, GV_LINE, Points, Cats); } /* leave this loop */ if (is_stack) { newp = stackp->next; G_free(stackp); stackp = newp; } break; } /* start direction */ cur_dir = 0; if (is_stack) cur_dir = stackp->dir; /* count paths going from current point and * get next direction as log2(direction) */ next_dir = -1; npaths = 0; for (i = cur_dir; i < 16; i++) { if ((direction & (1 << i)) != 0) { npaths++; if (next_dir < 0) next_dir = i; } } if (is_stack) { if (npaths > 0) { stackp->dir = next_dir; /* start a new path */ if (Out) { Vect_reset_line(Points); x = window->west + (next_col + 0.5) * window->ew_res; y = window->north - (next_row + 0.5) * window->ns_res; Vect_append_point(Points, x, y, 0.0); } if (pl) { value = stackp->value; pp.row = next_row; pp.col = next_col; pp.value = value; pl_add(pl, &pp); } npoints++; } else { /* stack point without path ? */ if (stackp->dir == 0) G_warning(_("No path from row %d, col %d"), stackp->row, stackp->col); /* drop this point from the stack */ G_debug(1, "drop point from stack"); newp = stackp->next; G_free(stackp); stackp = newp; break; } } else { /* not a stack point */ if (npaths == 0) { /* should not happen */ G_fatal_error(_("Invalid direction %d"), direction); } if (npaths > 1) { /* write out the result */ if (Out && Points->n_points > 1) { Vect_write_line(Out, GV_LINE, Points, Cats); } ngbr_rc.row = next_row; ngbr_rc.col = next_col; /* add it to the stack and leave this loop */ G_debug(1, "add point to stack: row %d, col %d, dir %d", next_row, next_col, next_dir); newp = (struct spoint *)G_malloc(sizeof(struct spoint)); newp->row = next_row; newp->col = next_col; newp->dir = next_dir - 1; newp->value = value; newp->next = stackp; stackp = newp; break; } } is_stack = 0; /* identify next downstream cell */ next_row += row_offset[next_dir]; next_col += col_offset[next_dir]; G_debug(1, "next cell at row %d, col %d", next_row, next_col); if (next_col >= 0 && next_col < window->cols && next_row >= 0 && next_row < window->rows) { /* add point */ if (Out) { x = window->west + (next_col + 0.5) * window->ew_res; y = window->north - (next_row + 0.5) * window->ns_res; Vect_append_point(Points, x, y, 0.0); } if (pl) { if (out_mode == OUT_CNT) value += 1; else if (out_mode == OUT_CPY || out_mode == OUT_ACC) { /* read input raster */ if (val_row != next_row) { lseek(val_fd, (off_t) next_row * window->cols * sizeof(DCELL), SEEK_SET); if (read(val_fd, val_buf, window->cols * sizeof(DCELL)) != window->cols * sizeof(DCELL)) { G_fatal_error(_("Unable to read from temp file")); } val_row = next_row; } if (out_mode == OUT_CPY) value = val_buf[next_col]; else value += val_buf[next_col]; } pp.row = next_row; pp.col = next_col; pp.value = value; pl_add(pl, &pp); } /* avoid circular paths * avoid tracing the same path segment several times: * a path can split and later on merge again, * then split again, merge again * path segments are identical from every merge point * to the next split point */ ngbr_rc.row = next_row; ngbr_rc.col = next_col; if (!pavlrc_insert(visited, &ngbr_rc)) { go = 1; npoints++; } else { /* reached a merge point */ /* write out the result */ if (Out && Points->n_points > 1) { Vect_write_line(Out, GV_LINE, Points, Cats); } } } else { G_warning(_("Path is leaving the current region")); } } /* end while */ } pavlrc_destroy(visited); G_free(dir_buf); if (val_buf) G_free(val_buf); if (Out) { Vect_destroy_line_struct(Points); Vect_destroy_cats_struct(Cats); } return (npoints > 1); } int dir_degree(int dir_fd, int val_fd, struct point *startp, struct Cell_head *window, struct Map_info *Out, struct point_list *pl, int out_mode) { /* * The idea is that each cell of the direction surface has a value representing * the direction towards the next cell in the path. The direction is read from * the input raster, and a simple case/switch is used to determine which cell to * read next. This is repeated via a while loop until a null direction is found. * * directions are degrees CCW from East with East = 360 */ int neighbour, next_row, next_col, go = 1; int npoints; int dir_row = -1, val_row = -1; DCELL direction; DCELL *dir_buf; DCELL *val_buf; double x, y; double value; struct ppoint pp; struct line_pnts *Points; struct line_cats *Cats; dir_buf = Rast_allocate_d_buf(); val_buf = NULL; next_row = startp->row; next_col = startp->col; value = startp->value; if (Out) { Points = Vect_new_line_struct(); Cats = Vect_new_cats_struct(); Vect_cat_set(Cats, 1, (int)value); x = window->west + (next_col + 0.5) * window->ew_res; y = window->north - (next_row + 0.5) * window->ns_res; Vect_append_point(Points, x, y, 0.0); } if (pl) { if (out_mode == OUT_CNT) value = 1; else if (out_mode == OUT_CPY || out_mode == OUT_ACC) { /* read input raster */ val_buf = Rast_allocate_d_buf(); if (val_row != next_row) { lseek(val_fd, (off_t) next_row * window->cols * sizeof(DCELL), SEEK_SET); if (read(val_fd, val_buf, window->cols * sizeof(DCELL)) != window->cols * sizeof(DCELL)) { G_fatal_error(_("Unable to read from temp file")); } val_row = next_row; } value = val_buf[next_col]; } pp.row = next_row; pp.col = next_col; pp.value = value; pl_add(pl, &pp); } npoints = 1; while (go) { go = 0; /* Directional algorithm * 1) read cell direction * 2) shift to cell in that direction */ /* find the direction recorded at row,col */ if (dir_row != next_row) { lseek(dir_fd, (off_t) next_row * window->cols * sizeof(DCELL), SEEK_SET); if (read(dir_fd, dir_buf, window->cols * sizeof(DCELL)) != window->cols * sizeof(DCELL)) { G_fatal_error(_("Unable to read from temp file")); } dir_row = next_row; } direction = *(dir_buf + next_col); neighbour = 0; if (!Rast_is_d_null_value(&direction)) { neighbour = direction * 10; G_debug(2, "direction read: %lf, neighbour found: %i", direction, neighbour); } switch (neighbour) { case 225: /* ENE */ next_row -= 1; next_col += 2; break; case 450: /* NE */ next_row -= 1; next_col += 1; break; case 675: /* NNE */ next_row -= 2; next_col += 1; break; case 900: /* N */ next_row -= 1; break; case 1125: /* NNW */ next_row -= 2; next_col -= 1; break; case 1350: /* NW */ next_col -= 1; next_row -= 1; break; case 1575: /* WNW */ next_col -= 2; next_row -= 1; break; case 1800: /* W*/ next_col -= 1; break; case 2025: /* WSW */ next_row += 1; next_col -= 2; break; case 2250: /* SW */ next_row += 1; next_col -= 1; break; case 2475: /* SSW */ next_row += 2; next_col -= 1; break; case 2700: /* S */ next_row += 1; break; case 2925: /* SSE */ next_row += 2; next_col += 1; break; case 3150: /* SE */ next_row += 1; next_col += 1; break; case 3375: /* ESE */ next_row += 1; next_col += 2; break; case 3600: /* E */ next_col += 1; break; default: /* end of path */ next_row = -1; next_col = -1; break; } /* end switch/case */ if (next_col >= 0 && next_col < window->cols && next_row >= 0 && next_row < window->rows) { if (Out) { x = window->west + (next_col + 0.5) * window->ew_res; y = window->north - (next_row + 0.5) * window->ns_res; Vect_append_point(Points, x, y, 0.0); } if (pl) { if (out_mode == OUT_CNT) value += 1; else if (out_mode == OUT_CPY || out_mode == OUT_ACC) { /* read input raster */ if (val_row != next_row) { lseek(val_fd, (off_t) next_row * window->cols * sizeof(DCELL), SEEK_SET); if (read(val_fd, val_buf, window->cols * sizeof(DCELL)) != window->cols * sizeof(DCELL)) { G_fatal_error(_("Unable to read from temp file")); } val_row = next_row; } if (out_mode == OUT_CPY) value = val_buf[next_col]; else value += val_buf[next_col]; } pp.row = next_row; pp.col = next_col; pp.value = value; pl_add(pl, &pp); } go = 1; npoints++; } } /* end while */ if (Out && Points->n_points > 1) { Vect_write_line(Out, GV_LINE, Points, Cats); } G_free(dir_buf); if (val_buf) G_free(val_buf); if (Out) { Vect_destroy_line_struct(Points); Vect_destroy_cats_struct(Cats); } return (npoints > 1); }