GRASS Programmer's Manual
6.4.1(2011)
|
00001 00016 #include <grass/config.h> 00017 #include <grass/gis.h> 00018 #include <grass/Vect.h> 00019 #include <grass/glocale.h> 00020 00021 #ifdef HAVE_GEOS 00022 00023 static GEOSGeometry *Vect__read_line_geos(struct Map_info *, long, int *); 00024 static GEOSCoordSequence *V1_read_line_geos(struct Map_info *, long, int *); 00025 static GEOSCoordSequence *V2_read_line_geos(struct Map_info *, int); 00026 static GEOSCoordSequence *read_polygon_points(struct Map_info *, int, int*); 00027 00046 GEOSGeometry *Vect_read_line_geos(struct Map_info *Map, int line, int *type) 00047 { 00048 P_LINE *Line; 00049 00050 G_debug(3, "Vect_read_line_geos(): line = %d", line); 00051 00052 if (!VECT_OPEN(Map)) 00053 G_fatal_error("Vect_read_line_geos(): %s", _("vector map is not opened")); 00054 00055 if (line < 1 || line > Map->plus.n_lines) 00056 G_fatal_error(_("Vect_read_line_geos(): feature id %d is not reasonable " 00057 "(max features in vector map <%s>: %d)"), 00058 line, Vect_get_full_name(Map), Map->plus.n_lines); 00059 00060 if (Map->format != GV_FORMAT_NATIVE) 00061 G_fatal_error("Vect_read_line_geos(): %s", _("only native format supported")); 00062 00063 Line = Map->plus.Line[line]; 00064 if (Line == NULL) 00065 G_fatal_error("Vect_read_line_geos(): %s %d", 00066 _("Attempt to read dead line"), line); 00067 00068 return Vect__read_line_geos(Map, Line->offset, type); 00069 } 00070 00082 GEOSGeometry *Vect_read_area_geos(struct Map_info * Map, int area) 00083 { 00084 int i, nholes, isle; 00085 GEOSGeometry *boundary, **holes; 00086 00087 G_debug(3, "Vect_read_area_geos(): area = %d", area); 00088 00089 boundary = GEOSGeom_createLinearRing(Vect_get_area_points_geos(Map, area)); 00090 if (!boundary) { 00091 G_fatal_error(_("Vect_read_area_geos(): unable to read area id %d"), 00092 area); 00093 } 00094 00095 nholes = Vect_get_area_num_isles(Map, area); 00096 holes = (GEOSGeometry **) G_malloc(nholes * sizeof(GEOSGeometry *)); 00097 for (i = 0; i < nholes; i++) { 00098 isle = Vect_get_area_isle(Map, area, i); 00099 if (isle < 1) { 00100 nholes--; 00101 continue; 00102 } 00103 holes[i] = GEOSGeom_createLinearRing(Vect_get_isle_points_geos(Map, isle)); 00104 if (!(holes[i])) 00105 G_fatal_error(_("Vect_read_area_geos(): unable to read isle id %d of area id %d"), 00106 isle, area); 00107 } 00108 00109 return GEOSGeom_createPolygon(boundary, holes, nholes); 00110 } 00111 00128 GEOSGeometry *Vect_line_to_geos(struct Map_info *Map, 00129 const struct line_pnts *points, int type) 00130 { 00131 int i, with_z; 00132 GEOSGeometry *geom; 00133 GEOSCoordSequence *pseq; 00134 00135 G_debug(3, "Vect_line_to_geos(): type = %d", type); 00136 00137 with_z = Vect_is_3d(Map); 00138 00139 /* read only points / lines / boundaries */ 00140 if (!(type & (GV_POINT | GV_LINES))) 00141 return NULL; 00142 00143 if (type == GV_POINT) { 00144 if (points->n_points != 1) 00145 /* point is not valid */ 00146 return NULL; 00147 } 00148 else { 00149 if (points->n_points < 2) 00150 /* line/boundary is not valid */ 00151 return NULL; 00152 } 00153 00154 pseq = GEOSCoordSeq_create(points->n_points, with_z ? 3 : 2); 00155 00156 for (i = 0; i < points->n_points; i++) { 00157 GEOSCoordSeq_setX(pseq, i, points->x[i]); 00158 GEOSCoordSeq_setY(pseq, i, points->y[i]); 00159 if (with_z) 00160 GEOSCoordSeq_setZ(pseq, i, points->z[i]); 00161 } 00162 00163 if (type == GV_POINT) 00164 geom = GEOSGeom_createPoint(pseq); 00165 else if (type == GV_LINE) 00166 geom = GEOSGeom_createLineString(pseq); 00167 else { /* boundary */ 00168 geom = GEOSGeom_createLineString(pseq); 00169 if (GEOSisRing(geom)) { 00170 /* GEOSGeom_destroy(geom); */ 00171 geom = GEOSGeom_createLinearRing(pseq); 00172 } 00173 } 00174 00175 /* GEOSCoordSeq_destroy(pseq); */ 00176 00177 return geom; 00178 } 00179 00194 GEOSGeometry *Vect__read_line_geos(struct Map_info *Map, long offset, int *type) 00195 { 00196 int ftype; 00197 00198 GEOSGeometry *geom; 00199 GEOSCoordSequence *pseq; 00200 00201 pseq = V1_read_line_geos(Map, offset, &ftype); 00202 if (!pseq) 00203 G_fatal_error(_("Unable to read line offset %ld"), offset); 00204 00205 if (ftype & GV_POINT) { 00206 G_debug(3, " geos_type = point"); 00207 geom = GEOSGeom_createPoint(pseq); 00208 } 00209 else if (ftype & GV_LINE) { 00210 G_debug(3, " geos_type = linestring"); 00211 geom = GEOSGeom_createLineString(pseq); 00212 } 00213 else { /* boundary */ 00214 geom = GEOSGeom_createLineString(pseq); 00215 if (GEOSisRing(geom)) { 00216 /* GEOSGeom_destroy(geom); */ 00217 geom = GEOSGeom_createLinearRing(pseq); 00218 G_debug(3, " geos_type = linearring"); 00219 } 00220 else { 00221 G_debug(3, " geos_type = linestring"); 00222 } 00223 } 00224 00225 /* GEOSCoordSeq_destroy(pseq); */ 00226 00227 if (type) 00228 *type = ftype; 00229 00230 return geom; 00231 } 00232 00245 GEOSCoordSequence *V2_read_line_geos(struct Map_info *Map, int line) 00246 { 00247 int ftype; 00248 P_LINE *Line; 00249 00250 G_debug(3, "V2_read_line_geos(): line = %d", line); 00251 00252 Line = Map->plus.Line[line]; 00253 00254 if (Line == NULL) 00255 G_fatal_error("V2_read_line_geos(): %s %d", 00256 _("Attempt to read dead line"), line); 00257 00258 return V1_read_line_geos(Map, Line->offset, &ftype); 00259 } 00260 00261 00278 GEOSCoordSequence *V1_read_line_geos(struct Map_info *Map, long offset, int *type) 00279 { 00280 int i, n_points; 00281 int do_cats, n_cats; 00282 char rhead, nc; 00283 long size; 00284 double *x, *y, *z; 00285 00286 GEOSCoordSequence *pseq; 00287 00288 G_debug(3, "V1_read_line_geos(): offset = %ld", offset); 00289 00290 Map->head.last_offset = offset; 00291 00292 /* reads must set in_head, but writes use default */ 00293 dig_set_cur_port(&(Map->head.port)); 00294 00295 dig_fseek(&(Map->dig_fp), offset, 0); 00296 00297 if (0 >= dig__fread_port_C(&rhead, 1, &(Map->dig_fp))) 00298 return NULL; /* end of file */ 00299 00300 if (!(rhead & 0x01)) /* dead line */ 00301 return GEOSCoordSeq_create(0, (Map->head.with_z) ? 3 : 2); 00302 00303 if (rhead & 0x02) /* categories exists */ 00304 do_cats = 1; /* do not return here let file offset moves forward to next */ 00305 else /* line */ 00306 do_cats = 0; 00307 00308 rhead >>= 2; 00309 *type = dig_type_from_store((int) rhead); 00310 00311 /* read only points / lines / boundaries */ 00312 if (!(*type & (GV_POINT | GV_LINES))) 00313 return GEOSCoordSeq_create(0, (Map->head.with_z) ? 3 : 2); 00314 00315 /* skip categories */ 00316 if (do_cats) { 00317 if (Map->head.Version_Minor == 1) { /* coor format 5.1 */ 00318 if (0 >= dig__fread_port_I(&n_cats, 1, &(Map->dig_fp))) 00319 return NULL; 00320 } 00321 else { /* coor format 5.0 */ 00322 if (0 >= dig__fread_port_C(&nc, 1, &(Map->dig_fp))) 00323 return NULL; 00324 n_cats = (int) nc; 00325 } 00326 G_debug(3, " n_cats = %d", n_cats); 00327 00328 if (Map->head.Version_Minor == 1) { /* coor format 5.1 */ 00329 size = (2 * PORT_INT) * n_cats; 00330 } 00331 else { /* coor format 5.0 */ 00332 size = (PORT_SHORT + PORT_INT) * n_cats; 00333 } 00334 dig_fseek(&(Map->dig_fp), size, SEEK_CUR); 00335 } 00336 00337 if (*type & GV_POINTS) { 00338 n_points = 1; 00339 } 00340 else { 00341 if (0 >= dig__fread_port_I(&n_points, 1, &(Map->dig_fp))) 00342 return NULL; 00343 } 00344 00345 G_debug(3, " n_points = %d dim = %d", n_points, (Map->head.with_z) ? 3 : 2); 00346 00347 pseq = GEOSCoordSeq_create(n_points, (Map->head.with_z) ? 3 : 2); 00348 00349 x = (double *) G_malloc(n_points * sizeof(double)); 00350 y = (double *) G_malloc(n_points * sizeof(double)); 00351 if (Map->head.with_z) 00352 z = (double *) G_malloc(n_points * sizeof(double)); 00353 else 00354 z = NULL; 00355 00356 if (0 >= dig__fread_port_D(x, n_points, &(Map->dig_fp))) 00357 return NULL; /* end of file */ 00358 00359 if (0 >= dig__fread_port_D(y, n_points, &(Map->dig_fp))) 00360 return NULL; /* end of file */ 00361 00362 if (Map->head.with_z) { 00363 if (0 >= dig__fread_port_D(z, n_points, &(Map->dig_fp))) 00364 return NULL; /* end of file */ 00365 00366 } 00367 00368 for (i = 0; i < n_points; i++) { 00369 GEOSCoordSeq_setX(pseq, i, x[i]); 00370 GEOSCoordSeq_setY(pseq, i, y[i]); 00371 if (Map->head.with_z) 00372 GEOSCoordSeq_setZ(pseq, i, z[i]); 00373 } 00374 00375 G_debug(3, " off = %ld", dig_ftell(&(Map->dig_fp))); 00376 00377 G_free((void *) x); 00378 G_free((void *) y); 00379 if (z) 00380 G_free((void *) z); 00381 00382 return pseq; 00383 } 00384 00399 GEOSCoordSequence *Vect_get_area_points_geos(struct Map_info *Map, int area) 00400 { 00401 struct Plus_head *Plus; 00402 P_AREA *Area; 00403 00404 G_debug(3, "Vect_get_area_points_geos(): area = %d", area); 00405 00406 Plus = &(Map->plus); 00407 Area = Plus->Area[area]; 00408 00409 if (Area == NULL) { /* dead area */ 00410 G_warning(_("Attempt to read points of nonexistent area id %d"), area); 00411 return NULL; /* error , because we should not read dead areas */ 00412 } 00413 00414 return read_polygon_points(Map, Area->n_lines, Area->lines); 00415 } 00416 00430 GEOSCoordSequence *Vect_get_isle_points_geos(struct Map_info *Map, int isle) 00431 { 00432 struct Plus_head *Plus; 00433 P_ISLE *Isle; 00434 00435 G_debug(3, "Vect_get_isle_points_geos(): isle = %d", isle); 00436 00437 Plus = &(Map->plus); 00438 Isle = Plus->Isle[isle]; 00439 00440 return read_polygon_points(Map, Isle->n_lines, Isle->lines); 00441 } 00442 00443 GEOSCoordSequence *read_polygon_points(struct Map_info *Map, int n_lines, int *lines) 00444 { 00445 int i, j, k; 00446 int line, aline; 00447 unsigned int n_points, n_points_shell; 00448 double x, y, z; 00449 int *dir; 00450 00451 GEOSCoordSequence **pseq, *pseq_shell; 00452 00453 G_debug(3, " n_lines = %d", n_lines); 00454 pseq = (GEOSCoordSequence **) G_malloc(n_lines * sizeof(GEOSCoordSequence *)); 00455 dir = (int*) G_malloc(n_lines * sizeof(int)); 00456 00457 n_points_shell = 0; 00458 for (i = 0; i < n_lines; i++) { 00459 line = lines[i]; 00460 aline = abs(line); 00461 G_debug(3, " append line(%d) = %d", i, line); 00462 00463 if (line > 0) 00464 dir[i] = GV_FORWARD; 00465 else 00466 dir[i] = GV_BACKWARD; 00467 00468 pseq[i] = V2_read_line_geos(Map, aline); 00469 if (!(pseq[i])) { 00470 G_fatal_error(_("Unable to read feature id %d"), aline); 00471 } 00472 00473 GEOSCoordSeq_getSize(pseq[i], &n_points); 00474 G_debug(3, " line n_points = %d", n_points); 00475 n_points_shell += n_points; 00476 } 00477 00478 /* create shell (outer ring) */ 00479 pseq_shell = GEOSCoordSeq_create(n_points_shell, Map->head.with_z ? 3 : 2); 00480 k = 0; 00481 for (i = 0; i < n_lines; i++) { 00482 GEOSCoordSeq_getSize(pseq[i], &n_points); 00483 if (dir[i] == GV_FORWARD) { 00484 for (j = 0; j < (int) n_points; j++, k++) { 00485 GEOSCoordSeq_getX(pseq[i], j, &x); 00486 GEOSCoordSeq_setX(pseq_shell, k, x); 00487 00488 GEOSCoordSeq_getY(pseq[i], j, &y); 00489 GEOSCoordSeq_setY(pseq_shell, k, y); 00490 00491 if (Map->head.with_z) { 00492 GEOSCoordSeq_getY(pseq[i], j, &z); 00493 GEOSCoordSeq_setZ(pseq_shell, k, z); 00494 } 00495 } 00496 } 00497 else { /* GV_BACKWARD */ 00498 for (j = (int) n_points - 1; j > -1; j--, k++) { 00499 GEOSCoordSeq_getX(pseq[i], j, &x); 00500 GEOSCoordSeq_setX(pseq_shell, k, x); 00501 00502 GEOSCoordSeq_getY(pseq[i], j, &y); 00503 GEOSCoordSeq_setY(pseq_shell, k, y); 00504 00505 if (Map->head.with_z) { 00506 GEOSCoordSeq_getY(pseq[i], j, &z); 00507 GEOSCoordSeq_setZ(pseq_shell, k, z); 00508 } 00509 } 00510 } 00511 } 00512 00513 G_free((void *) pseq); 00514 G_free((void *) dir); 00515 00516 return pseq_shell; 00517 } 00518 #endif /* HAVE_GEOS */