lwt_AddPolygon()

LWT_ELEMID* lwt_AddPolygon	(	LWT_TOPOLOGY *	topo,
		LWPOLY *	poly,
		double	tol,
		int *	nfaces
	)

Adds a polygon to the topology.

The boundary of the given polygon will snap to existing nodes or edges within given tolerance. Existing edges or faces may be split by the boundary of the polygon.

Parameters

topo	the topology to operate on
poly	the polygon to add
tol	snap tolerance, the topology tolerance will be used if 0
nfaces	output parameter, will be set to number of faces the polygon was split into, or -1 on error (liblwgeom error handler will be invoked with error message)

Returns

an array of <nfaces> face identifiers that sewed togheter will build up the input polygon (after snapping). Caller will need to free the array using lwfree(), if not null.

Definition at line 5875 of file lwgeom_topo.c.

References _LWT_MINTOLERANCE, _lwt_release_faces(), LWT_TOPOLOGY_T::be_iface, covers(), LWT_ISO_FACE::face_id, gbox_expand(), lwalloc(), LWDEBUGF, lwerror(), lwfree(), LWGEOM2GEOS(), lwgeom_free(), lwgeom_geos_errmsg, lwgeom_get_bbox(), lwline_construct(), lwline_free(), lwpoly_as_lwgeom(), lwt_AddLine(), lwt_be_getFaceWithinBox2D(), lwt_be_lastErrorMessage(), LWT_COL_FACE_ALL, lwt_GetFaceGeometry(), LWTFMT_ELEMID, ptarray_clone(), and LWT_TOPOLOGY_T::srid.

{
  int i;
  *nfaces = -1; /* error condition, by default */
  int num;
  LWT_ISO_FACE *faces;
  int nfacesinbox;
  LWT_ELEMID *ids = NULL;
  GBOX qbox;
  const GEOSPreparedGeometry *ppoly;
  GEOSGeometry *polyg;

  /* Get tolerance, if 0 was given */
  if ( ! tol ) tol = _LWT_MINTOLERANCE( topo, (LWGEOM*)poly );
  LWDEBUGF(1, "Working tolerance:%.15g", tol);

  /* Add each ring as an edge */
  for ( i=0; i<poly->nrings; ++i )
  {
    LWLINE *line;
    POINTARRAY *pa;
    LWT_ELEMID *eids;
    int nedges;

    pa = ptarray_clone(poly->rings[i]);
    line = lwline_construct(topo->srid, NULL, pa);
    eids = lwt_AddLine( topo, line, tol, &nedges );
    if ( nedges < 0 ) {
      /* probably too late as lwt_AddLine invoked lwerror */
      lwline_free(line);
      lwerror("Error adding ring %d of polygon", i);
      return NULL;
    }
    lwline_free(line);
    lwfree(eids);
  }

  /*
  -- Find faces covered by input polygon
  -- NOTE: potential snapping changed polygon edges
  */
  qbox = *lwgeom_get_bbox( lwpoly_as_lwgeom(poly) );
  gbox_expand(&qbox, tol);
  faces = lwt_be_getFaceWithinBox2D( topo, &qbox, &nfacesinbox,
                                     LWT_COL_FACE_ALL, 0 );
  if ( nfacesinbox == -1 )
  {
    lwfree(ids);
    lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
    return NULL;
  }

  num = 0;
  if ( nfacesinbox )
  {
    polyg = LWGEOM2GEOS(lwpoly_as_lwgeom(poly), 0);
    if ( ! polyg )
    {
      _lwt_release_faces(faces, nfacesinbox);
      lwerror("Could not convert poly geometry to GEOS: %s", lwgeom_geos_errmsg);
      return NULL;
    }
    ppoly = GEOSPrepare(polyg);
    ids = lwalloc(sizeof(LWT_ELEMID)*nfacesinbox);
    for ( i=0; i<nfacesinbox; ++i )
    {
      LWT_ISO_FACE *f = &(faces[i]);
      LWGEOM *fg;
      GEOSGeometry *fgg, *sp;
      int covers;

      /* check if a point on this face surface is covered by our polygon */
      fg = lwt_GetFaceGeometry( topo, f->face_id );
      if ( ! fg )
      {
        i = f->face_id; /* so we can destroy faces */
        GEOSPreparedGeom_destroy(ppoly);
        GEOSGeom_destroy(polyg);
        lwfree(ids);
        _lwt_release_faces(faces, nfacesinbox);
        lwerror("Could not get geometry of face %" LWTFMT_ELEMID, i);
        return NULL;
      }
      /* check if a point on this face's surface is covered by our polygon */
      fgg = LWGEOM2GEOS(fg, 0);
      lwgeom_free(fg);
      if ( ! fgg )
      {
        GEOSPreparedGeom_destroy(ppoly);
        GEOSGeom_destroy(polyg);
        _lwt_release_faces(faces, nfacesinbox);
        lwerror("Could not convert edge geometry to GEOS: %s", lwgeom_geos_errmsg);
        return NULL;
      }
      sp = GEOSPointOnSurface(fgg);
      GEOSGeom_destroy(fgg);
      if ( ! sp )
      {
        GEOSPreparedGeom_destroy(ppoly);
        GEOSGeom_destroy(polyg);
        _lwt_release_faces(faces, nfacesinbox);
        lwerror("Could not find point on face surface: %s", lwgeom_geos_errmsg);
        return NULL;
      }
      covers = GEOSPreparedCovers( ppoly, sp );
      GEOSGeom_destroy(sp);
      if (covers == 2)
      {
        GEOSPreparedGeom_destroy(ppoly);
        GEOSGeom_destroy(polyg);
        _lwt_release_faces(faces, nfacesinbox);
        lwerror("PreparedCovers error: %s", lwgeom_geos_errmsg);
        return NULL;
      }
      if ( ! covers )
      {
        continue; /* we're not composed by this face */
      }

      /* TODO: avoid duplicates ? */
      ids[num++] = f->face_id;
    }
    GEOSPreparedGeom_destroy(ppoly);
    GEOSGeom_destroy(polyg);
    _lwt_release_faces(faces, nfacesinbox);
  }

  /* possibly 0 if non face's surface point was found
   * to be covered by input polygon */
  *nfaces = num;

  return ids;
}

Here is the call graph for this function: