GRASS Programmer's Manual 6.4.1(2011)
|
00001 00019 #include <string.h> 00020 #include <stdlib.h> 00021 #include <stdio.h> 00022 #include <grass/glocale.h> 00023 #include <grass/gis.h> 00024 #include <grass/Vect.h> 00025 00037 int Vect_build_line_area(struct Map_info *Map, int iline, int side) 00038 { 00039 int j, area, isle, n_lines, line, type, direction; 00040 static int first = 1; 00041 long offset; 00042 struct Plus_head *plus; 00043 P_LINE *BLine; 00044 static struct line_pnts *Points, *APoints; 00045 plus_t *lines; 00046 double area_size; 00047 00048 plus = &(Map->plus); 00049 00050 G_debug(3, "Vect_build_line_area() line = %d, side = %d", iline, side); 00051 00052 if (first) { 00053 Points = Vect_new_line_struct(); 00054 APoints = Vect_new_line_struct(); 00055 first = 0; 00056 } 00057 00058 area = dig_line_get_area(plus, iline, side); 00059 if (area != 0) { 00060 G_debug(3, " area/isle = %d -> skip", area); 00061 return 0; 00062 } 00063 00064 n_lines = dig_build_area_with_line(plus, iline, side, &lines); 00065 G_debug(3, " n_lines = %d", n_lines); 00066 if (n_lines < 1) { 00067 return 0; 00068 } /* area was not built */ 00069 00070 /* Area or island ? */ 00071 Vect_reset_line(APoints); 00072 for (j = 0; j < n_lines; j++) { 00073 line = abs(lines[j]); 00074 BLine = plus->Line[line]; 00075 offset = BLine->offset; 00076 G_debug(3, " line[%d] = %d, offset = %ld", j, line, offset); 00077 type = Vect_read_line(Map, Points, NULL, line); 00078 if (lines[j] > 0) 00079 direction = GV_FORWARD; 00080 else 00081 direction = GV_BACKWARD; 00082 Vect_append_points(APoints, Points, direction); 00083 } 00084 00085 dig_find_area_poly(APoints, &area_size); 00086 G_debug(3, " area/isle size = %f", area_size); 00087 00088 if (area_size > 0) { /* area */ 00089 /* add area structure to plus */ 00090 area = dig_add_area(plus, n_lines, lines); 00091 if (area == -1) { /* error */ 00092 Vect_close(Map); 00093 G_fatal_error(_("Unable to add area (map closed, topo saved)")); 00094 } 00095 G_debug(3, " -> area %d", area); 00096 return area; 00097 } 00098 else if (area_size < 0) { /* island */ 00099 isle = dig_add_isle(plus, n_lines, lines); 00100 if (isle == -1) { /* error */ 00101 Vect_close(Map); 00102 G_fatal_error(_("Unable to add isle (map closed, topo saved)")); 00103 } 00104 G_debug(3, " -> isle %d", isle); 00105 return -isle; 00106 } 00107 else { 00108 /* TODO: What to do with such areas? Should be areas/isles of size 0 stored, 00109 * so that may be found and cleaned by some utility 00110 * Note: it would be useful for vertical closed polygons, but such would be added twice 00111 * as area */ 00112 G_warning(_("Area of size = 0.0 ignored")); 00113 } 00114 return 0; 00115 } 00116 00126 int Vect_isle_find_area(struct Map_info *Map, int isle) 00127 { 00128 int j, line, node, sel_area, first, area, poly; 00129 static int first_call = 1; 00130 struct Plus_head *plus; 00131 P_LINE *Line; 00132 P_NODE *Node; 00133 P_ISLE *Isle; 00134 P_AREA *Area; 00135 double size, cur_size; 00136 BOUND_BOX box, abox; 00137 static struct ilist *List; 00138 static struct line_pnts *APoints; 00139 00140 /* Note: We should check all isle points (at least) because if topology is not clean 00141 * and two areas overlap, isle which is not completely within area may be attached, 00142 * but it would take long time */ 00143 00144 G_debug(3, "Vect_isle_find_area () island = %d", isle); 00145 plus = &(Map->plus); 00146 00147 if (plus->Isle[isle] == NULL) { 00148 G_warning(_("Request to find area outside nonexistent isle")); 00149 return 0; 00150 } 00151 00152 if (first_call) { 00153 List = Vect_new_list(); 00154 APoints = Vect_new_line_struct(); 00155 first_call = 0; 00156 } 00157 00158 Isle = plus->Isle[isle]; 00159 line = abs(Isle->lines[0]); 00160 Line = plus->Line[line]; 00161 node = Line->N1; 00162 Node = plus->Node[node]; 00163 00164 /* select areas by box */ 00165 box.E = Node->x; 00166 box.W = Node->x; 00167 box.N = Node->y; 00168 box.S = Node->y; 00169 box.T = PORT_DOUBLE_MAX; 00170 box.B = -PORT_DOUBLE_MAX; 00171 Vect_select_areas_by_box(Map, &box, List); 00172 G_debug(3, "%d areas overlap island boundary point", List->n_values); 00173 00174 sel_area = 0; 00175 cur_size = -1; 00176 first = 1; 00177 Vect_get_isle_box(Map, isle, &box); 00178 for (j = 0; j < List->n_values; j++) { 00179 area = List->value[j]; 00180 G_debug(3, "area = %d", area); 00181 00182 Area = plus->Area[area]; 00183 00184 /* Before other tests, simply exclude those areas inside isolated isles formed by one boundary */ 00185 if (abs(Isle->lines[0]) == abs(Area->lines[0])) { 00186 G_debug(3, " area inside isolated isle"); 00187 continue; 00188 } 00189 00190 /* Check box */ 00191 /* Note: If build is run on large files of areas imported from nontopo format (shapefile) 00192 * attaching of isles takes very long time because each area is also isle and select by 00193 * box all overlapping areas selects all areas with box overlapping first node. 00194 * Then reading coordinates for all those areas would take a long time -> check first 00195 * if isle's box is completely within area box */ 00196 Vect_get_area_box(Map, area, &abox); 00197 if (box.E > abox.E || box.W < abox.W || box.N > abox.N || 00198 box.S < abox.S) { 00199 G_debug(3, " isle not completely inside area box"); 00200 continue; 00201 } 00202 00203 poly = Vect_point_in_area_outer_ring(Node->x, Node->y, Map, area); 00204 G_debug(3, " poly = %d", poly); 00205 00206 if (poly == 1) { /* pint in area, but node is not part of area inside isle (would be poly == 2) */ 00207 /* In rare case island is inside more areas in that case we have to calculate area 00208 * of outer ring and take the smaller */ 00209 if (sel_area == 0) { /* first */ 00210 sel_area = area; 00211 } 00212 else { /* is not first */ 00213 if (cur_size < 0) { /* second area */ 00214 /* This is slow, but should not be called often */ 00215 Vect_get_area_points(Map, sel_area, APoints); 00216 G_begin_polygon_area_calculations(); 00217 cur_size = 00218 G_area_of_polygon(APoints->x, APoints->y, 00219 APoints->n_points); 00220 G_debug(3, " first area size = %f (n points = %d)", 00221 cur_size, APoints->n_points); 00222 00223 } 00224 00225 Vect_get_area_points(Map, area, APoints); 00226 size = 00227 G_area_of_polygon(APoints->x, APoints->y, 00228 APoints->n_points); 00229 G_debug(3, " area size = %f (n points = %d)", cur_size, 00230 APoints->n_points); 00231 00232 if (size < cur_size) { 00233 sel_area = area; 00234 cur_size = size; 00235 } 00236 } 00237 G_debug(3, "sel_area = %d cur_size = %f", sel_area, cur_size); 00238 } 00239 } 00240 if (sel_area > 0) { 00241 G_debug(3, "Island %d in area %d", isle, sel_area); 00242 } 00243 else { 00244 G_debug(3, "Island %d is not in area", isle); 00245 } 00246 00247 return sel_area; 00248 } 00249 00258 int Vect_attach_isle(struct Map_info *Map, int isle) 00259 { 00260 int sel_area; 00261 P_ISLE *Isle; 00262 struct Plus_head *plus; 00263 00264 /* Note!: If topology is not clean and areas overlap, one island may fall to more areas 00265 * (partially or fully). Before isle is attached to area it must be check if it is not attached yet */ 00266 G_debug(3, "Vect_attach_isle (): isle = %d", isle); 00267 00268 plus = &(Map->plus); 00269 00270 sel_area = Vect_isle_find_area(Map, isle); 00271 G_debug(3, " isle = %d -> area outside = %d", isle, sel_area); 00272 if (sel_area > 0) { 00273 Isle = plus->Isle[isle]; 00274 if (Isle->area > 0) { 00275 G_debug(3, 00276 "Attempt to attach isle %d to more areas (=>topology is not clean)", 00277 isle); 00278 } 00279 else { 00280 Isle->area = sel_area; 00281 dig_area_add_isle(plus, sel_area, isle); 00282 } 00283 } 00284 return 0; 00285 } 00286 00295 int Vect_attach_isles(struct Map_info *Map, BOUND_BOX * box) 00296 { 00297 int i, isle; 00298 static int first = 1; 00299 static struct ilist *List; 00300 struct Plus_head *plus; 00301 00302 G_debug(3, "Vect_attach_isles ()"); 00303 00304 plus = &(Map->plus); 00305 00306 if (first) { 00307 List = Vect_new_list(); 00308 first = 0; 00309 } 00310 00311 Vect_select_isles_by_box(Map, box, List); 00312 G_debug(3, " number of isles to attach = %d", List->n_values); 00313 00314 for (i = 0; i < List->n_values; i++) { 00315 isle = List->value[i]; 00316 Vect_attach_isle(Map, isle); 00317 } 00318 return 0; 00319 } 00320 00329 int Vect_attach_centroids(struct Map_info *Map, BOUND_BOX * box) 00330 { 00331 int i, sel_area, centr; 00332 static int first = 1; 00333 static struct ilist *List; 00334 P_AREA *Area; 00335 P_LINE *Line; 00336 struct Plus_head *plus; 00337 00338 G_debug(3, "Vect_attach_centroids ()"); 00339 00340 plus = &(Map->plus); 00341 00342 if (first) { 00343 List = Vect_new_list(); 00344 first = 0; 00345 } 00346 00347 /* Warning: If map is updated on level2, it may happen that previously correct island 00348 * becomes incorrect. In that case, centroid of area forming the island is reattached 00349 * to outer area, because island polygon is not excluded. 00350 * 00351 * +-----------+ +-----------+ 00352 * | 1 | | 1 | 00353 * | +---+---+ | | +---+---+ | 00354 * | | 2 | 3 | | | | 2 | | 00355 * | | x | | | -> | | x | | 00356 * | | | | | | | | | 00357 * | +---+---+ | | +---+---+ | 00358 * | | | | 00359 * +-----------+ +-----------+ 00360 * centroid is centroid is 00361 * attached to 2 reattached to 1 00362 * 00363 * Because of this, when the centroid is reattached to another area, it is always necessary 00364 * to check if original area exist, unregister centroid from previous area. 00365 * To simplify code, this is implemented so that centroid is always firs unregistered 00366 * and if new area is found, it is registered again. 00367 * 00368 * This problem can be avoided altogether if properly attached centroids 00369 * are skipped 00370 * MM 2009 00371 */ 00372 00373 Vect_select_lines_by_box(Map, box, GV_CENTROID, List); 00374 G_debug(3, " number of centroids to reattach = %d", List->n_values); 00375 for (i = 0; i < List->n_values; i++) { 00376 int orig_area; 00377 00378 centr = List->value[i]; 00379 Line = plus->Line[centr]; 00380 00381 /* only attach unregistered and duplicate centroids because 00382 * 1) all properly attached centroids are properly attached, really! Don't touch. 00383 * 2) Vect_find_area() below does not always return the correct area 00384 * 3) it's faster 00385 */ 00386 if (Line->left > 0) 00387 continue; 00388 00389 orig_area = Line->left; 00390 00391 sel_area = Vect_find_area(Map, Line->E, Line->N); 00392 G_debug(3, " centroid %d is in area %d", centr, sel_area); 00393 if (sel_area > 0) { 00394 Area = plus->Area[sel_area]; 00395 if (Area->centroid == 0) { /* first centroid */ 00396 G_debug(3, " first centroid -> attach to area"); 00397 Area->centroid = centr; 00398 Line->left = sel_area; 00399 00400 if (sel_area != orig_area && plus->do_uplist) 00401 dig_line_add_updated(plus, centr); 00402 } 00403 else if (Area->centroid != centr) { /* duplicate centroid */ 00404 /* Note: it cannot happen that Area->centroid == centr, because the centroid 00405 * was not registered or a duplicate */ 00406 G_debug(3, " duplicate centroid -> do not attach to area"); 00407 Line->left = -sel_area; 00408 00409 if (-sel_area != orig_area && plus->do_uplist) 00410 dig_line_add_updated(plus, centr); 00411 } 00412 } 00413 00414 if (sel_area != orig_area && plus->do_uplist) 00415 dig_line_add_updated(plus, centr); 00416 } 00417 00418 return 0; 00419 } 00420 00430 int Vect_build_nat(struct Map_info *Map, int build) 00431 { 00432 struct Plus_head *plus; 00433 int i, s, type, lineid; 00434 long offset; 00435 int side, line, area; 00436 struct line_pnts *Points, *APoints; 00437 struct line_cats *Cats; 00438 P_LINE *Line; 00439 P_NODE *Node; 00440 P_AREA *Area; 00441 BOUND_BOX box; 00442 struct ilist *List; 00443 00444 G_debug(3, "Vect_build_nat() build = %d", build); 00445 00446 plus = &(Map->plus); 00447 00448 if (build == plus->built) 00449 return 1; /* Do nothing */ 00450 00451 /* Check if upgrade or downgrade */ 00452 if (build < plus->built) { /* lower level request */ 00453 00454 /* release old sources (this also initializes structures and numbers of elements) */ 00455 if (plus->built >= GV_BUILD_CENTROIDS && build < GV_BUILD_CENTROIDS) { 00456 /* reset info about areas stored for centroids */ 00457 int nlines = Vect_get_num_lines(Map); 00458 00459 for (line = 1; line <= nlines; line++) { 00460 Line = plus->Line[line]; 00461 if (Line && Line->type == GV_CENTROID) 00462 Line->left = 0; 00463 } 00464 dig_free_plus_areas(plus); 00465 dig_spidx_free_areas(plus); 00466 dig_free_plus_isles(plus); 00467 dig_spidx_free_isles(plus); 00468 } 00469 00470 00471 if (plus->built >= GV_BUILD_AREAS && build < GV_BUILD_AREAS) { 00472 /* reset info about areas stored for lines */ 00473 int nlines = Vect_get_num_lines(Map); 00474 00475 for (line = 1; line <= nlines; line++) { 00476 Line = plus->Line[line]; 00477 if (Line && Line->type == GV_BOUNDARY) { 00478 Line->left = 0; 00479 Line->right = 0; 00480 } 00481 } 00482 dig_free_plus_areas(plus); 00483 dig_spidx_free_areas(plus); 00484 dig_free_plus_isles(plus); 00485 dig_spidx_free_isles(plus); 00486 } 00487 if (plus->built >= GV_BUILD_BASE && build < GV_BUILD_BASE) { 00488 dig_free_plus_nodes(plus); 00489 dig_spidx_free_nodes(plus); 00490 dig_free_plus_lines(plus); 00491 dig_spidx_free_lines(plus); 00492 } 00493 00494 plus->built = build; 00495 return 1; 00496 } 00497 00498 Points = Vect_new_line_struct(); 00499 APoints = Vect_new_line_struct(); 00500 Cats = Vect_new_cats_struct(); 00501 List = Vect_new_list(); 00502 00503 if (plus->built < GV_BUILD_BASE) { 00504 int npoints, format; 00505 00506 format = G_info_format(); 00507 00508 /* 00509 * We shall go through all primitives in coor file and 00510 * add new node for each end point to nodes structure 00511 * if the node with the same coordinates doesn't exist yet. 00512 */ 00513 00514 /* register lines, create nodes */ 00515 Vect_rewind(Map); 00516 G_message(_("Registering primitives...")); 00517 i = 1; 00518 npoints = 0; 00519 while (1) { 00520 /* register line */ 00521 type = Vect_read_next_line(Map, Points, Cats); 00522 00523 /* Note: check for dead lines is not needed, because they are skipped by V1_read_next_line_nat() */ 00524 if (type == -1) { 00525 G_warning(_("Unable to read vector map")); 00526 return 0; 00527 } 00528 else if (type == -2) { 00529 break; 00530 } 00531 00532 npoints += Points->n_points; 00533 00534 offset = Map->head.last_offset; 00535 00536 G_debug(3, "Register line: offset = %ld", offset); 00537 lineid = dig_add_line(plus, type, Points, offset); 00538 dig_line_box(Points, &box); 00539 if (lineid == 1) 00540 Vect_box_copy(&(plus->box), &box); 00541 else 00542 Vect_box_extend(&(plus->box), &box); 00543 00544 /* Add all categories to category index */ 00545 if (build == GV_BUILD_ALL) { 00546 int c; 00547 00548 for (c = 0; c < Cats->n_cats; c++) { 00549 dig_cidx_add_cat(plus, Cats->field[c], Cats->cat[c], 00550 lineid, type); 00551 } 00552 if (Cats->n_cats == 0) /* add field 0, cat 0 */ 00553 dig_cidx_add_cat(plus, 0, 0, lineid, type); 00554 } 00555 00556 if (G_verbose() > G_verbose_min() && i % 1000 == 0) { 00557 if (format == G_INFO_FORMAT_PLAIN) 00558 fprintf(stderr, "%d..", i); 00559 else 00560 fprintf(stderr, "%9d\b\b\b\b\b\b\b\b\b", i); 00561 } 00562 00563 i++; 00564 } 00565 00566 if ( (G_verbose() > G_verbose_min() ) && format != G_INFO_FORMAT_PLAIN ) 00567 fprintf(stderr, "\r"); 00568 00569 G_message(_("%d primitives registered"), plus->n_lines); 00570 G_message(_("%d vertices registered"), npoints); 00571 00572 plus->built = GV_BUILD_BASE; 00573 } 00574 00575 if (build < GV_BUILD_AREAS) 00576 return 1; 00577 00578 if (plus->built < GV_BUILD_AREAS) { 00579 /* Build areas */ 00580 /* Go through all bundaries and try to build area for both sides */ 00581 G_important_message(_("Building areas...")); 00582 for (i = 1; i <= plus->n_lines; i++) { 00583 G_percent(i, plus->n_lines, 1); 00584 00585 /* build */ 00586 if (plus->Line[i] == NULL) { 00587 continue; 00588 } /* dead line */ 00589 Line = plus->Line[i]; 00590 if (Line->type != GV_BOUNDARY) { 00591 continue; 00592 } 00593 00594 for (s = 0; s < 2; s++) { 00595 if (s == 0) 00596 side = GV_LEFT; 00597 else 00598 side = GV_RIGHT; 00599 00600 G_debug(3, "Build area for line = %d, side = %d", i, side); 00601 Vect_build_line_area(Map, i, side); 00602 } 00603 } 00604 G_message(_("%d areas built"), plus->n_areas); 00605 G_message(_("%d isles built"), plus->n_isles); 00606 plus->built = GV_BUILD_AREAS; 00607 } 00608 00609 if (build < GV_BUILD_ATTACH_ISLES) 00610 return 1; 00611 00612 /* Attach isles to areas */ 00613 if (plus->built < GV_BUILD_ATTACH_ISLES) { 00614 G_important_message(_("Attaching islands...")); 00615 for (i = 1; i <= plus->n_isles; i++) { 00616 G_percent(i, plus->n_isles, 1); 00617 Vect_attach_isle(Map, i); 00618 } 00619 plus->built = GV_BUILD_ATTACH_ISLES; 00620 } 00621 00622 if (build < GV_BUILD_CENTROIDS) 00623 return 1; 00624 00625 /* Attach centroids to areas */ 00626 if (plus->built < GV_BUILD_CENTROIDS) { 00627 int nlines; 00628 00629 G_important_message(_("Attaching centroids...")); 00630 00631 nlines = Vect_get_num_lines(Map); 00632 for (line = 1; line <= nlines; line++) { 00633 G_percent(line, nlines, 1); 00634 00635 Line = plus->Line[line]; 00636 if (!Line) 00637 continue; /* Dead */ 00638 00639 if (Line->type != GV_CENTROID) 00640 continue; 00641 00642 Node = plus->Node[Line->N1]; 00643 00644 area = Vect_find_area(Map, Node->x, Node->y); 00645 00646 if (area > 0) { 00647 G_debug(3, "Centroid (line=%d) in area %d", line, area); 00648 00649 Area = plus->Area[area]; 00650 00651 if (Area->centroid == 0) { /* first */ 00652 Area->centroid = line; 00653 Line->left = area; 00654 } 00655 else { /* duplicate */ 00656 Line->left = -area; 00657 } 00658 } 00659 } 00660 plus->built = GV_BUILD_CENTROIDS; 00661 } 00662 00663 /* Add areas to category index */ 00664 for (area = 1; area <= plus->n_areas; area++) { 00665 int c; 00666 00667 if (plus->Area[area] == NULL) 00668 continue; 00669 00670 if (plus->Area[area]->centroid > 0) { 00671 Vect_read_line(Map, NULL, Cats, plus->Area[area]->centroid); 00672 00673 for (c = 0; c < Cats->n_cats; c++) { 00674 dig_cidx_add_cat(plus, Cats->field[c], Cats->cat[c], area, 00675 GV_AREA); 00676 } 00677 } 00678 00679 if (plus->Area[area]->centroid == 0 || Cats->n_cats == 0) /* no centroid or no cats */ 00680 dig_cidx_add_cat(plus, 0, 0, area, GV_AREA); 00681 } 00682 00683 return 1; 00684 }