diff --git a/raster/r.carve/Makefile b/raster/r.carve/Makefile index 4777cb00890..f452758a600 100644 --- a/raster/r.carve/Makefile +++ b/raster/r.carve/Makefile @@ -2,7 +2,7 @@ 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) diff --git a/raster/r.carve/enforce.h b/raster/r.carve/enforce.h index 3faeb7ce037..43a26e5f177 100644 --- a/raster/r.carve/enforce.h +++ b/raster/r.carve/enforce.h @@ -4,6 +4,8 @@ * * 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 @@ -43,17 +45,55 @@ typedef struct { } PointGrp; 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 *); @@ -68,6 +108,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 7a71dc97def..b148b2f1e15 100644 --- a/raster/r.carve/enforce_ds.c +++ b/raster/r.carve/enforce_ds.c @@ -4,6 +4,8 @@ * * 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 @@ -16,6 +18,7 @@ * ****************************************************************************/ +#include #include #include #include @@ -32,8 +35,8 @@ /* 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); @@ -47,21 +50,53 @@ 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) +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); @@ -77,32 +112,127 @@ int enforce_downstream(int infd, int outfd, struct Map_info *Map, 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 */ @@ -113,8 +243,8 @@ static int process_line(struct Map_info *Map, struct Map_info *outMap, if (!cats) 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()); @@ -123,7 +253,7 @@ static int process_line(struct Map_info *Map, struct Map_info *outMap, 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; @@ -168,7 +298,6 @@ static int process_line(struct Map_info *Map, struct Map_info *outMap, if (do_warn) { G_warning(_("Vect runs out of region and re-enters - " "this case is not yet implemented.")); - retval = 1; } /* now check to see if points go downslope(inorder) or upslope */ @@ -194,10 +323,7 @@ static int process_line(struct Map_info *Map, struct Map_info *outMap, /* 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; } static void clear_bitmap(struct BM *bm) @@ -491,3 +617,240 @@ static void process_line_segment(const int npts, void *rbuf, Point2 *pgxypts, Vect_destroy_line_struct(points); Vect_destroy_cats_struct(cats); } + +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 clause ("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 d9d534fc2be..b94e9863319 100644 --- a/raster/r.carve/lobf.c +++ b/raster/r.carve/lobf.c @@ -4,6 +4,8 @@ * * 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 diff --git a/raster/r.carve/main.c b/raster/r.carve/main.c index a242d3e6f16..45055e83af9 100644 --- a/raster/r.carve/main.c +++ b/raster/r.carve/main.c @@ -4,6 +4,8 @@ * * 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 @@ -58,10 +60,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]); @@ -90,6 +99,23 @@ int main(int argc, char **argv) parm.outvect->description = _("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"; width->type = TYPE_DOUBLE; @@ -105,6 +131,9 @@ 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); @@ -119,30 +148,28 @@ int main(int argc, char **argv) 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) { + if (sscanf(width->answer, "%lf", &parm.swidth) != 1 || + parm.swidth < 0.0) { 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); + adjust_swidth(&win, &parm.swidth); } } - if (depth->answer == NULL) - parm.sdepth = 0.0; + if (!(depth->answer)) { + adjust_sdepth(&parm.sdepth); + } else { - if (sscanf(depth->answer, "%lf", &parm.sdepth) != 1) { + if (sscanf(depth->answer, "%lf", &parm.sdepth) != 1 || + parm.sdepth < 0.0) { G_warning(_("Invalid depth value '%s' - using default."), depth->answer); - parm.sdepth = 0.0; + adjust_sdepth(&parm.sdepth); } } @@ -159,6 +186,53 @@ int main(int argc, char **argv) if ((rmapset = G_find_file2("cell", parm.inrast->answer, "")) == NULL) 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); parm.raster_type = Rast_get_map_type(infd); @@ -170,7 +244,8 @@ int main(int argc, char **argv) if (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); diff --git a/raster/r.carve/r.carve.html b/raster/r.carve/r.carve.html index 3d36aa45256..f81f3323cce 100644 --- a/raster/r.carve/r.carve.html +++ b/raster/r.carve/r.carve.html @@ -3,8 +3,17 @@

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 -the stream, so when and if the water gets into the stream it will +given width. It is possible to define different width and depth (column +in the vector map table) separately for each line segment (width_column, +depth_column parameter). 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(s) 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 m. 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 can then be combined with contours to interpolate a new DEM with @@ -109,6 +118,84 @@

EXAMPLE

+ +Another width_columns, depth_columns, width, depth parameters combinations: + +
+# g.region -p
+projection: 99 (Lambert Conformal Conic)
+zone:       0
+datum:      nad83
+ellipsoid:  a=6378137 es=0.006694380022900787
+north:      220750
+south:      220000
+west:       638300
+east:       639000
+nsres:      1
+ewres:      1
+rows:       750
+cols:       700
+cells:      525000
+
+# streams vector map with added width, depth columns with following values
+# db.select table=streams sql="SELECT cat, width, depth FROM streams WHERE cat in (95364, 95371, 95379, 95395, 95424, 95425)"
+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 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, which is constant value 15.0 m
+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

@@ -134,4 +221,9 @@

SEE ALSO

AUTHORS

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 + diff --git a/raster/r.carve/support.c b/raster/r.carve/support.c index dcb0b20265b..2c99d75c85d 100644 --- a/raster/r.carve/support.c +++ b/raster/r.carve/support.c @@ -4,6 +4,8 @@ * * 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 @@ -43,3 +45,56 @@ void update_rast_history(struct parms *parm) Rast_command_history(&hist); Rast_write_history(parm->outrast->answer, &hist); } + +void check_mem_alloc(struct ptr *pointer) +/* + * Function: check_mem_alloc + * ------------------------- + * Check memory allocation + * + * pointer: ptr struct + * ptr.type: pointer type + * ptr.int|double|char|dbString: pointer + * ptr.p_vect_id_cat_map: line id with matching cat ( + * vect_id_cat_map struct) + * + */ +{ + switch (pointer->type) { + case P_INT: + if (!pointer->p_int) { + G_free(pointer->p_int); + pointer->p_int = NULL; + G_fatal_error(_("Fail to allocate memory")); + } + break; + case P_DOUBLE: + if (!pointer->p_double) { + G_free(pointer->p_double); + pointer->p_double = NULL; + G_fatal_error(_("Fail to allocate memory")); + } + break; + case P_CHAR: + if (!pointer->p_char) { + G_free(pointer->p_char); + pointer->p_char = NULL; + G_fatal_error(_("Fail to allocate memory")); + } + break; + case P_DBSTRING: + if (!pointer->p_dbString) { + db_free_string(pointer->p_dbString); + pointer->p_dbString = NULL; + G_fatal_error(_("Fail to allocate memory")); + } + break; + case P_VECT_ID_CAT_MAP: + if (!pointer->p_vect_id_cat_map) { + G_free(pointer->p_vect_id_cat_map); + pointer->p_vect_id_cat_map = NULL; + G_fatal_error(_("Fail to allocate memory")); + } + break; + } +} diff --git a/raster/r.carve/vect.c b/raster/r.carve/vect.c index 6caba9ae9fb..f714184e022 100644 --- a/raster/r.carve/vect.c +++ b/raster/r.carve/vect.c @@ -4,6 +4,8 @@ * * 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