geos.c

Go to the documentation of this file.

#include <grass/config.h>
#include <grass/gis.h>
#include <grass/Vect.h>
#include <grass/glocale.h>

#ifdef HAVE_GEOS

static GEOSGeometry *Vect__read_line_geos(struct Map_info *, long, int *);
static GEOSCoordSequence *V1_read_line_geos(struct Map_info *, long, int *);
static GEOSCoordSequence *V2_read_line_geos(struct Map_info *, int);
static GEOSCoordSequence *read_polygon_points(struct Map_info *, int, int*);

GEOSGeometry *Vect_read_line_geos(struct Map_info *Map, int line, int *type)
{
    P_LINE *Line;
    
    G_debug(3, "Vect_read_line_geos(): line = %d", line);
    
    if (!VECT_OPEN(Map))
        G_fatal_error("Vect_read_line_geos(): %s", _("vector map is not opened"));
    
    if (line < 1 || line > Map->plus.n_lines)
        G_fatal_error(_("Vect_read_line_geos(): feature id %d is not reasonable "
                        "(max features in vector map <%s>: %d)"),
                      line, Vect_get_full_name(Map), Map->plus.n_lines);
    
    if (Map->format != GV_FORMAT_NATIVE)
        G_fatal_error("Vect_read_line_geos(): %s", _("only native format supported"));
    
    Line = Map->plus.Line[line];
    if (Line == NULL)
        G_fatal_error("Vect_read_line_geos(): %s %d",
                      _("Attempt to read dead line"), line);
    
    return Vect__read_line_geos(Map, Line->offset, type);
}

GEOSGeometry *Vect_read_area_geos(struct Map_info * Map, int area)
{
    int i, nholes, isle;
    GEOSGeometry *boundary, **holes;
    
    G_debug(3, "Vect_read_area_geos(): area = %d", area);

    boundary = GEOSGeom_createLinearRing(Vect_get_area_points_geos(Map, area));
    if (!boundary) {
        G_fatal_error(_("Vect_read_area_geos(): unable to read area id %d"),
                      area);
    }

    nholes = Vect_get_area_num_isles(Map, area);
    holes = (GEOSGeometry **) G_malloc(nholes * sizeof(GEOSGeometry *));
    for (i = 0; i < nholes; i++) {
        isle = Vect_get_area_isle(Map, area, i);
        if (isle < 1) {
            nholes--;
            continue;
        }
        holes[i] = GEOSGeom_createLinearRing(Vect_get_isle_points_geos(Map, isle));
        if (!(holes[i]))
            G_fatal_error(_("Vect_read_area_geos(): unable to read isle id %d of area id %d"),
                          isle, area);
    }
    
    return GEOSGeom_createPolygon(boundary, holes, nholes);
}

GEOSGeometry *Vect_line_to_geos(struct Map_info *Map,
                                const struct line_pnts *points, int type)
{
    int i, with_z;
    GEOSGeometry *geom;
    GEOSCoordSequence *pseq;

    G_debug(3, "Vect_line_to_geos(): type = %d", type);
    
    with_z = Vect_is_3d(Map);
    
    /* read only points / lines / boundaries */
    if (!(type & (GV_POINT | GV_LINES)))
        return NULL;

    if (type == GV_POINT) { 
        if (points->n_points != 1)
            /* point is not valid */
            return NULL;
    }
    else {                      
        if (points->n_points < 2)
            /* line/boundary is not valid */
            return NULL;
    }
    
    pseq = GEOSCoordSeq_create(points->n_points, with_z ? 3 : 2);
    
    for (i = 0; i < points->n_points; i++) {
        GEOSCoordSeq_setX(pseq, i, points->x[i]);
        GEOSCoordSeq_setY(pseq, i, points->y[i]);
        if (with_z)
            GEOSCoordSeq_setZ(pseq, i, points->z[i]);
    }

    if (type == GV_POINT)
        geom = GEOSGeom_createPoint(pseq);
    else if (type == GV_LINE)
        geom = GEOSGeom_createLineString(pseq);
    else { /* boundary */
        geom = GEOSGeom_createLineString(pseq);
        if (GEOSisRing(geom)) {
            /* GEOSGeom_destroy(geom); */
            geom = GEOSGeom_createLinearRing(pseq);
        }
    }
    
    /* GEOSCoordSeq_destroy(pseq); */

    return geom;
}

GEOSGeometry *Vect__read_line_geos(struct Map_info *Map, long offset, int *type)
{
    int ftype;
    
    GEOSGeometry *geom;
    GEOSCoordSequence *pseq;
    
    pseq = V1_read_line_geos(Map, offset, &ftype);
    if (!pseq)
        G_fatal_error(_("Unable to read line offset %ld"), offset);
    
    if (ftype & GV_POINT) {
        G_debug(3, "    geos_type = point");
        geom = GEOSGeom_createPoint(pseq);
    }
    else if (ftype & GV_LINE) {
        G_debug(3, "    geos_type = linestring");
        geom = GEOSGeom_createLineString(pseq);
    }
    else { /* boundary */
        geom = GEOSGeom_createLineString(pseq);
        if (GEOSisRing(geom)) {
            /* GEOSGeom_destroy(geom); */
            geom = GEOSGeom_createLinearRing(pseq);
            G_debug(3, "    geos_type = linearring");
        }
        else {
            G_debug(3, "    geos_type = linestring");
        }
    }
        
    /* GEOSCoordSeq_destroy(pseq); */
    
    if (type)
      *type = ftype;
    
    return geom;
}

GEOSCoordSequence *V2_read_line_geos(struct Map_info *Map, int line)
{
    int ftype;
    P_LINE *Line;
    
    G_debug(3, "V2_read_line_geos(): line = %d", line);
    
    Line = Map->plus.Line[line];

    if (Line == NULL)
        G_fatal_error("V2_read_line_geos(): %s %d",
                      _("Attempt to read dead line"), line);
    
    return V1_read_line_geos(Map, Line->offset, &ftype);
}


GEOSCoordSequence *V1_read_line_geos(struct Map_info *Map, long offset, int *type)
{
    int i, n_points;
    int do_cats, n_cats;
    char rhead, nc;
    long size;
    double *x, *y, *z;
    
    GEOSCoordSequence *pseq;
    
    G_debug(3, "V1_read_line_geos(): offset = %ld", offset);
    
    Map->head.last_offset = offset;
    
    /* reads must set in_head, but writes use default */
    dig_set_cur_port(&(Map->head.port));
    
    dig_fseek(&(Map->dig_fp), offset, 0);
    
    if (0 >= dig__fread_port_C(&rhead, 1, &(Map->dig_fp)))
        return NULL;            /* end of file */
    
    if (!(rhead & 0x01))        /* dead line */
        return GEOSCoordSeq_create(0, (Map->head.with_z) ? 3 : 2);

    if (rhead & 0x02)           /* categories exists */
        do_cats = 1;            /* do not return here let file offset moves forward to next */
    else                        /* line */
        do_cats = 0;
    
    rhead >>= 2;
    *type = dig_type_from_store((int) rhead);
    
    /* read only points / lines / boundaries */
    if (!(*type & (GV_POINT | GV_LINES)))
        return GEOSCoordSeq_create(0, (Map->head.with_z) ? 3 : 2);
 
    /* skip categories */
    if (do_cats) {
        if (Map->head.Version_Minor == 1) {     /* coor format 5.1 */
            if (0 >= dig__fread_port_I(&n_cats, 1, &(Map->dig_fp)))
                return NULL;
        }
        else {                                  /* coor format 5.0 */
            if (0 >= dig__fread_port_C(&nc, 1, &(Map->dig_fp)))
                return NULL;
            n_cats = (int) nc;
        }
        G_debug(3, "    n_cats = %d", n_cats);

        if (Map->head.Version_Minor == 1) {     /* coor format 5.1 */
            size = (2 * PORT_INT) * n_cats;
        }
        else {                          /* coor format 5.0 */
            size = (PORT_SHORT + PORT_INT) * n_cats;
        }
        dig_fseek(&(Map->dig_fp), size, SEEK_CUR);
    }

    if (*type & GV_POINTS) {
            n_points = 1;
    }
    else {
        if (0 >= dig__fread_port_I(&n_points, 1, &(Map->dig_fp)))
            return NULL;
    }
    
    G_debug(3, "    n_points = %d dim = %d", n_points, (Map->head.with_z) ? 3 : 2);
    
    pseq = GEOSCoordSeq_create(n_points, (Map->head.with_z) ? 3 : 2);
    
    x = (double *) G_malloc(n_points * sizeof(double));
    y = (double *) G_malloc(n_points * sizeof(double));
    if (Map->head.with_z)
        z = (double *) G_malloc(n_points * sizeof(double));
    else
        z = NULL;
    
    if (0 >= dig__fread_port_D(x, n_points, &(Map->dig_fp)))
        return NULL; /* end of file */

    if (0 >= dig__fread_port_D(y, n_points, &(Map->dig_fp)))
        return NULL; /* end of file */

    if (Map->head.with_z) {
        if (0 >= dig__fread_port_D(z, n_points, &(Map->dig_fp)))
            return NULL; /* end of file */

    }

    for (i = 0; i < n_points; i++) {
        GEOSCoordSeq_setX(pseq, i, x[i]);
        GEOSCoordSeq_setY(pseq, i, y[i]);
        if (Map->head.with_z)
            GEOSCoordSeq_setZ(pseq, i, z[i]);
    }
    
    G_debug(3, "    off = %ld", dig_ftell(&(Map->dig_fp)));
    
    G_free((void *) x);
    G_free((void *) y);
    if (z)
        G_free((void *) z);
    
    return pseq;
}

GEOSCoordSequence *Vect_get_area_points_geos(struct Map_info *Map, int area)
{
    struct Plus_head *Plus;
    P_AREA *Area;
    
    G_debug(3, "Vect_get_area_points_geos(): area = %d", area);
    
    Plus = &(Map->plus);
    Area = Plus->Area[area];

    if (Area == NULL) {         /* dead area */
        G_warning(_("Attempt to read points of nonexistent area id %d"), area);
        return NULL;            /* error , because we should not read dead areas */
    }
    
    return read_polygon_points(Map, Area->n_lines, Area->lines);
}

GEOSCoordSequence *Vect_get_isle_points_geos(struct Map_info *Map, int isle)
{
    struct Plus_head *Plus;
    P_ISLE *Isle;
    
    G_debug(3, "Vect_get_isle_points_geos(): isle = %d", isle);

    Plus = &(Map->plus);
    Isle = Plus->Isle[isle];

    return read_polygon_points(Map, Isle->n_lines, Isle->lines);
}

GEOSCoordSequence *read_polygon_points(struct Map_info *Map, int n_lines, int *lines)
{
    int i, j, k;
    int line, aline;
    unsigned int n_points, n_points_shell;
    double x, y, z;
    int *dir;
    
    GEOSCoordSequence **pseq, *pseq_shell;

    G_debug(3, "  n_lines = %d", n_lines);
    pseq = (GEOSCoordSequence **) G_malloc(n_lines * sizeof(GEOSCoordSequence *));
    dir  = (int*) G_malloc(n_lines * sizeof(int));

    n_points_shell = 0;
    for (i = 0; i < n_lines; i++) {
        line = lines[i];
        aline = abs(line);
        G_debug(3, "  append line(%d) = %d", i, line);

        if (line > 0)
            dir[i] = GV_FORWARD;
        else
            dir[i] = GV_BACKWARD;
        
        pseq[i] = V2_read_line_geos(Map, aline);
        if (!(pseq[i])) {
            G_fatal_error(_("Unable to read feature id %d"), aline);
        }
        
        GEOSCoordSeq_getSize(pseq[i], &n_points);
        G_debug(3, "  line n_points = %d", n_points);
        n_points_shell += n_points;
    }

    /* create shell (outer ring) */
    pseq_shell = GEOSCoordSeq_create(n_points_shell, Map->head.with_z ? 3 : 2);
    k = 0;
    for (i = 0; i < n_lines; i++) {
        GEOSCoordSeq_getSize(pseq[i], &n_points);
        if (dir[i] == GV_FORWARD) {
            for (j = 0; j < (int) n_points; j++, k++) {
                GEOSCoordSeq_getX(pseq[i], j, &x);
                GEOSCoordSeq_setX(pseq_shell, k, x);
                
                GEOSCoordSeq_getY(pseq[i], j, &y);
                GEOSCoordSeq_setY(pseq_shell, k, y);
                
                if (Map->head.with_z) {
                    GEOSCoordSeq_getY(pseq[i], j, &z);
                    GEOSCoordSeq_setZ(pseq_shell, k, z);
                }
            }
        }
        else { /* GV_BACKWARD */
            for (j = (int) n_points - 1; j > -1; j--, k++) {
                GEOSCoordSeq_getX(pseq[i], j, &x);
                GEOSCoordSeq_setX(pseq_shell, k, x);
                
                GEOSCoordSeq_getY(pseq[i], j, &y);
                GEOSCoordSeq_setY(pseq_shell, k, y);
                
                if (Map->head.with_z) {
                    GEOSCoordSeq_getY(pseq[i], j, &z);
                    GEOSCoordSeq_setZ(pseq_shell, k, z);
                }
            }
        }
    }
    
    G_free((void *) pseq);
    G_free((void *) dir);
    
    return pseq_shell;
}
#endif /* HAVE_GEOS */