GRASS Programmer's Manual  6.4.2(2012)
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
remove_areas.c
Go to the documentation of this file.
1 
20 #include <grass/gis.h>
21 #include <grass/Vect.h>
22 #include <grass/glocale.h>
23 
38 int
39 Vect_remove_small_areas(struct Map_info *Map, double thresh,
40  struct Map_info *Err, double *removed_area)
41 {
42  int area, nareas;
43  int nremoved = 0;
44  struct ilist *List;
45  struct ilist *AList;
46  struct line_pnts *Points;
47  struct line_cats *Cats;
48  double size_removed = 0.0;
49 
50  List = Vect_new_list();
51  AList = Vect_new_list();
52  Points = Vect_new_line_struct();
53  Cats = Vect_new_cats_struct();
54 
55  nareas = Vect_get_num_areas(Map);
56  for (area = 1; area <= nareas; area++) {
57  int i, j, centroid, dissolve_neighbour;
58  double length, size;
59 
60  G_percent(area, nareas, 1);
61  G_debug(3, "area = %d", area);
62  if (!Vect_area_alive(Map, area))
63  continue;
64 
65  size = Vect_get_area_area(Map, area);
66  if (size > thresh)
67  continue;
68  size_removed += size;
69 
70  /* The area is smaller than the limit -> remove */
71 
72  /* Remove centroid */
73  centroid = Vect_get_area_centroid(Map, area);
74  if (centroid > 0) {
75  if (Err) {
76  Vect_read_line(Map, Points, Cats, centroid);
77  Vect_write_line(Err, GV_CENTROID, Points, Cats);
78  }
79  Vect_delete_line(Map, centroid);
80  }
81 
82  /* Find the adjacent area with which the longest boundary is shared */
83 
84  Vect_get_area_boundaries(Map, area, List);
85 
86  /* Create a list of neighbour areas */
87  Vect_reset_list(AList);
88  for (i = 0; i < List->n_values; i++) {
89  int line, left, right, neighbour;
90 
91  line = List->value[i];
92 
93  if (!Vect_line_alive(Map, abs(line))) /* Should not happen */
94  G_fatal_error(_("Area is composed of dead boundary"));
95 
96  Vect_get_line_areas(Map, abs(line), &left, &right);
97  if (line > 0)
98  neighbour = left;
99  else
100  neighbour = right;
101 
102  G_debug(4, " line = %d left = %d right = %d neighbour = %d",
103  line, left, right, neighbour);
104 
105  Vect_list_append(AList, neighbour); /* this checks for duplicity */
106  }
107  G_debug(3, "num neighbours = %d", AList->n_values);
108 
109  /* Go through the list of neghours and find that with the longest boundary */
110  dissolve_neighbour = 0;
111  length = -1.0;
112  for (i = 0; i < AList->n_values; i++) {
113  int neighbour1;
114  double l = 0.0;
115 
116  neighbour1 = AList->value[i];
117  G_debug(4, " neighbour1 = %d", neighbour1);
118 
119  for (j = 0; j < List->n_values; j++) {
120  int line, left, right, neighbour2;
121 
122  line = List->value[j];
123  Vect_get_line_areas(Map, abs(line), &left, &right);
124  if (line > 0)
125  neighbour2 = left;
126  else
127  neighbour2 = right;
128 
129  if (neighbour2 == neighbour1) {
130  Vect_read_line(Map, Points, NULL, abs(line));
131  l += Vect_line_length(Points);
132  }
133  }
134  if (l > length) {
135  length = l;
136  dissolve_neighbour = neighbour1;
137  }
138  }
139 
140  G_debug(3, "dissolve_neighbour = %d", dissolve_neighbour);
141 
142  /* Make list of boundaries to be removed */
143  Vect_reset_list(AList);
144  for (i = 0; i < List->n_values; i++) {
145  int line, left, right, neighbour;
146 
147  line = List->value[i];
148  Vect_get_line_areas(Map, abs(line), &left, &right);
149  if (line > 0)
150  neighbour = left;
151  else
152  neighbour = right;
153 
154  G_debug(3, " neighbour = %d", neighbour);
155 
156  if (neighbour == dissolve_neighbour) {
157  Vect_list_append(AList, abs(line));
158  }
159  }
160 
161  /* Remove boundaries */
162  for (i = 0; i < AList->n_values; i++) {
163  int line;
164 
165  line = AList->value[i];
166 
167  if (Err) {
168  Vect_read_line(Map, Points, Cats, line);
169  Vect_write_line(Err, GV_BOUNDARY, Points, Cats);
170  }
171  Vect_delete_line(Map, line);
172  }
173 
174  nremoved++;
175  nareas = Vect_get_num_areas(Map);
176  }
177 
178  if (removed_area)
179  *removed_area = size_removed;
180 
181  return (nremoved);
182 }