GRASS Programmer's Manual 6.4.1(2011)
|
00001 00020 #include <stdlib.h> 00021 #include <stdio.h> 00022 #include <stdarg.h> 00023 #include <grass/glocale.h> 00024 #include <grass/gis.h> 00025 #include <grass/Vect.h> 00026 00027 00028 #ifndef HAVE_OGR 00029 static int format() 00030 { 00031 G_fatal_error(_("Requested format is not compiled in this version")); 00032 return 0; 00033 } 00034 #endif 00035 00036 static int (*Build_array[]) () = { 00037 Vect_build_nat 00038 #ifdef HAVE_OGR 00039 , Vect_build_ogr 00040 #else 00041 , format 00042 #endif 00043 }; 00044 00053 int Vect_build(struct Map_info *Map) 00054 { 00055 return Vect_build_partial(Map, GV_BUILD_ALL); 00056 } 00057 00058 00066 int Vect_get_built(struct Map_info *Map) 00067 { 00068 return (Map->plus.built); 00069 } 00070 00107 int Vect_build_partial(struct Map_info *Map, int build) 00108 { 00109 struct Plus_head *plus; 00110 int ret; 00111 00112 G_debug(3, "Vect_build(): build = %d", build); 00113 00114 /* If topology is already build (map on level2), set level to 1 so that lines will 00115 * be read by V1_read_ (all lines) */ 00116 Map->level = 1; /* may be not needed, because V1_read is used directly by Vect_build_ */ 00117 Map->support_updated = 1; 00118 00119 Map->plus.Spidx_built = 1; 00120 00121 plus = &(Map->plus); 00122 if (build > GV_BUILD_NONE) { 00123 G_message(_("Building topology for vector map <%s>..."), 00124 Vect_get_name(Map)); 00125 } 00126 plus->with_z = Map->head.with_z; 00127 plus->spidx_with_z = Map->head.with_z; 00128 00129 if (build == GV_BUILD_ALL) { 00130 dig_cidx_free(plus); /* free old (if any) category index) */ 00131 dig_cidx_init(plus); 00132 } 00133 00134 ret = ((*Build_array[Map->format]) (Map, build)); 00135 00136 if (ret == 0) { 00137 return 0; 00138 } 00139 00140 if (build > GV_BUILD_NONE) { 00141 G_verbose_message(_("Topology was built")); 00142 } 00143 00144 Map->level = LEVEL_2; 00145 plus->mode = GV_MODE_WRITE; 00146 00147 if (build == GV_BUILD_ALL) { 00148 plus->cidx_up_to_date = 1; /* category index was build */ 00149 dig_cidx_sort(plus); 00150 } 00151 00152 if (build > GV_BUILD_NONE) { 00153 G_message(_("Number of nodes: %d"), plus->n_nodes); 00154 G_message(_("Number of primitives: %d"), plus->n_lines); 00155 G_message(_("Number of points: %d"), plus->n_plines); 00156 G_message(_("Number of lines: %d"), plus->n_llines); 00157 G_message(_("Number of boundaries: %d"), plus->n_blines); 00158 G_message(_("Number of centroids: %d"), plus->n_clines); 00159 00160 if (plus->n_flines > 0) 00161 G_message(_("Number of faces: %d"), plus->n_flines); 00162 00163 if (plus->n_klines > 0) 00164 G_message(_("Number of kernels: %d"), plus->n_klines); 00165 } 00166 00167 if (plus->built >= GV_BUILD_AREAS) { 00168 int line, nlines, area, nareas, err_boundaries, err_centr_out, 00169 err_centr_dupl, err_nocentr; 00170 P_LINE *Line; 00171 struct Plus_head *Plus; 00172 00173 /* Count errors (it does not take much time comparing to build process) */ 00174 Plus = &(Map->plus); 00175 nlines = Vect_get_num_lines(Map); 00176 err_boundaries = err_centr_out = err_centr_dupl = 0; 00177 for (line = 1; line <= nlines; line++) { 00178 Line = Plus->Line[line]; 00179 if (!Line) 00180 continue; 00181 if (Line->type == GV_BOUNDARY && 00182 (Line->left == 0 || Line->right == 0)) { 00183 G_debug(3, "line = %d left = %d right = %d", line, Line->left, 00184 Line->right); 00185 err_boundaries++; 00186 } 00187 if (Line->type == GV_CENTROID) { 00188 if (Line->left == 0) 00189 err_centr_out++; 00190 else if (Line->left < 0) 00191 err_centr_dupl++; 00192 } 00193 } 00194 00195 err_nocentr = 0; 00196 nareas = Vect_get_num_areas(Map); 00197 for (area = 1; area <= nareas; area++) { 00198 if (!Vect_area_alive(Map, area)) 00199 continue; 00200 line = Vect_get_area_centroid(Map, area); 00201 if (line == 0) 00202 err_nocentr++; 00203 } 00204 00205 G_message(_("Number of areas: %d"), plus->n_areas); 00206 G_message(_("Number of isles: %d"), plus->n_isles); 00207 00208 if (err_boundaries) 00209 G_message(_("Number of incorrect boundaries: %d"), 00210 err_boundaries); 00211 00212 if (err_centr_out) 00213 G_message(_("Number of centroids outside area: %d"), 00214 err_centr_out); 00215 00216 if (err_centr_dupl) 00217 G_message(_("Number of duplicate centroids: %d"), 00218 err_centr_dupl); 00219 00220 if (err_nocentr) 00221 G_message(_("Number of areas without centroid: %d"), 00222 err_nocentr); 00223 00224 } 00225 else if (build > GV_BUILD_NONE) { 00226 G_message(_("Number of areas: -")); 00227 G_message(_("Number of isles: -")); 00228 } 00229 return 1; 00230 } 00231 00239 int Vect_save_topo(struct Map_info *Map) 00240 { 00241 struct Plus_head *plus; 00242 char fname[GPATH_MAX], buf[GPATH_MAX]; 00243 GVFILE fp; 00244 00245 G_debug(1, "Vect_save_topo()"); 00246 00247 plus = &(Map->plus); 00248 00249 /* write out all the accumulated info to the plus file */ 00250 sprintf(buf, "%s/%s", GRASS_VECT_DIRECTORY, Map->name); 00251 G__file_name(fname, buf, GV_TOPO_ELEMENT, Map->mapset); 00252 G_debug(1, "Open topo: %s", fname); 00253 dig_file_init(&fp); 00254 fp.file = fopen(fname, "w"); 00255 if (fp.file == NULL) { 00256 G_warning(_("Unable to open topo file for write <%s>"), fname); 00257 return 0; 00258 } 00259 00260 /* set portable info */ 00261 dig_init_portable(&(plus->port), dig__byte_order_out()); 00262 00263 if (0 > dig_write_plus_file(&fp, plus)) { 00264 G_warning(_("Error writing out topo file")); 00265 return 0; 00266 } 00267 00268 fclose(fp.file); 00269 00270 return 1; 00271 } 00272 00282 int Vect_topo_dump(struct Map_info *Map, FILE * out) 00283 { 00284 int i, j, line, isle; 00285 P_NODE *Node; 00286 P_LINE *Line; 00287 P_AREA *Area; 00288 P_ISLE *Isle; 00289 BOUND_BOX box; 00290 struct Plus_head *plus; 00291 00292 plus = &(Map->plus); 00293 00294 fprintf(out, "---------- TOPOLOGY DUMP ----------\n"); 00295 00296 /* box */ 00297 Vect_box_copy(&box, &(plus->box)); 00298 fprintf(out, "N,S,E,W,T,B: %f, %f, %f, %f, %f, %f\n", box.N, box.S, 00299 box.E, box.W, box.T, box.B); 00300 00301 /* nodes */ 00302 fprintf(out, "Nodes (%d nodes, alive + dead ):\n", plus->n_nodes); 00303 for (i = 1; i <= plus->n_nodes; i++) { 00304 if (plus->Node[i] == NULL) { 00305 continue; 00306 } 00307 Node = plus->Node[i]; 00308 fprintf(out, "node = %d, n_lines = %d, xyz = %f, %f, %f\n", i, 00309 Node->n_lines, Node->x, Node->y, Node->z); 00310 for (j = 0; j < Node->n_lines; j++) { 00311 line = Node->lines[j]; 00312 Line = plus->Line[abs(line)]; 00313 fprintf(out, " line = %3d, type = %d, angle = %f\n", line, 00314 Line->type, Node->angles[j]); 00315 } 00316 } 00317 00318 /* lines */ 00319 fprintf(out, "Lines (%d lines, alive + dead ):\n", plus->n_lines); 00320 for (i = 1; i <= plus->n_lines; i++) { 00321 if (plus->Line[i] == NULL) { 00322 continue; 00323 } 00324 Line = plus->Line[i]; 00325 fprintf(out, "line = %d, type = %d, offset = %ld n1 = %d, n2 = %d, " 00326 "left/area = %d, right = %d\n", 00327 i, Line->type, Line->offset, Line->N1, Line->N2, 00328 Line->left, Line->right); 00329 fprintf(out, "N,S,E,W,T,B: %f, %f, %f, %f, %f, %f\n", Line->N, 00330 Line->S, Line->E, Line->W, Line->T, Line->B); 00331 } 00332 00333 /* areas */ 00334 fprintf(out, "Areas (%d areas, alive + dead ):\n", plus->n_areas); 00335 for (i = 1; i <= plus->n_areas; i++) { 00336 if (plus->Area[i] == NULL) { 00337 continue; 00338 } 00339 Area = plus->Area[i]; 00340 00341 fprintf(out, "area = %d, n_lines = %d, n_isles = %d centroid = %d\n", 00342 i, Area->n_lines, Area->n_isles, Area->centroid); 00343 00344 fprintf(out, "N,S,E,W,T,B: %f, %f, %f, %f, %f, %f\n", Area->N, 00345 Area->S, Area->E, Area->W, Area->T, Area->B); 00346 00347 for (j = 0; j < Area->n_lines; j++) { 00348 line = Area->lines[j]; 00349 Line = plus->Line[abs(line)]; 00350 fprintf(out, " line = %3d\n", line); 00351 } 00352 for (j = 0; j < Area->n_isles; j++) { 00353 isle = Area->isles[j]; 00354 fprintf(out, " isle = %3d\n", isle); 00355 } 00356 } 00357 00358 /* isles */ 00359 fprintf(out, "Islands (%d islands, alive + dead ):\n", plus->n_isles); 00360 for (i = 1; i <= plus->n_isles; i++) { 00361 if (plus->Isle[i] == NULL) { 00362 continue; 00363 } 00364 Isle = plus->Isle[i]; 00365 00366 fprintf(out, "isle = %d, n_lines = %d area = %d\n", i, Isle->n_lines, 00367 Isle->area); 00368 00369 fprintf(out, "N,S,E,W,T,B: %f, %f, %f, %f, %f, %f\n", Isle->N, 00370 Isle->S, Isle->E, Isle->W, Isle->T, Isle->B); 00371 00372 for (j = 0; j < Isle->n_lines; j++) { 00373 line = Isle->lines[j]; 00374 Line = plus->Line[abs(line)]; 00375 fprintf(out, " line = %3d\n", line); 00376 } 00377 } 00378 00379 return 1; 00380 } 00381 00390 int Vect_save_spatial_index(struct Map_info *Map) 00391 { 00392 struct Plus_head *plus; 00393 char fname[GPATH_MAX], buf[GPATH_MAX]; 00394 GVFILE fp; 00395 00396 G_debug(1, "Vect_save_spatial_index()"); 00397 00398 plus = &(Map->plus); 00399 00400 /* write out rtrees to the sidx file */ 00401 sprintf(buf, "%s/%s", GRASS_VECT_DIRECTORY, Map->name); 00402 G__file_name(fname, buf, GV_SIDX_ELEMENT, Map->mapset); 00403 G_debug(1, "Open sidx: %s", fname); 00404 dig_file_init(&fp); 00405 fp.file = fopen(fname, "w"); 00406 if (fp.file == NULL) { 00407 G_warning(_("Unable open spatial index file for write <%s>"), fname); 00408 return 0; 00409 } 00410 00411 /* set portable info */ 00412 dig_init_portable(&(plus->spidx_port), dig__byte_order_out()); 00413 00414 if (0 > dig_write_spidx(&fp, plus)) { 00415 G_warning(_("Error writing out spatial index file")); 00416 return 0; 00417 } 00418 00419 fclose(fp.file); 00420 00421 return 1; 00422 } 00423 00433 int Vect_spatial_index_dump(struct Map_info *Map, FILE * out) 00434 { 00435 if (!(Map->plus.Spidx_built)) { 00436 Vect_build_sidx_from_topo(Map); 00437 } 00438 00439 fprintf(out, "---------- SPATIAL INDEX DUMP ----------\n"); 00440 00441 dig_dump_spidx(out, &(Map->plus)); 00442 00443 return 1; 00444 }