diff --git a/raster/r.carve/Makefile b/raster/r.carve/Makefile index 49b6ee9670d..f452758a600 100644 --- a/raster/r.carve/Makefile +++ b/raster/r.carve/Makefile @@ -2,13 +2,12 @@ MODULE_TOPDIR = ../.. PGM = r.carve -LIBES = $(VECTORLIB) $(BITMAPLIB) $(DIG2LIB) $(RASTERLIB) $(GISLIB) +LIBES = $(VECTORLIB) $(BITMAPLIB) $(DIG2LIB) $(RASTERLIB) $(GISLIB) $(DBMILIB) DEPENDENCIES = $(VECTORDEP) $(BITMAPDEP) $(DIG2DEP) $(RASTERDEP) $(GISDEP) EXTRA_INC = $(VECT_INC) EXTRA_CFLAGS = $(VECT_CFLAGS) - -include $(MODULE_TOPDIR)/include/Make/Module.make -default: cmd +include $(MODULE_TOPDIR)/include/Make/Module.make +default: cmd diff --git a/raster/r.carve/enforce.h b/raster/r.carve/enforce.h index 46d71804c09..fd448140b9a 100644 --- a/raster/r.carve/enforce.h +++ b/raster/r.carve/enforce.h @@ -1,10 +1,11 @@ - /**************************************************************************** * * MODULE: r.carve * * AUTHOR(S): Original author Bill Brown, UIUC GIS Laboratory * Brad Douglas + * Tomas Zigo (adding the option + * to read width, depth values from vector map table columns) * * PURPOSE: Takes vector stream data, converts it to 3D raster and * subtracts a specified depth @@ -49,18 +50,62 @@ typedef struct struct parms { - struct Option *inrast, *invect, *outrast, *outvect; + struct Option *inrast, *invect, *outrast, *outvect, *width_col, + *depth_col, *field; RASTER_MAP_TYPE raster_type; double swidth, sdepth; int wrap, noflat; }; +struct sql_statement +{ + dbString *sql; + int ncats; + struct vect_id_cat_map *id_cat_map; +}; + +struct vect_id_cat_map +{ + int id; + int cat; +}; + +struct ptr +{ + enum Type { + P_INT, + P_DOUBLE, + P_CHAR, + P_DBSTRING, + P_VECT_ID_CAT_MAP, + } type; + union { + int *p_int; + double *p_double; + char *p_char; + dbString *p_dbString; + struct vect_id_cat_map *p_vect_id_cat_map; + }; +}; + +typedef enum +{ + WIDTH, + DEPTH, +} value_type; /* enforce_ds.c */ -extern int enforce_downstream(int /*infd */ , int /*outfd */ , - struct Map_info * /*Map */ , - struct Map_info * /*outMap */ , - struct parms * /* parm */ ); +extern void enforce_downstream(int /*infd */ , int /*outfd */ , + struct Map_info * /*Map */ , + struct Map_info * /*outMap */ , + struct parms * /* parm */, + struct field_info * /* Fi */, + int * /* width_col_posw */, + int * /* depth_col_pos */, + char *[2] /* columns[2] */, + dbDriver * /* driver */); +extern void adjust_swidth(struct Cell_head *win, double *value); +extern void adjust_sdepth(double *value); /* lobf.c */ extern Point2 *pg_getpoints(PointGrp *); @@ -75,6 +120,7 @@ void *write_raster(void *, const int, const RASTER_MAP_TYPE); /* support.c */ extern void update_rast_history(struct parms *); +extern void check_mem_alloc(struct ptr *pointer); /* vect.c */ extern int open_new_vect(struct Map_info *, char *); diff --git a/raster/r.carve/enforce_ds.c b/raster/r.carve/enforce_ds.c index 50c0b0a7dc0..dd46d474632 100644 --- a/raster/r.carve/enforce_ds.c +++ b/raster/r.carve/enforce_ds.c @@ -1,10 +1,11 @@ - /**************************************************************************** * * MODULE: r.carve * * AUTHOR(S): Original author Bill Brown, UIUC GIS Laboratory * Brad Douglas + * Tomas Zigo (adding the option + * to read width, depth values from vector map table columns) * * PURPOSE: Takes vector stream data, converts it to 3D raster and * subtracts a specified depth @@ -17,6 +18,7 @@ * ****************************************************************************/ +#include #include #include #include @@ -31,179 +33,308 @@ # define MAX(a,b) ((a>b) ? a : b) #endif - /* function prototypes */ static void clear_bitmap(struct BM *bm); -static int process_line(struct Map_info *Map, struct Map_info *outMap, - void *rbuf, const int line, const struct parms *parm); +void process_line(struct Map_info *Map, struct Map_info *outMap, + const struct parms *parm, void *rbuf, const int *line); static void traverse_line_flat(Point2 * pgpts, const int pt, const int npts); static void traverse_line_noflat(Point2 * pgpts, const double depth, - const int pt, const int npts); + const int pt, const int npts); static void set_min_point(void *buf, const int col, const int row, - const double elev, const double depth, - const RASTER_MAP_TYPE rtype); + const double elev, const double depth, + const RASTER_MAP_TYPE rtype); static double lowest_cell_near_point(void *data, const RASTER_MAP_TYPE rtype, - const double px, const double py, - const double rad); + const double px, const double py, + const double rad); static void process_line_segment(const int npts, void *rbuf, Point2 * pgxypts, - Point2 * pgpts, struct BM *bm, - struct Map_info *outMap, - const struct parms *parm); - - -/****************************************************************** - * Returns 0 on success, -1 on error, 1 on warning, writing message - * to errbuf for -1 or 1 */ - -int enforce_downstream(int infd, int outfd, - struct Map_info *Map, struct Map_info *outMap, - struct parms *parm) + Point2 * pgpts, struct BM *bm, + struct Map_info *outMap, + const struct parms *parm); +double get_value(unsigned short int *ctype, dbColumn *col); +void set_value(dbColumn *col, unsigned short int *ctype, char *answer, + double *parm, double *def_value, dbTable *table, + struct Cell_head *wind, value_type type); +struct sql_statement create_select_sql_statement(struct Map_info *Map, + struct field_info *Fi, + struct boxlist *box_list, + char *columns[2], + unsigned short int *field, + char *keycol); + + +void enforce_downstream(int infd, int outfd, + struct Map_info *Map, struct Map_info *outMap, + struct parms *parm, struct field_info *Fi, + int *width_col_pos, int *depth_col_pos, + char *columns[2], dbDriver *driver) +/* + * Function: enforce_downstream + * ------------------------- + * Enforce downstream + * + * infd: input raster elevation map + * outfd: output raster elevation map with carved streams + * Map: input streams vector map + * outMap: output vector map for adjusted stream points + * parm: module parameters + * Fi: layer (old: field) information + * width_col_pos: width column position in the columns array + * depth_col_pos: depth column position in the columns array + * columns: array of width and depth column names + * driver: db driver + * + */ { struct Cell_head wind; - int retval = 0; - int line, nlines; void *rbuf = NULL; - - /* width used here is actually distance to center of stream */ - parm->swidth /= 2; + struct bound_box box; + struct boxlist *box_list; + int more; + unsigned int c; + unsigned short int field, width_col_type, depth_col_type; + // Width used here is actually distance to center of stream + unsigned short int const distance = 2; + double *def_depth = NULL, *def_width = NULL; + dbCursor cursor; + dbTable *table = NULL; + dbColumn *width_col = NULL, *depth_col = NULL, *cat_col = NULL; + struct field_info *finfo = NULL; + struct sql_statement sql; G_get_window(&wind); Vect_set_constraint_region(Map, wind.north, wind.south, wind.east, - wind.west, wind.top, wind.bottom); + wind.west, wind.top, wind.bottom); /* allocate and clear memory for entire raster */ - rbuf = - G_calloc(Rast_window_rows() * Rast_window_cols(), - Rast_cell_size(parm->raster_type)); + rbuf = G_calloc(Rast_window_rows() * Rast_window_cols(), + Rast_cell_size(parm->raster_type)); /* first read whole elevation file into buf */ read_raster(rbuf, infd, parm->raster_type); G_message(_("Processing lines... ")); - nlines = Vect_get_num_lines(Map); - for (line = 1; line <= nlines; line++) - retval = process_line(Map, outMap, rbuf, line, parm); + box_list = Vect_new_boxlist(0); + + box.N = wind.north; + box.E = wind.east; + box.S = wind.south; + box.W = wind.west; + box.B = wind.bottom; + box.T = wind.top; + + Vect_select_lines_by_box(Map, &box, GV_LINE, box_list); + field = Vect_get_field_number(Map, parm->field->answer); + + // Get key col + finfo = Vect_get_field(Map, field); + + if (parm->width_col->answer || parm->depth_col->answer) + { + sql = create_select_sql_statement(Map, Fi, box_list, columns, + &field, finfo->key); + + if (db_open_select_cursor(driver, sql.sql, + &cursor, DB_SEQUENTIAL) != DB_OK) + G_fatal_error(_("Unable to open select cursor")); + + def_width = G_malloc(sizeof(double)); + def_depth = G_malloc(sizeof(double)); + memcpy(def_width, &parm->swidth, sizeof(double)); + memcpy(def_depth, &parm->sdepth, sizeof(double)); + + table = db_get_cursor_table(&cursor); + G_debug(1, "Default width: %.2f, depth: %.2f", + *def_width / distance, *def_depth); + + cat_col = db_get_table_column_by_name(table, + finfo->key); + if (parm->width_col->answer) + { + width_col = db_get_table_column_by_name(table, + columns[*width_col_pos]); + width_col_type = db_get_column_sqltype(width_col); + } + if (parm->depth_col->answer) + { + depth_col = db_get_table_column_by_name(table, + columns[*depth_col_pos]); + depth_col_type = db_get_column_sqltype(depth_col); + } + + while (1) + { + if (db_fetch(&cursor, DB_NEXT, &more) != DB_OK) + G_fatal_error(_("Unable to fetch data from table <%s>"), + Fi->table); + + if (!more) + break; + + int cat = db_get_value_int(db_get_column_value(cat_col)); + + // Set width + set_value(width_col, &width_col_type, parm->width_col->answer, + &parm->swidth, def_width, table, &wind, WIDTH); + parm->swidth /= distance; + // Set depth + set_value(depth_col, &depth_col_type, parm->depth_col->answer, + &parm->sdepth, def_depth, table, &wind, DEPTH); + + for (c = 0; c < sql.ncats; c++) + { + // Line cat has an entry in the db table + if (cat == sql.id_cat_map[c].cat) + { + G_debug(3, "Process line with id: %d, cat: %d, " + "width: %.2f, depth: %.2f", + sql.id_cat_map[c].id, cat, parm->swidth, + parm->sdepth); + + process_line(Map, outMap, parm, rbuf, + &sql.id_cat_map[c].id); + } + } + } + db_free_column(cat_col); + if (parm->width_col->answer) + { + db_free_column(width_col); + } + if (parm->depth_col->answer) + { + db_free_column(depth_col); + } + if (db_close_cursor(&cursor) != DB_OK) + G_fatal_error(_("Unable to close select cursor")); + db_close_database_shutdown_driver(driver); + db_free_string(sql.sql); + G_free(sql.id_cat_map); + G_free(def_width); + G_free(def_depth); + } /* write output raster map */ write_raster(rbuf, outfd, parm->raster_type); G_free(rbuf); - - return retval; } - -static int process_line(struct Map_info *Map, struct Map_info *outMap, - void *rbuf, const int line, const struct parms *parm) +void process_line(struct Map_info *Map, struct Map_info *outMap, + const struct parms *parm, void *rbuf, const int *line) +/* + * Function: process_line + * ------------------------- + * Process line + * + * Map: input streams vector map + * outfd: output raster elevation map with carved streams + * outMap: output vector map for adjusted stream points + * parm: module parameters + * rbuf: raster elevation map buffer + * line: streams vector map line id + * + */ { - int i, retval = 0; - int do_warn = 0; - int npts = 0; - int in_out = 0; - int first_in = -1; + int i, do_warn = 0, first_in = -1, in_out = 0, npts = 0; double totdist = 0.; struct Cell_head wind; static struct line_pnts *points = NULL; static struct line_cats *cats = NULL; static struct BM *bm = NULL; - Point2 *pgpts, *pgxypts; + Point2 *pgpts = NULL, *pgxypts = NULL; PointGrp pg; - PointGrp pgxy; /* copy any points in region to this one */ + PointGrp pgxy; /* copy any points in region to this one */ G_get_window(&wind); if (!points) - points = Vect_new_line_struct(); + points = Vect_new_line_struct(); if (!cats) - cats = Vect_new_cats_struct(); + cats = Vect_new_cats_struct(); - if (!(Vect_read_line(Map, points, cats, line) & GV_LINE)) - return 0; + if (!(Vect_read_line(Map, points, cats, *line) & GV_LINE)) + return; if (!bm) - bm = BM_create(Rast_window_cols(), Rast_window_rows()); + bm = BM_create(Rast_window_cols(), Rast_window_rows()); clear_bitmap(bm); pg_init(&pg); pg_init(&pgxy); - G_percent(line, Vect_get_num_lines(Map), 10); + G_percent(*line, Vect_get_num_lines(Map), 10); for (i = 0; i < points->n_points; i++) { - Point2 pt, ptxy; - double elev; - int row = Rast_northing_to_row(points->y[i], &wind); - int col = Rast_easting_to_col(points->x[i], &wind); - - /* rough clipping */ - if (row < 0 || row > Rast_window_rows() - 1 || - col < 0 || col > Rast_window_cols() - 1) { - if (first_in != -1) - in_out = 1; - - G_debug(1, "outside region - row:%d col:%d", row, col); - - continue; - } - - if (first_in < 0) - first_in = i; - else if (in_out) - do_warn = 1; - - elev = lowest_cell_near_point(rbuf, parm->raster_type, points->x[i], - points->y[i], parm->swidth); - - ptxy[0] = points->x[i]; - ptxy[1] = points->y[i]; - pt[1] = elev; - - /* get distance from this point to previous point */ - if (i) - totdist += G_distance(points->x[i - 1], points->y[i - 1], - points->x[i], points->y[i]); - - pt[0] = totdist; - pg_addpt(&pg, pt); - pg_addpt(&pgxy, ptxy); - npts++; + Point2 pt, ptxy; + double elev; + int row = Rast_northing_to_row(points->y[i], &wind); + int col = Rast_easting_to_col(points->x[i], &wind); + + /* rough clipping */ + if (row < 0 || row > Rast_window_rows() - 1 || + col < 0 || col > Rast_window_cols() - 1) { + if (first_in != -1) + in_out = 1; + + G_debug(1, "outside region - row:%d col:%d", row, col); + + continue; + } + + if (first_in < 0) + first_in = i; + else if (in_out) + do_warn = 1; + + elev = lowest_cell_near_point(rbuf, parm->raster_type, points->x[i], + points->y[i], parm->swidth); + + ptxy[0] = points->x[i]; + ptxy[1] = points->y[i]; + pt[1] = elev; + + /* get distance from this point to previous point */ + if (i) + totdist += G_distance(points->x[i - 1], points->y[i - 1], + points->x[i], points->y[i]); + + pt[0] = totdist; + pg_addpt(&pg, pt); + pg_addpt(&pgxy, ptxy); + npts++; } if (do_warn) { - G_warning(_("Vect runs out of region and re-enters - " - "this case is not yet implemented.")); - retval = 1; + G_warning(_("Vect runs out of region and re-enters - " + "this case is not yet implemented.")); } /* now check to see if points go downslope(inorder) or upslope */ if (pg_y_from_x(&pg, 0.0) > pg_y_from_x(&pg, totdist)) { - pgpts = pg_getpoints(&pg); - pgxypts = pg_getpoints(&pgxy); + pgpts = pg_getpoints(&pg); + pgxypts = pg_getpoints(&pgxy); } else { - /* pgpts is now high to low */ - pgpts = pg_getpoints_reversed(&pg); + /* pgpts is now high to low */ + pgpts = pg_getpoints_reversed(&pg); - for (i = 0; i < npts; i++) - pgpts[i][0] = totdist - pgpts[i][0]; + for (i = 0; i < npts; i++) + pgpts[i][0] = totdist - pgpts[i][0]; - pgxypts = pg_getpoints_reversed(&pgxy); + pgxypts = pg_getpoints_reversed(&pgxy); } for (i = 0; i < (npts - 1); i++) { - if (parm->noflat) - /* make sure there are no flat segments in line */ - traverse_line_noflat(pgpts, parm->sdepth, i, npts); - else - /* ok to have flat segments in line */ - traverse_line_flat(pgpts, i, npts); + if (parm->noflat) + /* make sure there are no flat segments in line */ + traverse_line_noflat(pgpts, parm->sdepth, i, npts); + else + /* ok to have flat segments in line */ + traverse_line_flat(pgpts, i, npts); } - process_line_segment(npts, rbuf, pgxypts, pgpts, bm, outMap, parm); - - return retval; } @@ -212,8 +343,8 @@ static void clear_bitmap(struct BM *bm) int i, j; for (i = 0; i < Rast_window_rows(); i++) - for (j = 0; j < Rast_window_cols(); j++) - BM_set(bm, i, j, 0); + for (j = 0; j < Rast_window_cols(); j++) + BM_set(bm, i, j, 0); } @@ -222,91 +353,91 @@ static void traverse_line_flat(Point2 * pgpts, const int pt, const int npts) int j, k; if (pgpts[pt + 1][1] <= pgpts[pt][1]) - return; + return; for (j = (pt + 2); j < npts; j++) - if (pgpts[j][1] <= pgpts[pt][1]) - break; + if (pgpts[j][1] <= pgpts[pt][1]) + break; if (j == npts) { - /* if we got to the end, level it out */ - for (j = (pt + 1); j < npts; j++) - pgpts[j][1] = pgpts[pt][1]; + /* if we got to the end, level it out */ + for (j = (pt + 1); j < npts; j++) + pgpts[j][1] = pgpts[pt][1]; } else { - /* linear interp between point pt and the next < */ - for (k = (pt + 1); k < j; k++) - pgpts[k][1] = LINTERP(pgpts[j][1], pgpts[pt][1], - (pgpts[j][0] - pgpts[k][0]) / - (pgpts[j][0] - pgpts[pt][0])); + /* linear interp between point pt and the next < */ + for (k = (pt + 1); k < j; k++) + pgpts[k][1] = LINTERP(pgpts[j][1], pgpts[pt][1], + (pgpts[j][0] - pgpts[k][0]) / + (pgpts[j][0] - pgpts[pt][0])); } } static void traverse_line_noflat(Point2 * pgpts, const double depth, - const int pt, const int npts) + const int pt, const int npts) { int j, k; if (pgpts[pt + 1][1] < pgpts[pt][1]) - return; + return; for (j = (pt + 2); j < npts; j++) - if (pgpts[j][1] < pgpts[pt][1]) - break; + if (pgpts[j][1] < pgpts[pt][1]) + break; if (j == npts) { - /* if we got to the end, lower end by depth OR .01 */ - --j; - pgpts[j][1] = pgpts[pt][1] - (depth > 0 ? depth : 0.01); + /* if we got to the end, lower end by depth OR .01 */ + --j; + pgpts[j][1] = pgpts[pt][1] - (depth > 0 ? depth : 0.01); } /* linear interp between point pt and the next < */ for (k = (pt + 1); k < j; k++) - pgpts[k][1] = LINTERP(pgpts[j][1], pgpts[pt][1], - (pgpts[j][0] - pgpts[k][0]) / - (pgpts[j][0] - pgpts[pt][0])); + pgpts[k][1] = LINTERP(pgpts[j][1], pgpts[pt][1], + (pgpts[j][0] - pgpts[k][0]) / + (pgpts[j][0] - pgpts[pt][0])); } /* sets value for a cell */ static void set_min_point(void *data, const int col, const int row, - const double elev, const double depth, - const RASTER_MAP_TYPE rtype) + const double elev, const double depth, + const RASTER_MAP_TYPE rtype) { switch (rtype) { case CELL_TYPE: - { - CELL *cbuf = data; + { + CELL *cbuf = data; - cbuf[row * Rast_window_cols() + col] = - MIN(cbuf[row * Rast_window_cols() + col], elev) - (int)depth; - } - break; + cbuf[row * Rast_window_cols() + col] = + MIN(cbuf[row * Rast_window_cols() + col], elev) - (int)depth; + } + break; case FCELL_TYPE: - { - FCELL *fbuf = data; + { + FCELL *fbuf = data; - fbuf[row * Rast_window_cols() + col] = - MIN(fbuf[row * Rast_window_cols() + col], elev) - depth; - } - break; + fbuf[row * Rast_window_cols() + col] = + MIN(fbuf[row * Rast_window_cols() + col], elev) - depth; + } + break; case DCELL_TYPE: - { - DCELL *dbuf = data; + { + DCELL *dbuf = data; - dbuf[row * Rast_window_cols() + col] = - MIN(dbuf[row * Rast_window_cols() + col], elev) - depth; - } - break; + dbuf[row * Rast_window_cols() + col] = + MIN(dbuf[row * Rast_window_cols() + col], elev) - depth; + } + break; } } /* returns the lowest value cell within radius rad of px, py */ static double lowest_cell_near_point(void *data, const RASTER_MAP_TYPE rtype, - const double px, const double py, - const double rad) + const double px, const double py, + const double rad) { int r, row, col, row1, row2, col1, col2, rowoff, coloff; int rastcols, rastrows; @@ -334,88 +465,88 @@ static double lowest_cell_near_point(void *data, const RASTER_MAP_TYPE rtype, switch (rtype) { case CELL_TYPE: - { - CELL *cbuf = data; + { + CELL *cbuf = data; - if (!(Rast_is_c_null_value(&cbuf[row1 * rastcols + col1]))) - min = cbuf[row1 * rastcols + col1]; - } - break; + if (!(Rast_is_c_null_value(&cbuf[row1 * rastcols + col1]))) + min = cbuf[row1 * rastcols + col1]; + } + break; case FCELL_TYPE: - { - FCELL *fbuf = data; + { + FCELL *fbuf = data; - if (!(Rast_is_f_null_value(&fbuf[row1 * rastcols + col1]))) - min = fbuf[row1 * rastcols + col1]; - } - break; + if (!(Rast_is_f_null_value(&fbuf[row1 * rastcols + col1]))) + min = fbuf[row1 * rastcols + col1]; + } + break; case DCELL_TYPE: - { - DCELL *dbuf = data; + { + DCELL *dbuf = data; - if (!(Rast_is_d_null_value(&dbuf[row1 * rastcols + col1]))) - min = dbuf[row1 * rastcols + col1]; - } - break; + if (!(Rast_is_d_null_value(&dbuf[row1 * rastcols + col1]))) + min = dbuf[row1 * rastcols + col1]; + } + break; } for (r = row1; r < row2; r++) { - double cy = Rast_row_to_northing(r + 0.5, &wind); - int c; - - for (c = col1; c < col2; c++) { - double cx = Rast_col_to_easting(c + 0.5, &wind); - - if (G_distance(px, py, cx, cy) <= SQR(rad)) { - switch (rtype) { - case CELL_TYPE: - { - CELL *cbuf = data; - - if (Rast_is_d_null_value(&min)) { - if (!(Rast_is_c_null_value(&cbuf[r * rastcols + c]))) - min = cbuf[r * rastcols + c]; - } - else { - if (!(Rast_is_c_null_value(&cbuf[r * rastcols + c]))) - if (cbuf[r * rastcols + c] < min) - min = cbuf[r * rastcols + c]; - } - } - break; - case FCELL_TYPE: - { - FCELL *fbuf = data; - - if (Rast_is_d_null_value(&min)) { - if (!(Rast_is_f_null_value(&fbuf[r * rastcols + c]))) - min = fbuf[r * rastcols + c]; - } - else { - if (!(Rast_is_f_null_value(&fbuf[r * rastcols + c]))) - if (fbuf[r * rastcols + c] < min) - min = fbuf[r * rastcols + c]; - } - } - break; - case DCELL_TYPE: - { - DCELL *dbuf = data; - - if (Rast_is_d_null_value(&min)) { - if (!(Rast_is_d_null_value(&dbuf[r * rastcols + c]))) - min = dbuf[r * rastcols + c]; - } - else { - if (!(Rast_is_d_null_value(&dbuf[r * rastcols + c]))) - if (dbuf[r * rastcols + c] < min) - min = dbuf[r * rastcols + c]; - } - } - break; - } - } - } + double cy = Rast_row_to_northing(r + 0.5, &wind); + int c; + + for (c = col1; c < col2; c++) { + double cx = Rast_col_to_easting(c + 0.5, &wind); + + if (G_distance(px, py, cx, cy) <= SQR(rad)) { + switch (rtype) { + case CELL_TYPE: + { + CELL *cbuf = data; + + if (Rast_is_d_null_value(&min)) { + if (!(Rast_is_c_null_value(&cbuf[r * rastcols + c]))) + min = cbuf[r * rastcols + c]; + } + else { + if (!(Rast_is_c_null_value(&cbuf[r * rastcols + c]))) + if (cbuf[r * rastcols + c] < min) + min = cbuf[r * rastcols + c]; + } + } + break; + case FCELL_TYPE: + { + FCELL *fbuf = data; + + if (Rast_is_d_null_value(&min)) { + if (!(Rast_is_f_null_value(&fbuf[r * rastcols + c]))) + min = fbuf[r * rastcols + c]; + } + else { + if (!(Rast_is_f_null_value(&fbuf[r * rastcols + c]))) + if (fbuf[r * rastcols + c] < min) + min = fbuf[r * rastcols + c]; + } + } + break; + case DCELL_TYPE: + { + DCELL *dbuf = data; + + if (Rast_is_d_null_value(&min)) { + if (!(Rast_is_d_null_value(&dbuf[r * rastcols + c]))) + min = dbuf[r * rastcols + c]; + } + else { + if (!(Rast_is_d_null_value(&dbuf[r * rastcols + c]))) + if (dbuf[r * rastcols + c] < min) + min = dbuf[r * rastcols + c]; + } + } + break; + } + } + } } G_debug(3, "min:%.2lf", min); @@ -424,15 +555,15 @@ static double lowest_cell_near_point(void *data, const RASTER_MAP_TYPE rtype, } -/* Now for each segment in the line, use distance from segment +/* Now for each segment in the line, use distance from segment * to find beginning row from northernmost point, beginning - * col from easternmost, ending row & col, then loop through + * col from easternmost, ending row & col, then loop through * bounding box and use distance from segment to emboss * new elevations */ static void process_line_segment(const int npts, void *rbuf, - Point2 * pgxypts, Point2 * pgpts, - struct BM *bm, struct Map_info *outMap, - const struct parms *parm) + Point2 * pgxypts, Point2 * pgpts, + struct BM *bm, struct Map_info *outMap, + const struct parms *parm) { int i, row1, row2, col1, col2; int prevrow, prevcol; @@ -455,69 +586,325 @@ static void process_line_segment(const int npts, void *rbuf, prevcol = Rast_easting_to_col(pgxypts[0][0], &wind); for (i = 1; i < npts; i++) { - int c, r; - - int row = Rast_northing_to_row(pgxypts[i][1], &wind); - int col = Rast_easting_to_col(pgxypts[i][0], &wind); - - /* get bounding box of line segment */ - row1 = MAX(0, MIN(row, prevrow) - rowoff); - row2 = MIN(Rast_window_rows() - 1, MAX(row, prevrow) + rowoff); - col1 = MAX(0, MIN(col, prevcol) - coloff); - col2 = MIN(Rast_window_cols() - 1, MAX(col, prevcol) + coloff); - - for (r = row1; r <= row2; r++) { - cy = Rast_row_to_northing(r + 0.5, &wind); - - for (c = col1; c <= col2; c++) { - double distance; - - cellx = Rast_col_to_easting(c + 0.5, &wind); - celly = cy; /* gets written over in distance2... */ - - /* Thought about not going past endpoints (use - * status to check) but then pieces end up missing - * from outside corners - if it goes past ends, - * should probably do some interp or will get flats. - * Here we use a bitmap and only change cells once - * on the way down */ - - distance = sqrt(dig_distance2_point_to_line(cellx, celly, 0, - pgxypts[i - 1][0], - pgxypts[i - 1][1], 0, - pgxypts[i][0], pgxypts[i][1], - 0, 0, &cellx, &celly, NULL, - NULL, NULL)); - - if (distance <= parm->swidth && !BM_get(bm, c, r)) { - double dist, elev; - - Vect_reset_line(points); - - dist = G_distance(pgxypts[i][0], pgxypts[i][1], - cellx, celly); - - elev = LINTERP(pgpts[i][1], pgpts[i - 1][1], - (dist / (pgpts[i][0] - pgpts[i - 1][0]))); - - BM_set(bm, c, r, 1); - - /* TODO - may want to use a function for the - * cross section of stream */ - set_min_point(rbuf, c, r, elev, parm->sdepth, - parm->raster_type); - - /* Add point to output vector map */ - if (parm->outvect->answer) { - Vect_append_point(points, cellx, celly, - elev - parm->sdepth); - Vect_write_line(outMap, GV_POINT, points, cats); - } - } - } - } - - prevrow = row; - prevcol = col; + int c, r; + + int row = Rast_northing_to_row(pgxypts[i][1], &wind); + int col = Rast_easting_to_col(pgxypts[i][0], &wind); + + /* get bounding box of line segment */ + row1 = MAX(0, MIN(row, prevrow) - rowoff); + row2 = MIN(Rast_window_rows() - 1, MAX(row, prevrow) + rowoff); + col1 = MAX(0, MIN(col, prevcol) - coloff); + col2 = MIN(Rast_window_cols() - 1, MAX(col, prevcol) + coloff); + + for (r = row1; r <= row2; r++) { + cy = Rast_row_to_northing(r + 0.5, &wind); + + for (c = col1; c <= col2; c++) { + double distance; + + cellx = Rast_col_to_easting(c + 0.5, &wind); + celly = cy; /* gets written over in distance2... */ + + /* Thought about not going past endpoints (use + * status to check) but then pieces end up missing + * from outside corners - if it goes past ends, + * should probably do some interp or will get flats. + * Here we use a bitmap and only change cells once + * on the way down */ + + distance = sqrt(dig_distance2_point_to_line(cellx, celly, 0, + pgxypts[i - 1][0], + pgxypts[i - 1][1], 0, + pgxypts[i][0], pgxypts[i][1], + 0, 0, &cellx, &celly, NULL, + NULL, NULL)); + + if (distance <= parm->swidth && !BM_get(bm, c, r)) { + double dist, elev; + + Vect_reset_line(points); + + dist = G_distance(pgxypts[i][0], pgxypts[i][1], + cellx, celly); + + elev = LINTERP(pgpts[i][1], pgpts[i - 1][1], + (dist / (pgpts[i][0] - pgpts[i - 1][0]))); + + BM_set(bm, c, r, 1); + + /* TODO - may want to use a function for the + * cross section of stream */ + set_min_point(rbuf, c, r, elev, parm->sdepth, + parm->raster_type); + + /* Add point to output vector map */ + if (parm->outvect->answer) { + Vect_append_point(points, cellx, celly, + elev - parm->sdepth); + Vect_write_line(outMap, GV_POINT, points, cats); + } + } + } } + + prevrow = row; + prevcol = col; + } +} + +double get_value(unsigned short int *ctype, dbColumn *col) +/* + * Function: get_value + * -------------------- + * Get value from column + * + * ctype: column type + * col: column + * + * return: column value + */ +{ + if (*ctype == DB_SQL_TYPE_INTEGER) + return (double)(db_get_value_int(db_get_column_value(col))); + else + return db_get_value_double(db_get_column_value(col)); +} + +void set_value(dbColumn *col, unsigned short int *ctype, char *answer, + double *parm, double *def_value, dbTable *table, + struct Cell_head *wind, value_type type) +/* + * Function: set_value + * -------------------- + * Set correct column (with or depth) value + * + * col: db table column + * ctype: db table column type + * answer: column parameter answer + * parm: module parameters + * def_value: default column value + * table: streams vector map table + * wind: current region settings + * type: column type (width or depth) + */ +{ + double value; + + if (answer) + { + if (!db_get_column_value(col)->isNull) + { + value = get_value(ctype, col); + switch (type) + { + case WIDTH: + adjust_swidth(wind, &value); + break; + case DEPTH: + adjust_sdepth(&value); + break; + } + *parm = value; + } + else + *parm = *def_value; + } + else + *parm = *def_value; +} + + +void adjust_swidth(struct Cell_head *win, double *value) +/* + * Function: adjust_swidth + * ------------------------ + * Adjust width value + * + * win: current region settings + * value: input width value + */ +{ + double width = 0.0; + + if (*value <= width) + { + double def_width = G_distance((win->east + win->west) / 2, + (win->north + win->south) / 2, + ((win->east + win->west) / 2) + + win->ew_res, + (win->north + win->south) / 2); + *value = def_width; + } +} + + +void adjust_sdepth(double *value) +/* + * Function: adjust_sdepth + * ----------------------- + * Adjust depth value + * + * value: input depth value + */ +{ + double def_depth = 0.0; + + if (*value < def_depth) + *value = def_depth; +} + + +struct sql_statement create_select_sql_statement(struct Map_info *Map, + struct field_info *Fi, + struct boxlist *box_list, + char *columns[2], + unsigned short int *field, + char *keycol) +/* + * Function: create_select_sql_statement + * ------------------------------------- + * Creates sql string + * + * Map: input streams vector line map + * Fi: layer (old: field) information + * box_list: list of bounding boxes with id (line streams inside region) + * columns: width and depth column names + * field: layer (old: field) number) + * keycol: key column name + * + * return ret: sql_statement struct + * sql_statement.sql: sql string + * sql.sql: number of cats + * sql.id_cat_map: line id with matching cat ( + * vect_id_cat_map struct) + * + */ +{ + int cat_buf_len, size, line, ncats, next_cat, prev_line = -2; + bool no_cat = false; + char query[DB_SQL_MAX]; + char *cat_buf = NULL, *where_buf = NULL; + dbString sql; /* value_string; */ + struct sql_statement ret; + struct ptr *p = G_malloc(sizeof(struct ptr)); + + db_init_string(&sql); + + sprintf(query, "SELECT %s, %s, %s FROM ", keycol, columns[0], + columns[1]); + db_set_string(&sql, query); + if (db_append_string(&sql, Fi->table) != DB_OK) + G_fatal_error(_("Unable to append string")); + + cat_buf_len = 0; + ncats = 0; + + // Create WHERE clausule ("WHERE keycol in (value, ....)") + for (line = 0; line < box_list->n_values; line++) + { + int cat = Vect_get_line_cat(Map, box_list->id[line], *field); + // Filter out cat = -1 + if (cat < 0) + { + no_cat = true; + prev_line = line; + + if (line == box_list->n_values - 1 && !no_cat) + { + size = snprintf(NULL, 0, "%d)", cat); + cat_buf = G_realloc(cat_buf, strlen(cat_buf) + size + 1); + p->type = P_CHAR; + p->p_char = cat_buf; + check_mem_alloc(p); + cat_buf_len += sprintf(cat_buf + cat_buf_len, "%d)", cat); + } + continue; + } + else + { + ncats += 1; + if (ncats == 1) + ret.id_cat_map = G_malloc(sizeof(struct vect_id_cat_map)); + else + ret.id_cat_map = G_realloc(ret.id_cat_map, + sizeof(struct vect_id_cat_map) * + ncats); + + p->type = P_VECT_ID_CAT_MAP; + p->p_vect_id_cat_map = ret.id_cat_map; + check_mem_alloc(p); + // Saving id with corresponding cat value + ret.id_cat_map[ncats - 1].cat = cat; + ret.id_cat_map[ncats - 1].id = box_list->id[line]; + } + if (line == 0 || (prev_line == 0 && no_cat)) + { + size = snprintf(NULL, 0, "(%d, ", cat); + cat_buf = G_malloc(size + 1); + p->type = P_CHAR; + p->p_char = cat_buf; + check_mem_alloc(p); + cat_buf_len += sprintf(cat_buf + cat_buf_len, "(%d, ", cat); + + if (prev_line == 0) + prev_line = -2; + } + else if (line > 0 && line < box_list->n_values - 1) + { + next_cat = Vect_get_line_cat(Map, box_list->id[line + 1], *field); + if (next_cat < 0 && line + 1 == box_list->n_values - 1) + { + size = snprintf(NULL, 0, "%d)", cat); + cat_buf = G_realloc(cat_buf, strlen(cat_buf) + size + 1); + p->type = P_CHAR; + p->p_char = cat_buf; + check_mem_alloc(p); + cat_buf_len += sprintf(cat_buf + cat_buf_len, "%d)", cat); + } + else + { + size = snprintf(NULL, 0, "%d, ", cat); + cat_buf = G_realloc(cat_buf, strlen(cat_buf) + size + 1); + p->type = P_CHAR; + p->p_char = cat_buf; + check_mem_alloc(p); + cat_buf_len += sprintf(cat_buf + cat_buf_len, "%d, ", cat); + } + } + else if (line == box_list->n_values - 1) + { + size = snprintf(NULL, 0, "%d)", cat); + cat_buf = G_realloc(cat_buf, strlen(cat_buf) + size + 1); + p->type = P_CHAR; + p->p_char = cat_buf; + check_mem_alloc(p); + cat_buf_len += sprintf(cat_buf + cat_buf_len, "%d)", cat); + } + } + + size = snprintf(NULL, 0, " WHERE %s in %s", keycol, cat_buf); + where_buf = G_malloc(size + 1); + p->type = P_CHAR; + p->p_char = where_buf; + check_mem_alloc(p); + + sprintf(where_buf, " WHERE %s in %s", keycol, cat_buf); + if (db_append_string(&sql, where_buf) != DB_OK) + G_fatal_error(_("Unable to append string")); + G_free(cat_buf); + G_free(where_buf); + + ret.sql = G_malloc(sizeof(sql)); + p->type = P_DBSTRING; + check_mem_alloc(p); + p->p_dbString = ret.sql; + ret.sql = &sql; + + ret.ncats = ncats; + + G_debug(1, "Sql statement: %s", ret.sql->string); + G_debug(1, "Number of cats: %d", ret.ncats); + + return ret; } diff --git a/raster/r.carve/lobf.c b/raster/r.carve/lobf.c index 1d6cf96bafa..1110cd790ed 100644 --- a/raster/r.carve/lobf.c +++ b/raster/r.carve/lobf.c @@ -1,10 +1,11 @@ - /**************************************************************************** * * MODULE: r.carve * * AUTHOR(S): Original author Bill Brown, UIUC GIS Laboratory * Brad Douglas + * Tomas Zigo (adding the option + * to read width, depth values from vector map table columns) * * PURPOSE: Takes vector stream data, converts it to 3D raster and * subtracts a specified depth @@ -17,7 +18,7 @@ * ****************************************************************************/ -/* +/* * Routines to create a line of best fit when given a set of coord pairs. */ @@ -50,38 +51,38 @@ double pg_y_from_x(PointGrp * pg, const double x) void pg_addpt(PointGrp * pg, Point2 pt) { if (pg->npts < MAX_PTS - 1) { - double x = pt[0]; - double y = pt[1]; - - /* add point to group */ - pg->pnts[pg->npts][0] = x; - pg->pnts[pg->npts][1] = y; - pg->sum_x += x; - pg->sum_y += y; - pg->sum_xy += (x * y); - pg->sum_x_sq += SQR(x); - ++pg->npts; + double x = pt[0]; + double y = pt[1]; + + /* add point to group */ + pg->pnts[pg->npts][0] = x; + pg->pnts[pg->npts][1] = y; + pg->sum_x += x; + pg->sum_y += y; + pg->sum_xy += (x * y); + pg->sum_x_sq += SQR(x); + ++pg->npts; } if (pg->npts > 1) { - double denom; - - /* solve for x and y using Cramer's Rule */ - - /* check for divide by zero */ - if (0 == - (denom = DET2_2(pg->sum_x_sq, pg->sum_x, pg->sum_x, pg->npts))) { - G_warning(_("trying to divide by zero...no unique solution for " - "system...skipping...")); - pg->slope = pg->yinter = 0.0; - } - else { - pg->slope = DET2_2(pg->sum_xy, pg->sum_x, pg->sum_y, pg->npts) / - denom; - pg->yinter = - DET2_2(pg->sum_x_sq, pg->sum_xy, pg->sum_x, - pg->sum_y) / denom; - } + double denom; + + /* solve for x and y using Cramer's Rule */ + + /* check for divide by zero */ + if (0 == + (denom = DET2_2(pg->sum_x_sq, pg->sum_x, pg->sum_x, pg->npts))) { + G_warning(_("trying to divide by zero...no unique solution for " + "system...skipping...")); + pg->slope = pg->yinter = 0.0; + } + else { + pg->slope = DET2_2(pg->sum_xy, pg->sum_x, pg->sum_y, pg->npts) / + denom; + pg->yinter = + DET2_2(pg->sum_x_sq, pg->sum_xy, pg->sum_x, + pg->sum_y) / denom; + } } } @@ -105,13 +106,13 @@ Point2 *pg_getpoints_reversed(PointGrp * pg) double x, y; for (i = 0; i < iter; i++) { - /* swap points */ - x = pg->pnts[i][0]; - y = pg->pnts[i][1]; - pg->pnts[i][0] = pg->pnts[pg->npts - i - 1][0]; - pg->pnts[i][1] = pg->pnts[pg->npts - i - 1][1]; - pg->pnts[pg->npts - i - 1][0] = x; - pg->pnts[pg->npts - i - 1][1] = y; + /* swap points */ + x = pg->pnts[i][0]; + y = pg->pnts[i][1]; + pg->pnts[i][0] = pg->pnts[pg->npts - i - 1][0]; + pg->pnts[i][1] = pg->pnts[pg->npts - i - 1][1]; + pg->pnts[pg->npts - i - 1][0] = x; + pg->pnts[pg->npts - i - 1][1] = y; } return pg->pnts; diff --git a/raster/r.carve/main.c b/raster/r.carve/main.c index 0e8871ca557..633398d0f33 100644 --- a/raster/r.carve/main.c +++ b/raster/r.carve/main.c @@ -1,10 +1,11 @@ - /**************************************************************************** * * MODULE: r.carve * * AUTHOR(S): Original author Bill Brown, UIUC GIS Laboratory * Brad Douglas + * Tomas Zigo (adding the option + * to read width, depth values from vector map table columns) * * PURPOSE: Takes vector stream data, converts it to 3D raster and * subtracts a specified depth @@ -28,8 +29,8 @@ /* - Note: Use rast input type for rast output - Read vect file, + Note: Use rast input type for rast output + Read vect file, for each line, use a shadow line struct to represent stream profile, where x is distance along stream and y is elevation, @@ -39,13 +40,13 @@ when next pnt elev increases, find next point <= than last confirmed pt just use linear interp for now - write line to new raster + write line to new raster Use orig line struct for XYs, shadow struct Y for cell val if new raster already has value there, use lower? - actually probably want to use val for trunk stream + actually probably want to use val for trunk stream and then verify branch in reverse - for now maybe create a conflict map - After that's working, add width to lines? + After that's working, add width to lines? or use r.grow */ @@ -62,10 +63,17 @@ int main(int argc, char **argv) struct Flag *noflat; const char *vmapset, *rmapset; - int infd, outfd; + unsigned int i; + int clen, infd, outfd; + int width_col_pos = 0, depth_col_pos = 1; struct Map_info Map; struct Map_info outMap; struct Cell_head win; + struct field_info *Fi = NULL; + dbDriver *driver = NULL; + dbString dbstr; + dbColumn *column = NULL; + char *columns[2] = {0}; /* start GIS engine */ G_gisinit(argv[0]); @@ -75,7 +83,7 @@ int main(int argc, char **argv) G_add_keyword(_("hydrology")); module->label = _("Generates stream channels."); module->description = _("Takes vector stream data, transforms it " - "to raster and subtracts depth from the output DEM."); + "to raster and subtracts depth from the output DEM."); parm.inrast = G_define_standard_option(G_OPT_R_INPUT); parm.inrast->key = "raster"; @@ -84,7 +92,7 @@ int main(int argc, char **argv) parm.invect = G_define_standard_option(G_OPT_V_INPUT); parm.invect->key = "vector"; parm.invect->label = - _("Name of input vector map containing stream(s)"); + _("Name of input vector map containing stream(s)"); parm.outrast = G_define_standard_option(G_OPT_R_OUTPUT); @@ -92,7 +100,24 @@ int main(int argc, char **argv) parm.outvect->key = "points"; parm.outvect->required = NO; parm.outvect->description = - _("Name for output vector map for adjusted stream points"); + _("Name for output vector map for adjusted stream points"); + + parm.field = G_define_standard_option(G_OPT_V_FIELD); + parm.field->key = "field"; + parm.field->label = _("Layer number"); + parm.field->guisection = _("Optional"); + + parm.width_col = G_define_standard_option(G_OPT_DB_COLUMN); + parm.width_col->key = "width_column"; + parm.width_col->description = + _("Name of column for 'width' parameter (data type must be numeric)"); + parm.width_col->guisection = _("Optional"); + + parm.depth_col = G_define_standard_option(G_OPT_DB_COLUMN); + parm.depth_col->key = "depth_column"; + parm.depth_col->description = + _("Name of column for 'depth' parameter (data type must be numeric)"); + parm.depth_col->guisection = _("Optional"); width = G_define_option(); width->key = "width"; @@ -109,60 +134,116 @@ int main(int argc, char **argv) noflat->key = 'n'; noflat->description = _("No flat areas allowed in flow direction"); + /* GUI dependency */ + parm.invect->guidependency = G_store(parm.field->key); + /* parse options */ if (G_parser(argc, argv)) - exit(EXIT_FAILURE); + exit(EXIT_FAILURE); G_check_input_output_name(parm.inrast->answer, parm.outrast->answer, - G_FATAL_EXIT); + G_FATAL_EXIT); if (parm.outvect->answer) - Vect_check_input_output_name(parm.invect->answer, parm.outvect->answer, - G_FATAL_EXIT); + Vect_check_input_output_name(parm.invect->answer, parm.outvect->answer, + G_FATAL_EXIT); /* setup lat/lon projection and distance calculations */ init_projection(&win, &parm.wrap); /* default width - one cell at center */ - if (width->answer == NULL) { - parm.swidth = G_distance((win.east + win.west) / 2, - (win.north + win.south) / 2, - ((win.east + win.west) / 2) + win.ew_res, - (win.north + win.south) / 2); + if (!(width->answer)) + { + parm.swidth = 0.0; + adjust_swidth(&win, &parm.swidth); } - else { - if (sscanf(width->answer, "%lf", &parm.swidth) != 1) { - G_warning(_("Invalid width value '%s' - using default."), - width->answer); - parm.swidth = - G_distance((win.east + win.west) / 2, - (win.north + win.south) / 2, - ((win.east + win.west) / 2) + win.ew_res, - (win.north + win.south) / 2); - } + else + { + if (sscanf(width->answer, "%lf", &parm.swidth) != 1 || parm.swidth < 0.0) + { + G_warning(_("Invalid width value '%s' - using default."), + width->answer); + adjust_swidth(&win, &parm.swidth); + } } - if (depth->answer == NULL) - parm.sdepth = 0.0; - else { - if (sscanf(depth->answer, "%lf", &parm.sdepth) != 1) { - G_warning(_("Invalid depth value '%s' - using default."), - depth->answer); - parm.sdepth = 0.0; - } + if (!(depth->answer)) + { + adjust_sdepth(&parm.sdepth); } + else + { + if (sscanf(depth->answer, "%lf", &parm.sdepth) != 1 || parm.sdepth < 0.0) + { + G_warning(_("Invalid depth value '%s' - using default."), + depth->answer); + adjust_sdepth(&parm.sdepth); + } + } parm.noflat = noflat->answer; /* open input files */ if ((vmapset = G_find_vector2(parm.invect->answer, "")) == NULL) - G_fatal_error(_("Vector map <%s> not found"), parm.invect->answer); + G_fatal_error(_("Vector map <%s> not found"), parm.invect->answer); Vect_set_open_level(2); if (Vect_open_old(&Map, parm.invect->answer, vmapset) < 0) - G_fatal_error(_("Unable to open vector map <%s>"), parm.invect->answer); + G_fatal_error(_("Unable to open vector map <%s>"), parm.invect->answer); if ((rmapset = G_find_file2("cell", parm.inrast->answer, "")) == NULL) - G_fatal_error(_("Raster map <%s> not found"), parm.inrast->answer); + G_fatal_error(_("Raster map <%s> not found"), parm.inrast->answer); + + // Open database driver + db_init_string(&dbstr); + + if (parm.field->answer) + { + Fi = Vect_get_field2(&Map, parm.field->answer); + if (!(Fi)) + G_fatal_error(_("Database connection not defined for layer <%s>"), + parm.field->answer); + + driver = db_start_driver_open_database(Fi->driver, Fi->database); + if (!(driver)) + G_fatal_error(_("Unable to open database <%s> by driver <%s>"), + Fi->database, Fi->driver); + + columns[width_col_pos] = parm.width_col->answer; + columns[depth_col_pos] = parm.depth_col->answer; + + clen = sizeof(columns) / sizeof(columns[0]); + for (i = 0; i < clen; i++) + { + if (columns[i]) + { + // Check if to_column exists and get its SQL type + db_get_column(driver, Fi->table, columns[i], &column); + + if (column) + { + db_free_column(column); + column = NULL; + int ctype = db_column_Ctype(driver, Fi->table, columns[i]); + if (ctype != DB_C_TYPE_INT && ctype != DB_C_TYPE_DOUBLE) + { + // Close db connection + db_close_database_shutdown_driver(driver); + driver = NULL; + G_fatal_error(_("Incompatible column type for <%s> column"), + columns[i]); + } + } + else + { + // Close db connection + db_close_database_shutdown_driver(driver); + driver = NULL; + G_fatal_error(_("Column <%s> not found in table <%s>"), + columns[i], Fi->table); + } + } + } + } infd = Rast_open_old(parm.inrast->answer, rmapset); @@ -173,16 +254,17 @@ int main(int argc, char **argv) /* if specified, open vector for output */ if (parm.outvect->answer) - open_new_vect(&outMap, parm.outvect->answer); + open_new_vect(&outMap, parm.outvect->answer); - enforce_downstream(infd, outfd, &Map, &outMap, &parm); + enforce_downstream(infd, outfd, &Map, &outMap, &parm, Fi, + &width_col_pos, &depth_col_pos, columns, driver); Rast_close(infd); Rast_close(outfd); close_vect(&Map, 0); if (parm.outvect->answer) - close_vect(&outMap, 1); + close_vect(&outMap, 1); /* write command line to history file */ update_rast_history(&parm); @@ -200,22 +282,22 @@ static int init_projection(struct Cell_head *window, int *wrap_ncols) G_get_set_window(window); if (((window->west == (window->east - 360.0)) || - (window->east == (window->west - 360.0))) && - (G_projection() == PROJECTION_LL)) { + (window->east == (window->west - 360.0))) && + (G_projection() == PROJECTION_LL)) { #if 0 - G_get_ellipsoid_parameters(&a, &e2); - G_begin_geodesic_distance(a, e2); + G_get_ellipsoid_parameters(&a, &e2); + G_begin_geodesic_distance(a, e2); - /* add 1.1, not 1.0, to ensure that we round up */ - *wrap_ncols = - (360.0 - (window->east - window->west)) / window->ew_res + 1.1; + /* add 1.1, not 1.0, to ensure that we round up */ + *wrap_ncols = + (360.0 - (window->east - window->west)) / window->ew_res + 1.1; #else - G_fatal_error(_("Lat/Long location is not supported by %s. Please reproject map first."), - G_program_name()); + G_fatal_error(_("Lat/Long location is not supported by %s. Please reproject map first."), + G_program_name()); #endif } else { - *wrap_ncols = 0; + *wrap_ncols = 0; } G_begin_distance_calculations(); diff --git a/raster/r.carve/r.carve.html b/raster/r.carve/r.carve.html index a43ab747bbe..b597052f0f6 100644 --- a/raster/r.carve/r.carve.html +++ b/raster/r.carve/r.carve.html @@ -3,7 +3,15 @@

DESCRIPTION

r.carve accepts vector stream data as input, transforms them to raster, and subtracts a default-depth + additional-depth from a DEM. If the given width is more than 1 cell, it will carve the stream with the -given width. With the -n flag it should eliminate all flat cells within +given width. It is possible to define different width and depth (column +in the vector map table) separately for each line segment. Another way is to +combine the value from the vector map table column (represent width or depth) +and the width or depth parameter. If the width or depth value in the +vector map table is NULL, the value of the width or depth argument +parameter is used if defined, otherwise the default value for width is +used, which represents half the resolution of the region, and for depth +it is 0.0. If the width and depth values are less than 0.0 m, the default +value is used. With the -n flag it should eliminate all flat cells within the stream, so when and if the water gets into the stream it will flow. The points option generates x,y,z for points which define the stream with the z-value of the bottom of the carved-in stream. These points @@ -97,6 +105,68 @@

EXAMPLE

+
+Another width_columns, depth_columns, width, depth parameters combinations
+
+# clipped streams vector map (clipped by region)
+# added width, depth columns with following values
+# db.select table=streams_clipped sql="SELECT cat, width, depth FROM streams_clipped"
+cat|width|depth
+95364||
+95371||
+95379||
+95395|-99|
+95424|50|8
+95425|20|2
+
+# width and depth values are read from vector map table columns
+# if value in column is less than 0.0 m or NULL, default value is used (for width parameter default value is half of region resolution and for depth is 0.0 m)
+r.carve rast=elev_lid792_1m vect=streams_clipped out=carved_dem width_column=width depth_column=depth
+
+# debug info
+D0/0: Process line with id: 4, cat: 95364, width: 0.50, depth: 0.00
+D0/0: Process line with id: 3, cat: 95371, width: 0.50, depth: 0.00
+D0/0: Process line with id: 5, cat: 95379, width: 0.50, depth: 0.00
+D0/0: Process line with id: 6, cat: 95379, width: 0.50, depth: 0.00
+D0/0: Process line with id: 7, cat: 95395, width: 0.50, depth: 0.00
+D0/0: Process line with id: 1, cat: 95424, width: 25.00, depth: 8.00
+D0/0: Process line with id: 2, cat: 95425, width: 10.00, depth: 2.00
+
+width value = width value / 2.0
+
+# width and depth values are read from vector map table columns
+# if value in column is less than 0.0 m (for width parameter default value is half of region resolution and for depth is 0.0 m)
+# if value in width column is NULL, width, depth parameter argument value is used
+r.carve rast=elev_lid792_1m vect=streams out=carved_dem width_column=width depth_column=depth width=15.0 depth=6.0
+
+# debug info
+D0/0: Process line with id: 4, cat: 95364, width: 7.50, depth: 6.00
+D0/0: Process line with id: 3, cat: 95371, width: 7.50, depth: 6.00
+D0/0: Process line with id: 5, cat: 95379, width: 7.50, depth: 6.00
+D0/0: Process line with id: 6, cat: 95379, width: 7.50, depth: 6.00
+D0/0: Process line with id: 7, cat: 95395, width: 0.50, depth: 6.00
+D0/0: Process line with id: 1, cat: 95424, width: 25.00, depth: 8.00
+D0/0: Process line with id: 2, cat: 95425, width: 10.00, depth: 2.00
+
+width value = width value / 2.0
+
+# width values are read from vector map table columns, and depth is constant value 6.0 m (depth parameter argument)
+# if width value in column is less than 0.0 m, default value is used which is half of region resolution
+# if value in width column is NULL, width parameter argument value is used
+r.carve rast=elev_lid792_1m vect=streams out=carved_dem width_columns=width depth_column=depth width=15.0 depth=6.0
+
+# debug info
+D0/0: Process line with id: 4, cat: 95364, width: 7.50, depth: 6.00
+D0/0: Process line with id: 3, cat: 95371, width: 7.50, depth: 6.00
+D0/0: Process line with id: 5, cat: 95379, width: 7.50, depth: 6.00
+D0/0: Process line with id: 6, cat: 95379, width: 7.50, depth: 6.00
+D0/0: Process line with id: 7, cat: 95395, width: 0.50, depth: 6.00
+D0/0: Process line with id: 1, cat: 95424, width: 25.00, depth: 6.00
+D0/0: Process line with id: 2, cat: 95425, width: 10.00, depth: 6.00
+
+width value = width value / 2.0
+
+

KNOWN ISSUES

@@ -117,12 +187,13 @@

SEE ALSO

r.flow, r.fill.dir, r.watershed - +

AUTHOR

Bill Brown (GMSL)
-GRASS 6 update: Brad Douglas - +GRASS 6 update: Brad Douglas
+Adding the option to read width, depth values from vector map table columns +Tomas Zigo