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..5a4f8ea56fa 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 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 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..76cf2951992 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..e1ec17b0da8 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 @@ -98,6 +107,83 @@

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

The module does not operate yet in latitude-longitude locations. It @@ -117,12 +203,12 @@

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