GRASS Programmer's Manual  6.4.1(2011)
remove_areas.c
Go to the documentation of this file.
00001 
00020 #include <grass/gis.h>
00021 #include <grass/Vect.h>
00022 #include <grass/glocale.h>
00023 
00038 int
00039 Vect_remove_small_areas(struct Map_info *Map, double thresh,
00040                         struct Map_info *Err, double *removed_area)
00041 {
00042     int area, nareas;
00043     int nremoved = 0;
00044     struct ilist *List;
00045     struct ilist *AList;
00046     struct line_pnts *Points;
00047     struct line_cats *Cats;
00048     double size_removed = 0.0;
00049 
00050     List = Vect_new_list();
00051     AList = Vect_new_list();
00052     Points = Vect_new_line_struct();
00053     Cats = Vect_new_cats_struct();
00054 
00055     nareas = Vect_get_num_areas(Map);
00056     for (area = 1; area <= nareas; area++) {
00057         int i, j, centroid, dissolve_neighbour;
00058         double length, size;
00059 
00060         G_percent(area, nareas, 1);
00061         G_debug(3, "area = %d", area);
00062         if (!Vect_area_alive(Map, area))
00063             continue;
00064 
00065         size = Vect_get_area_area(Map, area);
00066         if (size > thresh)
00067             continue;
00068         size_removed += size;
00069 
00070         /* The area is smaller than the limit -> remove */
00071 
00072         /* Remove centroid */
00073         centroid = Vect_get_area_centroid(Map, area);
00074         if (centroid > 0) {
00075             if (Err) {
00076                 Vect_read_line(Map, Points, Cats, centroid);
00077                 Vect_write_line(Err, GV_CENTROID, Points, Cats);
00078             }
00079             Vect_delete_line(Map, centroid);
00080         }
00081 
00082         /* Find the adjacent area with which the longest boundary is shared */
00083 
00084         Vect_get_area_boundaries(Map, area, List);
00085 
00086         /* Create a list of neighbour areas */
00087         Vect_reset_list(AList);
00088         for (i = 0; i < List->n_values; i++) {
00089             int line, left, right, neighbour;
00090 
00091             line = List->value[i];
00092 
00093             if (!Vect_line_alive(Map, abs(line)))       /* Should not happen */
00094                 G_fatal_error(_("Area is composed of dead boundary"));
00095 
00096             Vect_get_line_areas(Map, abs(line), &left, &right);
00097             if (line > 0)
00098                 neighbour = left;
00099             else
00100                 neighbour = right;
00101 
00102             G_debug(4, "  line = %d left = %d right = %d neighbour = %d",
00103                     line, left, right, neighbour);
00104 
00105             Vect_list_append(AList, neighbour); /* this checks for duplicity */
00106         }
00107         G_debug(3, "num neighbours = %d", AList->n_values);
00108 
00109         /* Go through the list of neghours and find that with the longest boundary */
00110         dissolve_neighbour = 0;
00111         length = -1.0;
00112         for (i = 0; i < AList->n_values; i++) {
00113             int neighbour1;
00114             double l = 0.0;
00115 
00116             neighbour1 = AList->value[i];
00117             G_debug(4, "   neighbour1 = %d", neighbour1);
00118 
00119             for (j = 0; j < List->n_values; j++) {
00120                 int line, left, right, neighbour2;
00121 
00122                 line = List->value[j];
00123                 Vect_get_line_areas(Map, abs(line), &left, &right);
00124                 if (line > 0)
00125                     neighbour2 = left;
00126                 else
00127                     neighbour2 = right;
00128 
00129                 if (neighbour2 == neighbour1) {
00130                     Vect_read_line(Map, Points, NULL, abs(line));
00131                     l += Vect_line_length(Points);
00132                 }
00133             }
00134             if (l > length) {
00135                 length = l;
00136                 dissolve_neighbour = neighbour1;
00137             }
00138         }
00139 
00140         G_debug(3, "dissolve_neighbour = %d", dissolve_neighbour);
00141 
00142         /* Make list of boundaries to be removed */
00143         Vect_reset_list(AList);
00144         for (i = 0; i < List->n_values; i++) {
00145             int line, left, right, neighbour;
00146 
00147             line = List->value[i];
00148             Vect_get_line_areas(Map, abs(line), &left, &right);
00149             if (line > 0)
00150                 neighbour = left;
00151             else
00152                 neighbour = right;
00153 
00154             G_debug(3, "   neighbour = %d", neighbour);
00155 
00156             if (neighbour == dissolve_neighbour) {
00157                 Vect_list_append(AList, abs(line));
00158             }
00159         }
00160 
00161         /* Remove boundaries */
00162         for (i = 0; i < AList->n_values; i++) {
00163             int line;
00164 
00165             line = AList->value[i];
00166 
00167             if (Err) {
00168                 Vect_read_line(Map, Points, Cats, line);
00169                 Vect_write_line(Err, GV_BOUNDARY, Points, Cats);
00170             }
00171             Vect_delete_line(Map, line);
00172         }
00173 
00174         nremoved++;
00175         nareas = Vect_get_num_areas(Map);
00176     }
00177 
00178     if (removed_area)
00179         *removed_area = size_removed;
00180 
00181     return (nremoved);
00182 }
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Defines