GRASS Programmer's Manual
6.4.1(2011)
|
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 }