|
|
@@ -1,752 +0,0 @@
|
|
|
-
|
|
|
-/****************************************************************
|
|
|
- *
|
|
|
- * MODULE: v.buffer
|
|
|
- *
|
|
|
- * AUTHOR(S): Radim Blazek
|
|
|
- *
|
|
|
- * PURPOSE: Vector buffer
|
|
|
- *
|
|
|
- * COPYRIGHT: (C) 2001 by the GRASS Development Team
|
|
|
- *
|
|
|
- * This program is free software under the
|
|
|
- * GNU General Public License (>=v2).
|
|
|
- * Read the file COPYING that comes with GRASS
|
|
|
- * for details.
|
|
|
- *
|
|
|
- **************************************************************/
|
|
|
-#include <stdlib.h>
|
|
|
-#include <string.h>
|
|
|
-#include <math.h>
|
|
|
-#include <grass/gis.h>
|
|
|
-#include <grass/Vect.h>
|
|
|
-#include <grass/dbmi.h>
|
|
|
-#include <grass/glocale.h>
|
|
|
-
|
|
|
-#define DEBUG_NONE 0
|
|
|
-#define DEBUG_BUFFER 1
|
|
|
-#define DEBUG_CLEAN 2
|
|
|
-
|
|
|
-#define PI M_PI
|
|
|
-
|
|
|
-/* TODO: look at RET value and use, is it OK? */
|
|
|
-#define RET 0.000000001 /* Representation error tolerance */
|
|
|
-
|
|
|
-/* returns shortest distance to nearest input feature within the buffer
|
|
|
- * or a number greater than buffer (2*buffer) if not found */
|
|
|
-double input_distance(struct Map_info *In, int type, double buffer, double x,
|
|
|
- double y)
|
|
|
-{
|
|
|
- int i;
|
|
|
- static struct ilist *List = NULL;
|
|
|
- static struct line_pnts *Points = NULL;
|
|
|
- double current_distance;
|
|
|
-
|
|
|
- BOUND_BOX box;
|
|
|
-
|
|
|
- if (List == NULL)
|
|
|
- List = Vect_new_list();
|
|
|
- else
|
|
|
- Vect_reset_list(List);
|
|
|
-
|
|
|
- if (Points == NULL)
|
|
|
- Points = Vect_new_line_struct();
|
|
|
-
|
|
|
-
|
|
|
- /* Inside area ? */
|
|
|
- if (type & GV_AREA) {
|
|
|
- int area, centroid;
|
|
|
-
|
|
|
- /* inside */
|
|
|
- area = Vect_find_area(In, x, y);
|
|
|
- centroid = 0;
|
|
|
- if (area)
|
|
|
- centroid = Vect_get_area_centroid(In, area);
|
|
|
-
|
|
|
- G_debug(3, " area = %d centroid = %d", area, centroid);
|
|
|
- if (centroid)
|
|
|
- return 0.0;
|
|
|
- }
|
|
|
-
|
|
|
- /* ouside area, within buffer? */
|
|
|
- /* The centroid is in buffer if at least one point/line/boundary is in buffer distance */
|
|
|
- box.E = x + buffer;
|
|
|
- box.W = x - buffer;
|
|
|
- box.N = y + buffer;
|
|
|
- box.S = y - buffer;
|
|
|
- box.T = PORT_DOUBLE_MAX;
|
|
|
- box.B = -PORT_DOUBLE_MAX;
|
|
|
-
|
|
|
- Vect_select_lines_by_box(In, &box, GV_POINTS | GV_LINES, List);
|
|
|
- G_debug(3, " %d lines selected by box", List->n_values);
|
|
|
-
|
|
|
- current_distance = 2 * buffer;
|
|
|
- for (i = 0; i < List->n_values; i++) {
|
|
|
- int line, ltype;
|
|
|
- double dist;
|
|
|
-
|
|
|
- line = List->value[i];
|
|
|
-
|
|
|
- G_debug(3, " line[%d] = %d", i, line);
|
|
|
-
|
|
|
- ltype = Vect_read_line(In, Points, NULL, line);
|
|
|
-
|
|
|
- Vect_line_distance(Points, x, y, 0., 0, NULL, NULL, NULL, &dist, NULL,
|
|
|
- NULL);
|
|
|
- G_debug(3, " dist = %f", dist);
|
|
|
- if (dist > buffer)
|
|
|
- continue;
|
|
|
-
|
|
|
- /* lines */
|
|
|
- if (type & ltype) {
|
|
|
- if (dist < current_distance)
|
|
|
- current_distance = dist;
|
|
|
- }
|
|
|
-
|
|
|
- /* Areas */
|
|
|
- if ((type & GV_AREA) && ltype == GV_BOUNDARY) {
|
|
|
- int j, side[2], centr[2], area_in;
|
|
|
-
|
|
|
- Vect_get_line_areas(In, line, &side[0], &side[1]);
|
|
|
-
|
|
|
- for (j = 0; j < 2; j++) {
|
|
|
- centr[j] = 0;
|
|
|
-
|
|
|
- if (side[j] > 0)
|
|
|
- area_in = side[j];
|
|
|
- else /* island */
|
|
|
- area_in = Vect_get_isle_area(In, abs(side[j]));
|
|
|
-
|
|
|
- if (area_in > 0)
|
|
|
- centr[j] = Vect_get_area_centroid(In, area_in);
|
|
|
- }
|
|
|
-
|
|
|
- if (centr[0] || centr[1]) {
|
|
|
- if (dist < current_distance)
|
|
|
- current_distance = dist;
|
|
|
- }
|
|
|
- }
|
|
|
- }
|
|
|
- return current_distance;
|
|
|
-}
|
|
|
-
|
|
|
-/*
|
|
|
- * Test if area in Out is in buffer.
|
|
|
- * Return: 1 in buffer
|
|
|
- * 0 outside buffer
|
|
|
- */
|
|
|
-int area_in_buffer(struct Map_info *In, struct Map_info *Out,
|
|
|
- int area, int type, double buffer, double tolerance)
|
|
|
-{
|
|
|
- double dist;
|
|
|
- int ret, i;
|
|
|
- static struct ilist *List = NULL;
|
|
|
- static struct line_pnts *Points = NULL;
|
|
|
- double x, y;
|
|
|
-
|
|
|
- G_debug(3, " area_in_buffer(): area = %d", area);
|
|
|
-
|
|
|
- if (List == NULL)
|
|
|
- List = Vect_new_list();
|
|
|
-
|
|
|
- if (Points == NULL)
|
|
|
- Points = Vect_new_line_struct();
|
|
|
-
|
|
|
- /* Warning: because of tolerance, centroid in area calculated by Vect_get_point_in_area()
|
|
|
- * may fall into buffer (arc interpolated by polygon) even if area is outside!!! */
|
|
|
-
|
|
|
- /* Because of finite number of decimal places in double representation, RET is used. */
|
|
|
-
|
|
|
- /* Test1: Centroid (the only reliable, I think).
|
|
|
- * There are 3 possible cases for the distance of centroid to nearest input feature:
|
|
|
- * 1) < (buffer - tolerance) => area inside the buffer
|
|
|
- * 2) > buffer => area outside the buffer
|
|
|
- * 3) > (buffer - tolerance) and < buffer (that means within the tolerance) => uncertain */
|
|
|
- ret = Vect_get_point_in_area(Out, area, &x, &y);
|
|
|
- if (ret == 0) { /* centroid OK */
|
|
|
- /* test the distance to nearest input feature */
|
|
|
- dist = input_distance(In, type, buffer, x, y);
|
|
|
- G_debug(3, " centroid %f, %f dist = %f", x, y, dist);
|
|
|
-
|
|
|
- if (dist < (buffer - tolerance - RET)) {
|
|
|
- G_debug(3, " centroid in buffer -> area in buffer");
|
|
|
- return 1;
|
|
|
- }
|
|
|
- else if (dist > (buffer + RET)) {
|
|
|
- G_debug(3, " centroid outside buffer -> area outside buffer");
|
|
|
- return 0;
|
|
|
- }
|
|
|
- }
|
|
|
-
|
|
|
- /* Test2: counterclockwise (ccw) boundary
|
|
|
- * Boundaries around features are written in ccw order, that means area inside buffer is on the
|
|
|
- * left side. Area bundaries are listed in cw order (ccw boundaries as negative).
|
|
|
- * So if any boundary in area list is negative, i.e. in cww order, i.e. buffer inside area,
|
|
|
- * whole area _must_ be in buffer. But, also areas without such boundary may be in buffer,
|
|
|
- * see below. */
|
|
|
- /* TODO: The problem!!!: unfortunately it is not true. After snap, break and remove duplicate
|
|
|
- * it may happen that ccw line appeareas in the area outside the buffer!!! */
|
|
|
- Vect_get_area_boundaries(Out, area, List);
|
|
|
-
|
|
|
- for (i = 0; i < List->n_values; i++) {
|
|
|
- if (List->value[i] < 0) {
|
|
|
- G_debug(3, " negative boundary -> area in buffer");
|
|
|
- return 1;
|
|
|
- }
|
|
|
- }
|
|
|
-
|
|
|
- /* Test3: Input feature within area.
|
|
|
- * I believe, that if no negative boundary was found and area is in buffer,
|
|
|
- * the distance of the nearest input feature for at least one area vertex must be < (buffer-tolerance)
|
|
|
- * This test is left as last because it is slow (check each vertex). */
|
|
|
- Vect_get_area_points(Out, area, Points);
|
|
|
- for (i = 0; i < Points->n_points - 1; i++) {
|
|
|
- dist = input_distance(In, type, buffer, Points->x[i], Points->y[i]);
|
|
|
- G_debug(3, " vertex dist = %f", dist);
|
|
|
-
|
|
|
- if (dist < (buffer - tolerance - RET)) {
|
|
|
- G_debug(3, " vertex in buffer -> area in buffer");
|
|
|
- return 1;
|
|
|
- }
|
|
|
- }
|
|
|
-
|
|
|
- /* Test4?: It may be that that Test3 does not cover all possible cases. If so, somebody must
|
|
|
- * find such example and send it to me */
|
|
|
-
|
|
|
- G_debug(3, " -> area outside buffer");
|
|
|
- return 0;
|
|
|
-}
|
|
|
-
|
|
|
-void stop(struct Map_info *In, struct Map_info *Out)
|
|
|
-{
|
|
|
- Vect_close(In);
|
|
|
-
|
|
|
- G_message(_("Rebuilding topology..."));
|
|
|
- Vect_build_partial(Out, GV_BUILD_NONE);
|
|
|
- Vect_build(Out);
|
|
|
- Vect_close(Out);
|
|
|
-}
|
|
|
-
|
|
|
-int main(int argc, char *argv[])
|
|
|
-{
|
|
|
- struct Map_info In, Out;
|
|
|
- struct line_pnts *Points, *BPoints;
|
|
|
- struct line_cats *Cats, *BCats;
|
|
|
- struct GModule *module;
|
|
|
- struct Option *in_opt, *out_opt, *type_opt, *buffer_opt, *tolerance_opt,
|
|
|
- *bufcol_opt, *scale_opt, *debug_opt, *field_opt;
|
|
|
- double buffer, tolerance, dtmp;
|
|
|
- int type, debug;
|
|
|
- int ret, nareas, area, nlines, line;
|
|
|
- char *Areas, *Lines;
|
|
|
- int field;
|
|
|
-
|
|
|
- /* Attributes if sizecol is used */
|
|
|
- int i, nrec, ctype;
|
|
|
- struct field_info *Fi;
|
|
|
- dbDriver *Driver;
|
|
|
- dbCatValArray cvarr;
|
|
|
- int size_val_int;
|
|
|
- double size_val, scale, orig_tolerance;
|
|
|
-
|
|
|
-
|
|
|
- module = G_define_module();
|
|
|
- module->keywords = _("vector");
|
|
|
- module->description =
|
|
|
- _("Creates a buffer around features of given type (areas must contain centroid).");
|
|
|
-
|
|
|
- in_opt = G_define_standard_option(G_OPT_V_INPUT);
|
|
|
- out_opt = G_define_standard_option(G_OPT_V_OUTPUT);
|
|
|
-
|
|
|
- type_opt = G_define_standard_option(G_OPT_V_TYPE);
|
|
|
- type_opt->options = "point,line,boundary,centroid,area";
|
|
|
- type_opt->answer = "point,line,area";
|
|
|
-
|
|
|
- field_opt = G_define_standard_option(G_OPT_V_FIELD);
|
|
|
-
|
|
|
- buffer_opt = G_define_option();
|
|
|
- buffer_opt->key = "buffer";
|
|
|
- buffer_opt->type = TYPE_DOUBLE;
|
|
|
- buffer_opt->required = NO;
|
|
|
- buffer_opt->description = _("Buffer distance in map units");
|
|
|
-
|
|
|
- bufcol_opt = G_define_option();
|
|
|
- bufcol_opt->key = "bufcol";
|
|
|
- bufcol_opt->type = TYPE_STRING;
|
|
|
- bufcol_opt->required = NO;
|
|
|
- bufcol_opt->description =
|
|
|
- _("Attribute column to use for buffer distances");
|
|
|
- bufcol_opt->guisection = _("Advanced");
|
|
|
-
|
|
|
- scale_opt = G_define_option();
|
|
|
- scale_opt->key = "scale";
|
|
|
- scale_opt->type = TYPE_DOUBLE;
|
|
|
- scale_opt->required = NO;
|
|
|
- scale_opt->answer = "1.0";
|
|
|
- scale_opt->description = _("Scaling factor for attribute column values");
|
|
|
- scale_opt->guisection = _("Advanced");
|
|
|
-
|
|
|
- tolerance_opt = G_define_option();
|
|
|
- tolerance_opt->key = "tolerance";
|
|
|
- tolerance_opt->type = TYPE_DOUBLE;
|
|
|
- tolerance_opt->required = NO;
|
|
|
- tolerance_opt->answer = "0.01";
|
|
|
- tolerance_opt->guisection = _("Advanced");
|
|
|
- tolerance_opt->description =
|
|
|
- _("Maximum distance between theoretical arc and polygon segments "
|
|
|
- "as multiple of buffer");
|
|
|
-
|
|
|
- debug_opt = G_define_option();
|
|
|
- debug_opt->key = "debug";
|
|
|
- debug_opt->type = TYPE_STRING;
|
|
|
- debug_opt->required = NO;
|
|
|
- debug_opt->options = "buffer,clean";
|
|
|
- debug_opt->guisection = _("Advanced");
|
|
|
- debug_opt->description = _("Stop the process at a certain stage");
|
|
|
-
|
|
|
- G_gisinit(argv[0]);
|
|
|
- if (G_parser(argc, argv))
|
|
|
- exit(EXIT_FAILURE);
|
|
|
-
|
|
|
- type = Vect_option_to_types(type_opt);
|
|
|
- field = atoi(field_opt->answer);
|
|
|
-
|
|
|
- if ((buffer_opt->answer && bufcol_opt->answer) ||
|
|
|
- (!(buffer_opt->answer || bufcol_opt->answer)))
|
|
|
- G_fatal_error("Select a buffer distance or column, but not both.");
|
|
|
-
|
|
|
- if (bufcol_opt->answer)
|
|
|
- G_warning(_("The bufcol option may contain bugs during the cleaning "
|
|
|
- "step. If you encounter problems, use the debug "
|
|
|
- "option or clean manually with v.clean tool=break; "
|
|
|
- "v.category step=0; v.extract -d type=area"));
|
|
|
-
|
|
|
- orig_tolerance = atof(tolerance_opt->answer);
|
|
|
- tolerance = orig_tolerance;
|
|
|
-
|
|
|
- scale = atof(scale_opt->answer);
|
|
|
- if (scale <= 0.0)
|
|
|
- G_fatal_error("Illegal scale value");
|
|
|
-
|
|
|
- if (buffer_opt->answer) {
|
|
|
- buffer = fabs(atof(buffer_opt->answer));
|
|
|
-
|
|
|
- tolerance *= buffer;
|
|
|
- G_message(_("The tolerance in map units: %g"), tolerance);
|
|
|
-
|
|
|
- /* At least 8 points for circle. */
|
|
|
- dtmp = 0.999 * buffer * (1 - cos(2 * PI / 8 / 2));
|
|
|
- G_debug(3, "Minimum tolerance = %f", dtmp);
|
|
|
- if (tolerance > dtmp) {
|
|
|
- tolerance = dtmp;
|
|
|
- G_warning(_("The tolerance was reset to %g (map units)"),
|
|
|
- tolerance);
|
|
|
- }
|
|
|
- }
|
|
|
-
|
|
|
- debug = DEBUG_NONE;
|
|
|
- if (debug_opt->answer) {
|
|
|
- if (debug_opt->answer[0] == 'b')
|
|
|
- debug = DEBUG_BUFFER;
|
|
|
- else if (debug_opt->answer[0] == 'c')
|
|
|
- debug = DEBUG_CLEAN;
|
|
|
- }
|
|
|
-
|
|
|
- Vect_check_input_output_name(in_opt->answer, out_opt->answer,
|
|
|
- GV_FATAL_EXIT);
|
|
|
-
|
|
|
- Points = Vect_new_line_struct();
|
|
|
- BPoints = Vect_new_line_struct();
|
|
|
- Cats = Vect_new_cats_struct();
|
|
|
- BCats = Vect_new_cats_struct();
|
|
|
-
|
|
|
- /* open input vector */
|
|
|
- Vect_set_open_level(2);
|
|
|
- Vect_open_old(&In, in_opt->answer, "");
|
|
|
-
|
|
|
- Vect_set_fatal_error(GV_FATAL_PRINT);
|
|
|
- if (0 > Vect_open_new(&Out, out_opt->answer, 0)) {
|
|
|
- Vect_close(&In);
|
|
|
- exit(EXIT_FAILURE);
|
|
|
- }
|
|
|
-
|
|
|
- /* check and load attribute column data */
|
|
|
- if (bufcol_opt->answer) {
|
|
|
- db_CatValArray_init(&cvarr);
|
|
|
-
|
|
|
- Fi = Vect_get_field(&In, field);
|
|
|
- if (Fi == NULL)
|
|
|
- G_fatal_error(_("Unable to get layer info for vector map"));
|
|
|
-
|
|
|
- Driver = db_start_driver_open_database(Fi->driver, Fi->database);
|
|
|
- if (Driver == NULL)
|
|
|
- G_fatal_error(_("Unable to open database <%s> by driver <%s>"),
|
|
|
- Fi->database, Fi->driver);
|
|
|
-
|
|
|
- /* Note do not check if the column exists in the table because it may be expression */
|
|
|
-
|
|
|
- /* TODO: only select values we need instead of all in column */
|
|
|
- nrec = db_select_CatValArray(Driver, Fi->table, Fi->key,
|
|
|
- bufcol_opt->answer, NULL, &cvarr);
|
|
|
-
|
|
|
- if (nrec < 0)
|
|
|
- G_fatal_error(_("Unable to select data from table"));
|
|
|
- G_debug(2, "%d records selected from table", nrec);
|
|
|
-
|
|
|
- ctype = cvarr.ctype;
|
|
|
- if (ctype != DB_C_TYPE_INT && ctype != DB_C_TYPE_DOUBLE)
|
|
|
- G_fatal_error(_("Column type not supported"));
|
|
|
-
|
|
|
- db_close_database_shutdown_driver(Driver);
|
|
|
-
|
|
|
- for (i = 0; i < cvarr.n_values; i++) {
|
|
|
- if (ctype == DB_C_TYPE_INT) {
|
|
|
- G_debug(4, "cat = %d val = %d", cvarr.value[i].cat,
|
|
|
- cvarr.value[i].val.i);
|
|
|
- }
|
|
|
- else if (ctype == DB_C_TYPE_DOUBLE) {
|
|
|
- G_debug(4, "cat = %d val = %f", cvarr.value[i].cat,
|
|
|
- cvarr.value[i].val.d);
|
|
|
- }
|
|
|
- }
|
|
|
-
|
|
|
- }
|
|
|
-
|
|
|
-
|
|
|
- Vect_copy_head_data(&In, &Out);
|
|
|
- Vect_hist_copy(&In, &Out);
|
|
|
- Vect_hist_command(&Out);
|
|
|
-
|
|
|
- /* Create buffers' boundaries */
|
|
|
-
|
|
|
- /* Lines (and Points) */
|
|
|
- if ((type & GV_POINTS) || (type & GV_LINES)) {
|
|
|
- int nlines, line, ltype;
|
|
|
-
|
|
|
- nlines = Vect_get_num_lines(&In);
|
|
|
-
|
|
|
- G_message(_("Lines buffers... "));
|
|
|
- for (line = 1; line <= nlines; line++) {
|
|
|
- int cat;
|
|
|
-
|
|
|
- G_debug(3, "line = %d", line);
|
|
|
- G_percent(line, nlines, 2);
|
|
|
-
|
|
|
- ltype = Vect_read_line(&In, Points, Cats, line);
|
|
|
- if (!(ltype & type))
|
|
|
- continue;
|
|
|
-
|
|
|
- if (!Vect_cat_get(Cats, field, &cat))
|
|
|
- continue;
|
|
|
-
|
|
|
- if (bufcol_opt->answer) {
|
|
|
- /* get value from sizecol column */
|
|
|
- /* should probably be put in a function */
|
|
|
- if (ctype == DB_C_TYPE_INT) {
|
|
|
- ret =
|
|
|
- db_CatValArray_get_value_int(&cvarr, cat,
|
|
|
- &size_val_int);
|
|
|
- if (ret != DB_OK) {
|
|
|
- G_warning(_("No record for category %d in table <%s>"),
|
|
|
- cat, Fi->table);
|
|
|
- continue;
|
|
|
- }
|
|
|
- size_val = (double)size_val_int;
|
|
|
- }
|
|
|
-
|
|
|
- if (ctype == DB_C_TYPE_DOUBLE) {
|
|
|
- ret =
|
|
|
- db_CatValArray_get_value_double(&cvarr, cat,
|
|
|
- &size_val);
|
|
|
- if (ret != DB_OK) {
|
|
|
- G_warning(_("No record for category %d in table <%s>"),
|
|
|
- cat, Fi->table);
|
|
|
- continue;
|
|
|
- }
|
|
|
- }
|
|
|
-
|
|
|
- if (size_val < 0.0) {
|
|
|
- G_warning(_("Attribute is of invalid size (%.3f) for category %d"),
|
|
|
- size_val, cat);
|
|
|
- continue;
|
|
|
- }
|
|
|
-
|
|
|
- if (size_val == 0.0)
|
|
|
- continue;
|
|
|
-
|
|
|
- buffer = size_val * scale;
|
|
|
- G_debug(2, " dynamic buffer size = %.2f", buffer);
|
|
|
-
|
|
|
- tolerance = orig_tolerance * buffer;
|
|
|
- G_debug(2, _("The tolerance in map units: %g"), tolerance);
|
|
|
-
|
|
|
- /* At least 8 points for circle. */
|
|
|
- dtmp = 0.999 * buffer * (1 - cos(2 * PI / 8 / 2));
|
|
|
- G_debug(2, "Minimum tolerance = %f", dtmp);
|
|
|
- if (tolerance > dtmp) {
|
|
|
- tolerance = dtmp;
|
|
|
- G_warning(_("The tolerance was reset to %g (map units). [category %d]"),
|
|
|
- tolerance, cat);
|
|
|
- }
|
|
|
- }
|
|
|
-
|
|
|
-
|
|
|
- Vect_line_buffer(Points, buffer, tolerance, BPoints);
|
|
|
- Vect_write_line(&Out, GV_BOUNDARY, BPoints, BCats);
|
|
|
- }
|
|
|
- }
|
|
|
-
|
|
|
- /* Areas */
|
|
|
- if (type & GV_AREA) {
|
|
|
- int i, nareas, area, centroid, nisles, isle;
|
|
|
-
|
|
|
- nareas = Vect_get_num_areas(&In);
|
|
|
-
|
|
|
- G_message(_("Areas buffers... "));
|
|
|
- for (area = 1; area <= nareas; area++) {
|
|
|
- int cat;
|
|
|
-
|
|
|
- G_percent(area, nareas, 2);
|
|
|
-
|
|
|
- centroid = Vect_get_area_centroid(&In, area);
|
|
|
- if (centroid == 0)
|
|
|
- continue;
|
|
|
-
|
|
|
- Vect_read_line(&In, NULL, Cats, centroid);
|
|
|
- if (!Vect_cat_get(Cats, field, &cat))
|
|
|
- continue;
|
|
|
-
|
|
|
- if (bufcol_opt->answer) {
|
|
|
- /* get value from sizecol column */
|
|
|
-
|
|
|
- if (ctype == DB_C_TYPE_INT) {
|
|
|
- ret =
|
|
|
- db_CatValArray_get_value_int(&cvarr, cat,
|
|
|
- &size_val_int);
|
|
|
- if (ret != DB_OK) {
|
|
|
- G_warning(_("No record for category %d in table <%s>"),
|
|
|
- cat, Fi->table);
|
|
|
- continue;
|
|
|
- }
|
|
|
- size_val = (double)size_val_int;
|
|
|
- }
|
|
|
-
|
|
|
- if (ctype == DB_C_TYPE_DOUBLE) {
|
|
|
- ret =
|
|
|
- db_CatValArray_get_value_double(&cvarr, cat,
|
|
|
- &size_val);
|
|
|
- if (ret != DB_OK) {
|
|
|
- G_warning(_("No record for category %d in table <%s>"),
|
|
|
- cat, Fi->table);
|
|
|
- continue;
|
|
|
- }
|
|
|
- }
|
|
|
-
|
|
|
- if (size_val < 0.0) {
|
|
|
- G_warning(_("Attribute is of invalid size (%.3f) for category %d"),
|
|
|
- size_val, cat);
|
|
|
- continue;
|
|
|
- }
|
|
|
-
|
|
|
- if (size_val == 0.0)
|
|
|
- continue;
|
|
|
-
|
|
|
- buffer = size_val * scale;
|
|
|
- G_debug(2, " dynamic buffer size = %.2f", buffer);
|
|
|
-
|
|
|
- tolerance = orig_tolerance * buffer;
|
|
|
- G_debug(2, _("The tolerance in map units: %g"), tolerance);
|
|
|
-
|
|
|
- /* At least 8 points for circle. */
|
|
|
- dtmp = 0.999 * buffer * (1 - cos(2 * PI / 8 / 2));
|
|
|
- G_debug(2, "Minimum tolerance = %f", dtmp);
|
|
|
- if (tolerance > dtmp) {
|
|
|
- tolerance = dtmp;
|
|
|
- G_warning(_("The tolerance was reset to %g (map units). [category %d]"),
|
|
|
- tolerance, cat);
|
|
|
- }
|
|
|
- }
|
|
|
-
|
|
|
-
|
|
|
- /* outer ring */
|
|
|
- Vect_get_area_points(&In, area, Points);
|
|
|
- Vect_line_buffer(Points, buffer, tolerance, BPoints);
|
|
|
- Vect_write_line(&Out, GV_BOUNDARY, BPoints, BCats);
|
|
|
-
|
|
|
- /* islands */
|
|
|
- nisles = Vect_get_area_num_isles(&In, area);
|
|
|
- for (i = 0; i < nisles; i++) {
|
|
|
- double l;
|
|
|
-
|
|
|
- isle = Vect_get_area_isle(&In, area, i);
|
|
|
- Vect_get_isle_points(&In, isle, Points);
|
|
|
-
|
|
|
- /* Check if the isle is big enough */
|
|
|
- l = Vect_line_length(Points);
|
|
|
- if (l / 2 < 2 * buffer)
|
|
|
- continue;
|
|
|
-
|
|
|
- Vect_line_buffer(Points, buffer, tolerance, BPoints);
|
|
|
- Vect_write_line(&Out, GV_BOUNDARY, BPoints, BCats);
|
|
|
- }
|
|
|
- }
|
|
|
- }
|
|
|
-
|
|
|
- if (debug == DEBUG_BUFFER) {
|
|
|
- stop(&In, &Out);
|
|
|
- exit(EXIT_SUCCESS);
|
|
|
- }
|
|
|
-
|
|
|
- /* Create areas */
|
|
|
-
|
|
|
- /* Break lines */
|
|
|
- G_message(_("Building parts of topology..."));
|
|
|
- Vect_build_partial(&Out, GV_BUILD_BASE);
|
|
|
-
|
|
|
- /* Warning: snapping must be done, otherwise colinear boundaries are not broken and
|
|
|
- * topology cannot be built (the same angle). But snapping distance must be very, very
|
|
|
- * small, otherwise counterclockwise boundaries can appear in areas outside the buffer.
|
|
|
- * I have done some tests on real data (projected) and threshold 1e-8 was not enough,
|
|
|
- * Snapping threshold 1e-7 seems to work. Don't increase until we find example
|
|
|
- * where it is not sufficient. */
|
|
|
- /* TODO: look at snapping threshold better, calculate some theoretical value to avoid
|
|
|
- * the same angles of lines at nodes, don't forget about LongLat data, probably
|
|
|
- * calculate different threshold for each map, depending on map's bounding box */
|
|
|
- G_message(_("Snapping boundaries..."));
|
|
|
- Vect_snap_lines(&Out, GV_BOUNDARY, 1e-7, NULL);
|
|
|
-
|
|
|
- G_message(_("Breaking boundaries..."));
|
|
|
- Vect_break_lines(&Out, GV_BOUNDARY, NULL);
|
|
|
-
|
|
|
- G_message(_("Removing duplicates..."));
|
|
|
- Vect_remove_duplicates(&Out, GV_BOUNDARY, NULL);
|
|
|
-
|
|
|
- /* Dangles and bridges don't seem to be necessary if snapping is small enough. */
|
|
|
- /*
|
|
|
- G_message ( "Removing dangles..." );
|
|
|
- Vect_remove_dangles ( &Out, GV_BOUNDARY, -1, NULL, stderr );
|
|
|
-
|
|
|
- G_message ( "Removing bridges..." );
|
|
|
- Vect_remove_bridges ( &Out, NULL, stderr );
|
|
|
- */
|
|
|
-
|
|
|
- G_message(_("Attaching islands..."));
|
|
|
- Vect_build_partial(&Out, GV_BUILD_ATTACH_ISLES);
|
|
|
-
|
|
|
- /* Calculate new centroids for all areas */
|
|
|
- nareas = Vect_get_num_areas(&Out);
|
|
|
- Areas = (char *)G_calloc(nareas + 1, sizeof(char));
|
|
|
-
|
|
|
- for (area = 1; area <= nareas; area++) {
|
|
|
- G_debug(3, "area = %d", area);
|
|
|
-
|
|
|
-/*** BUG *** if dynamic bufcol was used, "buffer" will only hold last value ***/
|
|
|
- ret = area_in_buffer(&In, &Out, area, type, buffer, tolerance);
|
|
|
-
|
|
|
- if (ret) {
|
|
|
- G_debug(3, " -> in buffer");
|
|
|
- Areas[area] = 1;
|
|
|
- }
|
|
|
-
|
|
|
- /* Write out centroid (before check if correct, so that it isd visible for debug) */
|
|
|
- if (debug == DEBUG_CLEAN) {
|
|
|
- double x, y;
|
|
|
-
|
|
|
- ret = Vect_get_point_in_area(&Out, area, &x, &y);
|
|
|
- if (ret < 0) {
|
|
|
- G_warning(_("Cannot calculate area centroid"));
|
|
|
- continue;
|
|
|
- }
|
|
|
- Vect_reset_cats(Cats);
|
|
|
- if (Areas[area])
|
|
|
- Vect_cat_set(Cats, 1, 1);
|
|
|
-
|
|
|
- Vect_reset_line(Points);
|
|
|
- Vect_append_point(Points, x, y, 0.);
|
|
|
- Vect_write_line(&Out, GV_CENTROID, Points, Cats);
|
|
|
- }
|
|
|
- }
|
|
|
-
|
|
|
- if (debug == DEBUG_CLEAN) {
|
|
|
- stop(&In, &Out);
|
|
|
- exit(EXIT_SUCCESS);
|
|
|
- }
|
|
|
-
|
|
|
- /* Make a list of boundaries to be deleted (both sides inside) */
|
|
|
- nlines = Vect_get_num_lines(&Out);
|
|
|
- G_debug(3, "nlines = %d", nlines);
|
|
|
- Lines = (char *)G_calloc(nlines + 1, sizeof(char));
|
|
|
-
|
|
|
- for (line = 1; line <= nlines; line++) {
|
|
|
- int j, side[2], areas[2];
|
|
|
-
|
|
|
- G_debug(3, "line = %d", line);
|
|
|
-
|
|
|
- if (!Vect_line_alive(&Out, line))
|
|
|
- continue;
|
|
|
-
|
|
|
- Vect_get_line_areas(&Out, line, &side[0], &side[1]);
|
|
|
-
|
|
|
- for (j = 0; j < 2; j++) {
|
|
|
- if (side[j] == 0) { /* area/isle not build */
|
|
|
- areas[j] = 0;
|
|
|
- }
|
|
|
- else if (side[j] > 0) { /* area */
|
|
|
- areas[j] = side[j];
|
|
|
- }
|
|
|
- else { /* < 0 -> island */
|
|
|
- areas[j] = Vect_get_isle_area(&Out, abs(side[j]));
|
|
|
- }
|
|
|
- }
|
|
|
-
|
|
|
- G_debug(3, " areas = %d , %d -> Areas = %d, %d", areas[0], areas[1],
|
|
|
- Areas[areas[0]], Areas[areas[1]]);
|
|
|
- if (Areas[areas[0]] && Areas[areas[1]])
|
|
|
- Lines[line] = 1;
|
|
|
- }
|
|
|
- G_free(Areas);
|
|
|
-
|
|
|
- /* Delete boundaries */
|
|
|
- for (line = 1; line <= nlines; line++) {
|
|
|
- if (Lines[line]) {
|
|
|
- G_debug(3, " delete line %d", line);
|
|
|
- Vect_delete_line(&Out, line);
|
|
|
- }
|
|
|
- }
|
|
|
-
|
|
|
- G_free(Lines);
|
|
|
-
|
|
|
- /* Create new centroids */
|
|
|
- Vect_reset_cats(Cats);
|
|
|
- Vect_cat_set(Cats, 1, 1);
|
|
|
- nareas = Vect_get_num_areas(&Out);
|
|
|
- for (area = 1; area <= nareas; area++) {
|
|
|
- double x, y;
|
|
|
-
|
|
|
- G_debug(3, "area = %d", area);
|
|
|
-
|
|
|
- if (!Vect_area_alive(&Out, area))
|
|
|
- continue;
|
|
|
-
|
|
|
-/*** BUG *** if dynamic bufcol was used, "buffer" will only hold last value ***/
|
|
|
- ret = area_in_buffer(&In, &Out, area, type, buffer, tolerance);
|
|
|
-
|
|
|
- if (ret) {
|
|
|
- ret = Vect_get_point_in_area(&Out, area, &x, &y);
|
|
|
- if (ret < 0) {
|
|
|
- G_warning(_("Cannot calculate area centroid"));
|
|
|
- continue;
|
|
|
- }
|
|
|
-
|
|
|
- Vect_reset_line(Points);
|
|
|
- Vect_append_point(Points, x, y, 0.);
|
|
|
- Vect_write_line(&Out, GV_CENTROID, Points, Cats);
|
|
|
- }
|
|
|
- }
|
|
|
-
|
|
|
- G_message(_("Attaching centroids..."));
|
|
|
- Vect_build_partial(&Out, GV_BUILD_CENTROIDS);
|
|
|
-
|
|
|
- stop(&In, &Out);
|
|
|
- exit(EXIT_SUCCESS);
|
|
|
-}
|
Hamish Bowman