GRASS Programmer's Manual 6.4.1(2011)
geos.c
Go to the documentation of this file.
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 */
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Defines