lwgeom_in_gml.c

Go to the documentation of this file.

/**********************************************************************
 *
 * PostGIS - Spatial Types for PostgreSQL
 * http://postgis.net
 *
 * PostGIS is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 2 of the License, or
 * (at your option) any later version.
 *
 * PostGIS is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with PostGIS.  If not, see <http://www.gnu.org/licenses/>.
 *
 **********************************************************************
 *
 * Copyright 2009 - 2010 Oslandia
 *
 **********************************************************************/
 
 
#include <libxml/tree.h>
#include <libxml/parser.h>
#include <libxml/xpath.h>
#include <libxml/xpathInternals.h>
 
#include "postgres.h"
#include "executor/spi.h"
#include "utils/builtins.h"
 
#include "../postgis_config.h"
#include "lwgeom_pg.h"
#include "liblwgeom.h"
#include "lwgeom_transform.h"
 
 
Datum geom_from_gml(PG_FUNCTION_ARGS);
static LWGEOM *lwgeom_from_gml(const char *wkt, int xml_size);
static LWGEOM* parse_gml(xmlNodePtr xnode, bool *hasz, int *root_srid);
 
typedef struct struct_gmlSrs
{
        int32_t srid;
        bool reverse_axis;
}
gmlSrs;
 
#define XLINK_NS        ((char *) "http://www.w3.org/1999/xlink")
#define GML_NS          ((char *) "http://www.opengis.net/gml")
#define GML32_NS        ((char *) "http://www.opengis.net/gml/3.2")
 
 
 
static void gml_lwpgerror(char *msg, __attribute__((__unused__)) int error_code)
{
        POSTGIS_DEBUGF(3, "ST_GeomFromGML ERROR %i", error_code);
        lwpgerror("%s", msg);
}
 
PG_FUNCTION_INFO_V1(geom_from_gml);
Datum geom_from_gml(PG_FUNCTION_ARGS)
{
        GSERIALIZED *geom;
        text *xml_input;
        LWGEOM *lwgeom;
        char *xml;
        int root_srid=SRID_UNKNOWN;
        int xml_size;
 
        /* Get the GML stream */
        if (PG_ARGISNULL(0)) PG_RETURN_NULL();
        xml_input = PG_GETARG_TEXT_P(0);
        xml = text_to_cstring(xml_input);
        xml_size = VARSIZE_ANY_EXHDR(xml_input);
 
        /* Zero for undefined */
        root_srid = PG_GETARG_INT32(1);
 
#if POSTGIS_PROJ_VERSION < 60
        /* Internally lwgeom_from_gml calls gml_reproject_pa which, for PROJ before 6, called GetProj4String.
         * That function requires access to spatial_ref_sys, so in order to have it ready we need to ensure
         * the internal cache is initialized
         */
        postgis_initialize_cache(fcinfo);
#endif
        lwgeom = lwgeom_from_gml(xml, xml_size);
        if ( root_srid != SRID_UNKNOWN )
                lwgeom->srid = root_srid;
 
        geom = geometry_serialize(lwgeom);
        lwgeom_free(lwgeom);
 
        PG_RETURN_POINTER(geom);
}
 
 
static inline bool
is_gml_element(xmlNodePtr xn, const char *gml_name)
{
        char *colon_pos;
        char *node_name;
 
        /* Not an element node, can't do anything */
        if (!xn || xn->type != XML_ELEMENT_NODE)
                return false;
 
        /* If there's a colon in the element name, */
        /* move past it before checking for equality with */
        /* the element name we are looking for */
        node_name = (char*)xn->name;
        colon_pos = strchr(node_name, ':');
        if (colon_pos)
                node_name = colon_pos + 1;
 
        return strcmp(node_name, gml_name) == 0;
}
 
 
static bool is_gml_namespace(xmlNodePtr xnode, bool is_strict)
{
        xmlNsPtr *ns, *p;
 
        ns = xmlGetNsList(xnode->doc, xnode);
        /*
         * If no namespace is available we could return true anyway
         * (because we work only on GML fragment, we don't want to
         *  'oblige' to add namespace on the geometry root node)
         */
        if (ns == NULL) { return !is_strict; }
 
        /*
         * Handle namespaces:
         *  - http://www.opengis.net/gml      (GML 3.1.1 and priors)
         *  - http://www.opengis.net/gml/3.2  (GML 3.2.1)
         */
        for (p=ns ; *p ; p++)
        {
                if ((*p)->href == NULL || (*p)->prefix == NULL ||
                     xnode->ns == NULL || xnode->ns->prefix == NULL) continue;
 
                if (!xmlStrcmp(xnode->ns->prefix, (*p)->prefix))
                {
                        if (    !strcmp((char *) (*p)->href, GML_NS)
                             || !strcmp((char *) (*p)->href, GML32_NS))
                        {
                                xmlFree(ns);
                                return true;
                        } else {
                                xmlFree(ns);
                                return false;
                        }
                }
        }
 
        xmlFree(ns);
        return !is_strict; /* Same reason here to not return false */
}
 
 
static xmlChar *gmlGetProp(xmlNodePtr xnode, const char *charProp)
{
        xmlChar *value;
        xmlChar *prop = (xmlChar*)charProp;
 
        if (!is_gml_namespace(xnode, true))
                return xmlGetProp(xnode, prop);
        /*
         * Handle namespaces:
         *  - http://www.opengis.net/gml      (GML 3.1.1 and priors)
         *  - http://www.opengis.net/gml/3.2  (GML 3.2.1)
         */
        value = xmlGetNsProp(xnode, prop, (xmlChar *) GML_NS);
        if (value == NULL) value = xmlGetNsProp(xnode, prop, (xmlChar *) GML32_NS);
 
        /* In last case try without explicit namespace */
        if (value == NULL) value = xmlGetNoNsProp(xnode, prop);
 
        return value;
}
 
 
static bool is_xlink(xmlNodePtr node)
{
        xmlChar *prop;
 
        prop = xmlGetNsProp(node, (xmlChar *)"type", (xmlChar *) XLINK_NS);
        if (prop == NULL) return false;
        if (strcmp((char *) prop, "simple"))
        {
                xmlFree(prop);
                return false;
        }
 
        prop = xmlGetNsProp(node, (xmlChar *)"href", (xmlChar *) XLINK_NS);
        if (prop == NULL) return false;
        if (prop[0] != '#')
        {
                xmlFree(prop);
                return false;
        }
        xmlFree(prop);
 
        return true;
}
 
 
static xmlNodePtr get_xlink_node(xmlNodePtr xnode)
{
        char *id;
        xmlNsPtr *ns, *n;
        xmlXPathContext *ctx;
        xmlXPathObject *xpath;
        xmlNodePtr node, ret_node;
        xmlChar *href, *p, *node_id;
 
        href = xmlGetNsProp(xnode, (xmlChar *)"href", (xmlChar *) XLINK_NS);
        id = lwalloc((xmlStrlen(xnode->ns->prefix) * 2 + xmlStrlen(xnode->name)
                      + xmlStrlen(href) + sizeof("//:[@:id='']") + 1));
        p = href;
        p++; /* ignore '#' first char */
 
        /* XPath pattern look like: //gml:point[@gml:id='p1'] */
        sprintf(id, "//%s:%s[@%s:id='%s']",     (char *) xnode->ns->prefix,
                (char *) xnode->name,
                (char *) xnode->ns->prefix,
                (char *) p);
 
        ctx = xmlXPathNewContext(xnode->doc);
        if (ctx == NULL)
        {
                xmlFree(href);
                lwfree(id);
                return NULL;
        }
 
        /* Handle namespaces */
        ns = xmlGetNsList(xnode->doc, xnode);
        for (n=ns ; *n; n++) xmlXPathRegisterNs(ctx, (*n)->prefix, (*n)->href);
        xmlFree(ns);
 
        /* Execute XPath expression */
        xpath = xmlXPathEvalExpression((xmlChar *) id, ctx);
        lwfree(id);
        if (xpath == NULL || xpath->nodesetval == NULL || xpath->nodesetval->nodeNr != 1)
        {
                xmlFree(href);
                xmlXPathFreeObject(xpath);
                xmlXPathFreeContext(ctx);
                return NULL;
        }
        ret_node = xpath->nodesetval->nodeTab[0];
        xmlXPathFreeObject(xpath);
        xmlXPathFreeContext(ctx);
 
        /* Protection against circular calls */
        for (node = xnode ; node != NULL ; node = node->parent)
        {
                if (node->type != XML_ELEMENT_NODE) continue;
                node_id = gmlGetProp(node, "id");
                if (node_id != NULL)
                {
                        if (!xmlStrcmp(node_id, p))
                                gml_lwpgerror("invalid GML representation", 2);
                        xmlFree(node_id);
                }
        }
 
        xmlFree(href);
        return ret_node;
}
 
 
 
#if POSTGIS_PROJ_VERSION < 60
 
static POINTARRAY *
gml_reproject_pa(POINTARRAY *pa, int32_t srid_in, int32_t srid_out)
{
        PJ pj;
        char *text_in, *text_out;
 
        if (srid_in == SRID_UNKNOWN) return pa; /* nothing to do */
        if (srid_out == SRID_UNKNOWN) gml_lwpgerror("invalid GML representation", 3);
 
        text_in = GetProj4String(srid_in);
        text_out = GetProj4String(srid_out);
 
        pj.pj_from = projpj_from_string(text_in);
        pj.pj_to = projpj_from_string(text_out);
 
        lwfree(text_in);
        lwfree(text_out);
 
        if ( ptarray_transform(pa, &pj) == LW_FAILURE )
        {
                elog(ERROR, "gml_reproject_pa: reprojection failed");
        }
 
        pj_free(pj.pj_from);
        pj_free(pj.pj_to);
 
        return pa;
}
#else
/*
 * TODO: rework GML projection handling to skip the spatial_ref_sys
 * lookups, and use the Proj 6+ EPSG catalogue and built-in SRID
 * lookups directly. Drop this ugly hack.
 */
static POINTARRAY *
gml_reproject_pa(POINTARRAY *pa, int32_t epsg_in, int32_t epsg_out)
{
        PJ *pj;
        LWPROJ *lwp;
        char text_in[16];
        char text_out[16];
 
        if (epsg_in == SRID_UNKNOWN)
                return pa; /* nothing to do */
 
        if (epsg_out == SRID_UNKNOWN)
        {
                gml_lwpgerror("invalid GML representation", 3);
                return NULL;
        }
 
        snprintf(text_in, 16, "EPSG:%d", epsg_in);
        snprintf(text_out, 16, "EPSG:%d", epsg_out);
        pj = proj_create_crs_to_crs(NULL, text_in, text_out, NULL);
 
        lwp = lwproj_from_PJ(pj, LW_FALSE);
        if (!lwp)
        {
                proj_destroy(pj);
                gml_lwpgerror("Could not create LWPROJ*", 57);
                return NULL;
        }
 
        if (ptarray_transform(pa, lwp) == LW_FAILURE)
        {
                proj_destroy(pj);
                elog(ERROR, "gml_reproject_pa: reprojection failed");
                return NULL;
        }
        proj_destroy(pj);
        pfree(lwp);
 
        return pa;
}
#endif /* POSTGIS_PROJ_VERSION */
 
 
static int
gml_is_srid_planar(int32_t srid)
{
        char *result;
        char query[256];
        int is_planar, err;
 
        if (SPI_OK_CONNECT != SPI_connect ())
                lwpgerror("gml_is_srid_planar: could not connect to SPI manager");
 
        /* A way to find if this projection is planar or geocentric */
        sprintf(query, "SELECT position('+units=m ' in proj4text) \
                        FROM spatial_ref_sys WHERE srid='%d'", srid);
 
        err = SPI_exec(query, 1);
        if (err < 0) lwpgerror("gml_is_srid_planar: error executing query %d", err);
 
        /* No entry in spatial_ref_sys */
        if (SPI_processed <= 0)
        {
                SPI_finish();
                return -1;
        }
 
        result = SPI_getvalue(SPI_tuptable->vals[0], SPI_tuptable->tupdesc, 1);
        is_planar = atoi(result);
        SPI_finish();
 
        return is_planar;
}
 
 
static void parse_gml_srs(xmlNodePtr xnode, gmlSrs *srs)
{
        char *p;
        int is_planar;
        xmlNodePtr node;
        xmlChar *srsname;
        bool latlon = false;
        char sep = ':';
 
        node = xnode;
        srsname = gmlGetProp(node, "srsName");
        /*printf("srsname %s\n",srsname);*/
        if (!srsname)
        {
                if (node->parent == NULL)
                {
                        srs->srid = SRID_UNKNOWN;
                        srs->reverse_axis = false;
                        return;
                }
                parse_gml_srs(node->parent, srs);
        }
        else
        {
                /* Severals     srsName formats are available...
                 *  cf WFS 1.1.0 -> 9.2 (p36)
                 *  cf ISO 19142:2009 -> 7.9.2.4.4 (p34)
                 *  cf RFC 5165 <http://tools.ietf.org/html/rfc5165>
                 *  cf CITE WFS-1.1 (GetFeature-tc17.2)
                 */
 
                /* SRS pattern like:    EPSG:4326
                                        urn:EPSG:geographicCRS:4326
                                        urn:ogc:def:crs:EPSG:4326
                                        urn:ogc:def:crs:EPSG::4326
                                        urn:ogc:def:crs:EPSG:6.6:4326
                                        urn:x-ogc:def:crs:EPSG:6.6:4326
                                        http://www.opengis.net/gml/srs/epsg.xml#4326
                                        http://www.epsg.org/6.11.2/4326
                */
 
                if (!strncmp((char *) srsname, "EPSG:", 5))
                {
                        sep = ':';
                        latlon = false;
                }
                else if (!strncmp((char *) srsname, "urn:ogc:def:crs:EPSG:", 21)
                         || !strncmp((char *) srsname, "urn:x-ogc:def:crs:EPSG:", 23)
                         || !strncmp((char *) srsname, "urn:EPSG:geographicCRS:", 23))
                {
                        sep = ':';
                        latlon = true;
                }
                else if (!strncmp((char *) srsname,
                                  "http://www.opengis.net/gml/srs/epsg.xml#", 40))
                {
                        sep = '#';
                        latlon = false;
                }
                else gml_lwpgerror("unknown spatial reference system", 4);
 
                /* retrieve the last ':' or '#' char */
                for (p = (char *) srsname ; *p ; p++);
                for (--p ; *p != sep ; p--)
                        if (!isdigit(*p)) gml_lwpgerror("unknown spatial reference system", 5);
 
                srs->srid = atoi(++p);
 
                /* Check into spatial_ref_sys that this SRID really exist */
                is_planar = gml_is_srid_planar(srs->srid);
                if (srs->srid == SRID_UNKNOWN || is_planar == -1)
                        gml_lwpgerror("unknown spatial reference system", 6);
 
                /* About lat/lon issue, Cf: http://tinyurl.com/yjpr55z */
                srs->reverse_axis = !is_planar && latlon;
 
                xmlFree(srsname);
                return;
        }
}
 
 
static double parse_gml_double(char *d, bool space_before, bool space_after)
{
        char *p;
        int st;
        enum states
        {
                INIT            = 0,
                NEED_DIG        = 1,
                DIG             = 2,
                NEED_DIG_DEC    = 3,
                DIG_DEC         = 4,
                EXP             = 5,
                NEED_DIG_EXP    = 6,
                DIG_EXP         = 7,
                END             = 8
        };
 
        /*
         * Double pattern
         * [-|\+]?[0-9]+(\.)?([0-9]+)?([Ee](\+|-)?[0-9]+)?
         * We could also meet spaces before and/or after
         * this pattern upon parameters
         */
 
        if (space_before) while (isspace(*d)) d++;
        for (st = INIT, p = d ; *p ; p++)
        {
 
                if (isdigit(*p))
                {
                        if (st == INIT || st == NEED_DIG)               st = DIG;
                        else if (st == NEED_DIG_DEC)                    st = DIG_DEC;
                        else if (st == NEED_DIG_EXP || st == EXP)       st = DIG_EXP;
                        else if (st == DIG || st == DIG_DEC || st == DIG_EXP);
                        else gml_lwpgerror("invalid GML representation", 7);
                }
                else if (*p == '.')
                {
                        if      (st == DIG)                             st = NEED_DIG_DEC;
                        else    gml_lwpgerror("invalid GML representation", 8);
                }
                else if (*p == '-' || *p == '+')
                {
                        if      (st == INIT)                            st = NEED_DIG;
                        else if (st == EXP)                             st = NEED_DIG_EXP;
                        else    gml_lwpgerror("invalid GML representation", 9);
                }
                else if (*p == 'e' || *p == 'E')
                {
                        if      (st == DIG || st == DIG_DEC)            st = EXP;
                        else    gml_lwpgerror("invalid GML representation", 10);
                }
                else if (isspace(*p))
                {
                        if (!space_after) gml_lwpgerror("invalid GML representation", 11);
                        if (st == DIG || st == DIG_DEC || st == DIG_EXP)st = END;
                        else if (st == NEED_DIG_DEC)                    st = END;
                        else if (st == END);
                        else    gml_lwpgerror("invalid GML representation", 12);
                }
                else  gml_lwpgerror("invalid GML representation", 13);
        }
 
        if (st != DIG && st != NEED_DIG_DEC && st != DIG_DEC && st != DIG_EXP && st != END)
                gml_lwpgerror("invalid GML representation", 14);
 
        return atof(d);
}
 
 
static POINTARRAY* parse_gml_coordinates(xmlNodePtr xnode, bool *hasz)
{
        xmlChar *gml_coord, *gml_ts, *gml_cs, *gml_dec;
        char cs, ts, dec;
        POINTARRAY *dpa;
        int gml_dims;
        char *p, *q;
        bool digit;
        POINT4D pt;
 
        /* We begin to retrieve coordinates string */
        gml_coord = xmlNodeGetContent(xnode);
        p = (char *) gml_coord;
 
        /* Default GML coordinates pattern:     x1,y1 x2,y2
         *                                      x1,y1,z1 x2,y2,z2
         *
         * Cf GML 2.1.2 -> 4.3.1 (p18)
         */
 
        /* Retrieve separator between coordinates tuples */
        gml_ts = gmlGetProp(xnode, "ts");
        if (gml_ts == NULL) ts = ' ';
        else
        {
                if (xmlStrlen(gml_ts) > 1 || isdigit(gml_ts[0]))
                        gml_lwpgerror("invalid GML representation", 15);
                ts = gml_ts[0];
                xmlFree(gml_ts);
        }
 
        /* Retrieve separator between each coordinate */
        gml_cs = gmlGetProp(xnode, "cs");
        if (gml_cs == NULL) cs = ',';
        else
        {
                if (xmlStrlen(gml_cs) > 1 || isdigit(gml_cs[0]))
                        gml_lwpgerror("invalid GML representation", 16);
                cs = gml_cs[0];
                xmlFree(gml_cs);
        }
 
        /* Retrieve decimal separator */
        gml_dec = gmlGetProp(xnode, "decimal");
        if (gml_dec == NULL) dec = '.';
        else
        {
                if (xmlStrlen(gml_dec) > 1 || isdigit(gml_dec[0]))
                        gml_lwpgerror("invalid GML representation", 17);
                dec = gml_dec[0];
                xmlFree(gml_dec);
        }
 
        if (cs == ts || cs == dec || ts == dec)
                gml_lwpgerror("invalid GML representation", 18);
 
        /* HasZ, !HasM, 1 Point */
        dpa = ptarray_construct_empty(1, 0, 1);
 
        while (isspace(*p)) p++;                /* Eat extra whitespaces if any */
        for (q = p, gml_dims=0, digit = false ; *p ; p++)
        {
 
                if (isdigit(*p)) digit = true;  /* One state parser */
 
                /* Coordinate Separator */
                if (*p == cs)
                {
                        *p = '\0';
                        gml_dims++;
 
                        if (*(p+1) == '\0') gml_lwpgerror("invalid GML representation", 19);
 
                        if      (gml_dims == 1) pt.x = parse_gml_double(q, false, true);
                        else if (gml_dims == 2) pt.y = parse_gml_double(q, false, true);
 
                        q = p+1;
 
                        /* Tuple Separator (or end string) */
                }
                else if (digit && (*p == ts || *(p+1) == '\0'))
                {
                        if (*p == ts) *p = '\0';
                        gml_dims++;
 
                        if (gml_dims < 2 || gml_dims > 3)
                                gml_lwpgerror("invalid GML representation", 20);
 
                        if (gml_dims == 3)
                                pt.z = parse_gml_double(q, false, true);
                        else
                        {
                                pt.y = parse_gml_double(q, false, true);
                                *hasz = false;
                        }
 
                        ptarray_append_point(dpa, &pt, LW_TRUE);
                        digit = false;
 
                        q = p+1;
                        gml_dims = 0;
 
                        /* Need to put standard decimal separator to atof handle */
                }
                else if (*p == dec && dec != '.') *p = '.';
        }
 
        xmlFree(gml_coord);
 
        return dpa; /* ptarray_clone_deep(dpa); */
}
 
 
static POINTARRAY* parse_gml_coord(xmlNodePtr xnode, bool *hasz)
{
        xmlNodePtr xyz;
        POINTARRAY *dpa;
        bool x,y,z;
        xmlChar *c;
        POINT4D p;
 
        /* HasZ?, !HasM, 1 Point */
        dpa = ptarray_construct_empty(1, 0, 1);
 
        x = y = z = false;
        for (xyz = xnode->children ; xyz != NULL ; xyz = xyz->next)
        {
                if (xyz->type != XML_ELEMENT_NODE) continue;
                if (!is_gml_namespace(xyz, false)) continue;
 
                if (is_gml_element(xyz, "X"))
                {
                        if (x) gml_lwpgerror("invalid GML representation", 21);
                        c = xmlNodeGetContent(xyz);
                        p.x = parse_gml_double((char *) c, true, true);
                        x = true;
                        xmlFree(c);
                }
                else  if (is_gml_element(xyz, "Y"))
                {
                        if (y) gml_lwpgerror("invalid GML representation", 22);
                        c = xmlNodeGetContent(xyz);
                        p.y = parse_gml_double((char *) c, true, true);
                        y = true;
                        xmlFree(c);
                }
                else if (is_gml_element(xyz, "Z"))
                {
                        if (z) gml_lwpgerror("invalid GML representation", 23);
                        c = xmlNodeGetContent(xyz);
                        p.z = parse_gml_double((char *) c, true, true);
                        z = true;
                        xmlFree(c);
                }
        }
        /* Check dimension consistancy */
        if (!x || !y) gml_lwpgerror("invalid GML representation", 24);
        if (!z) *hasz = false;
 
        ptarray_append_point(dpa, &p, LW_FALSE);
 
        return dpa; /* ptarray_clone_deep(dpa); */
}
 
 
static POINTARRAY* parse_gml_pos(xmlNodePtr xnode, bool *hasz)
{
        xmlChar *dimension, *gmlpos;
        int dim, gml_dim;
        POINTARRAY *dpa;
        char *pos, *p;
        bool digit;
        POINT4D pt = {0, 0, 0, 0};
 
        /* HasZ, !HasM, 1 Point */
        dpa = ptarray_construct_empty(1, 0, 1);
 
    dimension = gmlGetProp(xnode, "srsDimension");
    if (dimension == NULL) /* in GML 3.0.0 it was dimension */
        dimension = gmlGetProp(xnode, "dimension");
    if (dimension == NULL) dim = 2;     /* We assume that we are in 2D */
    else
    {
        dim = atoi((char *) dimension);
        xmlFree(dimension);
        if (dim < 2 || dim > 3)
            gml_lwpgerror("invalid GML representation", 25);
    }
    if (dim == 2) *hasz = false;
 
    /* We retrieve gml:pos string */
    gmlpos = xmlNodeGetContent(xnode);
    pos = (char *) gmlpos;
    while (isspace(*pos)) pos++;        /* Eat extra whitespaces if any */
 
    /* gml:pos pattern:         x1 y1
        *                       x1 y1 z1
        */
    for (p=pos, gml_dim=0, digit=false ; *pos ; pos++)
    {
        if (isdigit(*pos)) digit = true;
        if (digit && (*pos == ' ' || *(pos+1) == '\0'))
        {
            if (*pos == ' ') *pos = '\0';
            gml_dim++;
            if  (gml_dim == 1)
                pt.x = parse_gml_double(p, true, true);
            else if (gml_dim == 2)
                pt.y = parse_gml_double(p, true, true);
            else if (gml_dim == 3)
                pt.z = parse_gml_double(p, true, true);
 
            p = pos+1;
            digit = false;
        }
    }
    xmlFree(gmlpos);
 
    /* Test again coherent dimensions on each coord */
    if (gml_dim == 2) *hasz = false;
    if (gml_dim < 2 || gml_dim > 3 || gml_dim != dim)
        gml_lwpgerror("invalid GML representation", 26);
 
    ptarray_append_point(dpa, &pt, LW_FALSE);
 
        return dpa; /* ptarray_clone_deep(dpa); */
}
 
 
static POINTARRAY* parse_gml_poslist(xmlNodePtr xnode, bool *hasz)
{
        xmlChar *dimension, *gmlposlist;
        char *poslist, *p;
        int dim, gml_dim;
        POINTARRAY *dpa;
        POINT4D pt = {0, 0, 0, 0};
        bool digit;
 
        /* Retrieve gml:srsDimension attribute if any */
        dimension = gmlGetProp(xnode, "srsDimension");
        if (dimension == NULL) /* in GML 3.0.0 it was dimension */
                dimension = gmlGetProp(xnode, "dimension");
        if (dimension == NULL) dim = 2; /* We assume that we are in common 2D */
        else
        {
                dim = atoi((char *) dimension);
                xmlFree(dimension);
                if (dim < 2 || dim > 3) gml_lwpgerror("invalid GML representation", 27);
        }
        if (dim == 2) *hasz = false;
 
        /* Retrieve gml:posList string */
        gmlposlist = xmlNodeGetContent(xnode);
        poslist = (char *) gmlposlist;
 
        /* HasZ?, !HasM, 1 point */
        dpa = ptarray_construct_empty(1, 0, 1);
 
        /* gml:posList pattern:         x1 y1 x2 y2
         *                              x1 y1 z1 x2 y2 z2
         */
        while (isspace(*poslist)) poslist++;    /* Eat extra whitespaces if any */
        for (p=poslist, gml_dim=0, digit=false ; *poslist ; poslist++)
        {
                if (isdigit(*poslist)) digit = true;
                if (digit && (*poslist == ' ' || *(poslist+1) == '\0'))
                {
                        if (*poslist == ' ') *poslist = '\0';
 
                        gml_dim++;
                        if      (gml_dim == 1) pt.x = parse_gml_double(p, true, true);
                        else if (gml_dim == 2) pt.y = parse_gml_double(p, true, true);
                        else if (gml_dim == 3) pt.z = parse_gml_double(p, true, true);
 
                        if (gml_dim == dim)
                        {
                                ptarray_append_point(dpa, &pt, LW_FALSE);
                                pt.x = pt.y = pt.z = pt.m = 0.0;
                                gml_dim = 0;
                        }
                        else if (*(poslist+1) == '\0')
                                gml_lwpgerror("invalid GML representation", 28);
 
                        p = poslist+1;
                        digit = false;
                }
        }
 
        xmlFree(gmlposlist);
 
        return dpa; /* ptarray_clone_deep(dpa); */
}
 
 
static POINTARRAY* parse_gml_data(xmlNodePtr xnode, bool *hasz, int *root_srid)
{
        POINTARRAY *pa = 0, *tmp_pa = 0;
        xmlNodePtr xa, xb;
        gmlSrs srs;
        bool found;
 
        pa = NULL;
 
        for (xa = xnode ; xa != NULL ; xa = xa->next)
        {
                if (xa->type != XML_ELEMENT_NODE) continue;
                if (!is_gml_namespace(xa, false)) continue;
                if (xa->name == NULL) continue;
 
                if (is_gml_element(xa, "pos"))
                {
                        tmp_pa = parse_gml_pos(xa, hasz);
                        if (pa == NULL) pa = tmp_pa;
                        else pa = ptarray_merge(pa, tmp_pa);
 
                }
                else if (is_gml_element(xa, "posList"))
                {
                        tmp_pa = parse_gml_poslist(xa, hasz);
                        if (pa == NULL) pa = tmp_pa;
                        else pa = ptarray_merge(pa, tmp_pa);
 
                }
                else if (is_gml_element(xa, "coordinates"))
                {
                        tmp_pa = parse_gml_coordinates(xa, hasz);
                        if (pa == NULL) pa = tmp_pa;
                        else pa = ptarray_merge(pa, tmp_pa);
 
                }
                else if (is_gml_element(xa, "coord"))
                {
                        tmp_pa = parse_gml_coord(xa, hasz);
                        if (pa == NULL) pa = tmp_pa;
                        else pa = ptarray_merge(pa, tmp_pa);
 
                }
                else if (is_gml_element(xa, "pointRep") ||
                         is_gml_element(xa, "pointProperty"))
                {
                        found = false;
                        for (xb = xa->children ; xb != NULL ; xb = xb->next)
                        {
                                if (xb->type != XML_ELEMENT_NODE) continue;
                                if (!is_gml_namespace(xb, false)) continue;
                                if (is_gml_element(xb, "Point"))
                                {
                                        found = true;
                                        break;
                                }
                        }
                        if (!found || xb == NULL)
                                gml_lwpgerror("invalid GML representation", 29);
 
                        if (is_xlink(xb)) xb = get_xlink_node(xb);
                        if (xb == NULL || xb->children == NULL)
                                gml_lwpgerror("invalid GML representation", 30);
 
                        tmp_pa = parse_gml_data(xb->children, hasz, root_srid);
                        if (tmp_pa->npoints != 1)
                                gml_lwpgerror("invalid GML representation", 31);
 
                        parse_gml_srs(xb, &srs);
                        if (srs.reverse_axis) tmp_pa = ptarray_flip_coordinates(tmp_pa);
                        if (*root_srid == SRID_UNKNOWN) *root_srid = srs.srid;
                        else if (srs.srid != *root_srid)
                                gml_reproject_pa(tmp_pa, srs.srid, *root_srid);
                        if (pa == NULL) pa = tmp_pa;
                        else pa = ptarray_merge(pa, tmp_pa);
                }
        }
 
        if (pa == NULL) gml_lwpgerror("invalid GML representation", 32);
 
        return pa;
}
 
 
static LWGEOM* parse_gml_point(xmlNodePtr xnode, bool *hasz, int *root_srid)
{
        gmlSrs srs;
        LWGEOM *geom;
        POINTARRAY *pa;
 
        if (is_xlink(xnode)) xnode = get_xlink_node(xnode);
 
        if (xnode->children == NULL)
                return lwpoint_as_lwgeom(lwpoint_construct_empty(*root_srid, 0, 0));
 
        pa = parse_gml_data(xnode->children, hasz, root_srid);
        if (pa->npoints != 1) gml_lwpgerror("invalid GML representation", 34);
 
        parse_gml_srs(xnode, &srs);
        if (srs.reverse_axis) pa = ptarray_flip_coordinates(pa);
        if (!*root_srid)
        {
                *root_srid = srs.srid;
                geom = (LWGEOM *) lwpoint_construct(*root_srid, NULL, pa);
        }
        else
        {
                if (srs.srid != *root_srid)
                        gml_reproject_pa(pa, srs.srid, *root_srid);
                geom = (LWGEOM *) lwpoint_construct(SRID_UNKNOWN, NULL, pa);
        }
 
        return geom;
}
 
 
static LWGEOM* parse_gml_line(xmlNodePtr xnode, bool *hasz, int *root_srid)
{
        gmlSrs srs;
        LWGEOM *geom;
        POINTARRAY *pa;
 
        if (is_xlink(xnode)) xnode = get_xlink_node(xnode);
 
        if (xnode->children == NULL)
                return lwline_as_lwgeom(lwline_construct_empty(*root_srid, 0, 0));
 
        pa = parse_gml_data(xnode->children, hasz, root_srid);
        if (pa->npoints < 2) gml_lwpgerror("invalid GML representation", 36);
 
        parse_gml_srs(xnode, &srs);
        if (srs.reverse_axis) pa = ptarray_flip_coordinates(pa);
        if (!*root_srid)
        {
                *root_srid = srs.srid;
                geom = (LWGEOM *) lwline_construct(*root_srid, NULL, pa);
        }
        else
        {
                if (srs.srid != *root_srid)
                        gml_reproject_pa(pa, srs.srid, *root_srid);
                geom = (LWGEOM *) lwline_construct(SRID_UNKNOWN, NULL, pa);
        }
 
        return geom;
}
 
 
static LWGEOM* parse_gml_curve(xmlNodePtr xnode, bool *hasz, int *root_srid)
{
        xmlNodePtr xa;
        size_t lss;
        bool found=false;
        gmlSrs srs;
        LWGEOM *geom=NULL;
        POINTARRAY *pa=NULL;
        POINTARRAY **ppa=NULL;
        uint32 npoints=0;
        xmlChar *interpolation=NULL;
 
        if (is_xlink(xnode)) xnode = get_xlink_node(xnode);
 
        /* Looking for gml:segments */
        for (xa = xnode->children ; xa != NULL ; xa = xa->next)
        {
                if (xa->type != XML_ELEMENT_NODE) continue;
                if (!is_gml_namespace(xa, false)) continue;
                if (is_gml_element(xa, "segments"))
                {
                        found = true;
                        break;
                }
        }
        if (!found) gml_lwpgerror("invalid GML representation", 37);
 
        ppa = (POINTARRAY**) lwalloc(sizeof(POINTARRAY*));
 
        /* Processing each gml:LineStringSegment */
        for (xa = xa->children, lss=0; xa != NULL ; xa = xa->next)
        {
                if (xa->type != XML_ELEMENT_NODE) continue;
                if (!is_gml_namespace(xa, false)) continue;
 
                if (!is_gml_element(xa, "LineStringSegment")) continue;
 
                /* GML SF is restricted to linear interpolation  */
                interpolation = gmlGetProp(xa, "interpolation");
                if (interpolation != NULL)
                {
                        if (strcmp((char *) interpolation, "linear"))
                                gml_lwpgerror("invalid GML representation", 38);
                        xmlFree(interpolation);
                }
 
                if (lss > 0) ppa = (POINTARRAY**) lwrealloc((POINTARRAY *) ppa,
                                           sizeof(POINTARRAY*) * (lss + 1));
 
                ppa[lss] = parse_gml_data(xa->children, hasz, root_srid);
                npoints += ppa[lss]->npoints;
                if (ppa[lss]->npoints < 2)
                        gml_lwpgerror("invalid GML representation", 39);
                lss++;
        }
        if (lss == 0) gml_lwpgerror("invalid GML representation", 40);
 
        /* Most common case, a single segment */
        if (lss == 1) pa = ppa[0];
 
        if (lss > 1)
        {
                /*
                 * "The curve segments are connected to one another, with the end point
                 *  of each segment except the last being the start point of the next
                 *  segment"  from  ISO 19107:2003 -> 6.3.16.1 (p43)
                 *
                 * So we must aggregate all the segments into a single one and avoid
                 * to copy the redundant points
                 */
                size_t cp_point_size = sizeof(POINT3D); /* All internals are done with 3D */
                size_t final_point_size = *hasz ? sizeof(POINT3D) : sizeof(POINT2D);
                pa = ptarray_construct(1, 0, npoints - lss + 1);
 
                /* Copy the first linestring fully */
                memcpy(getPoint_internal(pa, 0), getPoint_internal(ppa[0], 0), cp_point_size * (ppa[0]->npoints));
                npoints = ppa[0]->npoints;
                lwfree(ppa[0]);
 
                /* For the rest of linestrings, ensure the first point matches the
                 * last point of the previous one, and copy all points except the
                 * first one (since it'd be repeated)
                 */
                for (size_t i = 1; i < lss; i++)
                {
                        if (memcmp(getPoint_internal(pa, npoints - 1), getPoint_internal(ppa[i], 0), final_point_size))
                                gml_lwpgerror("invalid GML representation", 41);
 
                        memcpy(getPoint_internal(pa, npoints),
                               getPoint_internal(ppa[i], 1),
                               cp_point_size * (ppa[i]->npoints - 1));
 
                        npoints += ppa[i]->npoints - 1;
                        lwfree(ppa[i]);
                }
        }
 
        lwfree(ppa);
 
        parse_gml_srs(xnode, &srs);
        if (srs.reverse_axis) pa = ptarray_flip_coordinates(pa);
        if (srs.srid != *root_srid && *root_srid != SRID_UNKNOWN)
                gml_reproject_pa(pa, srs.srid, *root_srid);
        geom = (LWGEOM *) lwline_construct(*root_srid, NULL, pa);
 
        return geom;
}
 
 
static LWGEOM* parse_gml_linearring(xmlNodePtr xnode, bool *hasz, int *root_srid)
{
        gmlSrs srs;
        LWGEOM *geom;
        POINTARRAY **ppa = NULL;
 
        if (is_xlink(xnode)) xnode = get_xlink_node(xnode);
        parse_gml_srs(xnode, &srs);
 
        ppa = (POINTARRAY**) lwalloc(sizeof(POINTARRAY*));
        ppa[0] = parse_gml_data(xnode->children, hasz, root_srid);
 
        if (ppa[0]->npoints < 4
            || (!*hasz && !ptarray_is_closed_2d(ppa[0]))
            ||  (*hasz && !ptarray_is_closed_3d(ppa[0])))
            gml_lwpgerror("invalid GML representation", 42);
 
        if (srs.reverse_axis)
                ppa[0] = ptarray_flip_coordinates(ppa[0]);
 
        if (srs.srid != *root_srid && *root_srid != SRID_UNKNOWN)
                gml_reproject_pa(ppa[0], srs.srid, *root_srid);
 
        geom = (LWGEOM *) lwpoly_construct(*root_srid, NULL, 1, ppa);
 
        return geom;
}
 
 
static LWGEOM* parse_gml_polygon(xmlNodePtr xnode, bool *hasz, int *root_srid)
{
        gmlSrs srs;
        int i, ring;
        LWGEOM *geom;
        xmlNodePtr xa, xb;
        POINTARRAY **ppa = NULL;
 
        if (is_xlink(xnode)) xnode = get_xlink_node(xnode);
 
        if (xnode->children == NULL)
                return lwpoly_as_lwgeom(lwpoly_construct_empty(*root_srid, 0, 0));
 
        parse_gml_srs(xnode, &srs);
 
        for (xa = xnode->children ; xa != NULL ; xa = xa->next)
        {
                /* Polygon/outerBoundaryIs -> GML 2.1.2 */
                /* Polygon/exterior        -> GML 3.1.1 */
                if (xa->type != XML_ELEMENT_NODE) continue;
                if (!is_gml_namespace(xa, false)) continue;
                if (!(is_gml_element(xa, "outerBoundaryIs") ||
                      is_gml_element(xa, "exterior")))
                        continue;
 
                for (xb = xa->children ; xb != NULL ; xb = xb->next)
                {
                        if (xb->type != XML_ELEMENT_NODE) continue;
                        if (!is_gml_namespace(xb, false)) continue;
                        if (!is_gml_element(xb, "LinearRing")) continue;
 
                        ppa = (POINTARRAY**) lwalloc(sizeof(POINTARRAY*));
                        ppa[0] = parse_gml_data(xb->children, hasz, root_srid);
 
                        if (ppa[0]->npoints < 4
                                || (!*hasz && !ptarray_is_closed_2d(ppa[0]))
                                ||  (*hasz && !ptarray_is_closed_3d(ppa[0])))
                                gml_lwpgerror("invalid GML representation", 43);
 
                        if (srs.reverse_axis) ppa[0] = ptarray_flip_coordinates(ppa[0]);
                }
        }
 
        /* Found an <exterior> or <outerBoundaryIs> but no rings?!? We're outa here! */
        if ( ! ppa )
                gml_lwpgerror("invalid GML representation", 43);
 
        for (ring=1, xa = xnode->children ; xa != NULL ; xa = xa->next)
        {
                /* Polygon/innerBoundaryIs -> GML 2.1.2 */
                /* Polygon/interior        -> GML 3.1.1 */
                if (xa->type != XML_ELEMENT_NODE) continue;
                if (!is_gml_namespace(xa, false)) continue;
                if (!(is_gml_element(xa, "innerBoundaryIs") ||
                      is_gml_element(xa, "interior")))
                        continue;
 
                for (xb = xa->children ; xb != NULL ; xb = xb->next)
                {
                        if (xb->type != XML_ELEMENT_NODE) continue;
                        if (!is_gml_namespace(xb, false)) continue;
                        if (!is_gml_element(xb, "LinearRing")) continue;
 
                        ppa = (POINTARRAY**) lwrealloc((POINTARRAY *) ppa,
                                                       sizeof(POINTARRAY*) * (ring + 1));
                        ppa[ring] = parse_gml_data(xb->children, hasz, root_srid);
 
                        if (ppa[ring]->npoints < 4
                                || (!*hasz && !ptarray_is_closed_2d(ppa[ring]))
                                ||  (*hasz && !ptarray_is_closed_3d(ppa[ring])))
                                gml_lwpgerror("invalid GML representation", 43);
 
                        if (srs.reverse_axis) ppa[ring] = ptarray_flip_coordinates(ppa[ring]);
                        ring++;
                }
        }
 
        /* Exterior Ring is mandatory */
        if (ppa == NULL || ppa[0] == NULL) gml_lwpgerror("invalid GML representation", 44);
 
        if (srs.srid != *root_srid && *root_srid != SRID_UNKNOWN)
        {
                for (i=0 ; i < ring ; i++)
                        gml_reproject_pa(ppa[i], srs.srid, *root_srid);
        }
        geom = (LWGEOM *) lwpoly_construct(*root_srid, NULL, ring, ppa);
 
        return geom;
}
 
 
static LWGEOM* parse_gml_triangle(xmlNodePtr xnode, bool *hasz, int *root_srid)
{
        gmlSrs srs;
        LWGEOM *geom;
        xmlNodePtr xa, xb;
        POINTARRAY *pa = NULL;
        xmlChar *interpolation=NULL;
 
        if (is_xlink(xnode)) xnode = get_xlink_node(xnode);
 
        if (xnode->children == NULL)
                return lwtriangle_as_lwgeom(lwtriangle_construct_empty(*root_srid, 0, 0));
 
        /* GML SF is restricted to planar interpolation
               NOTA: I know Triangle is not part of SF, but
               we have to be consistent with other surfaces */
        interpolation = gmlGetProp(xnode, "interpolation");
        if (interpolation != NULL)
        {
                if (strcmp((char *) interpolation, "planar"))
                        gml_lwpgerror("invalid GML representation", 45);
                xmlFree(interpolation);
        }
 
        parse_gml_srs(xnode, &srs);
 
        for (xa = xnode->children ; xa != NULL ; xa = xa->next)
        {
                /* Triangle/exterior */
                if (xa->type != XML_ELEMENT_NODE) continue;
                if (!is_gml_namespace(xa, false)) continue;
                if (!is_gml_element(xa, "exterior")) continue;
 
                for (xb = xa->children ; xb != NULL ; xb = xb->next)
                {
                        /* Triangle/exterior/LinearRing */
                        if (xb->type != XML_ELEMENT_NODE) continue;
                        if (!is_gml_namespace(xb, false)) continue;
                        if (!is_gml_element(xb, "LinearRing")) continue;
 
                        pa = (POINTARRAY*) lwalloc(sizeof(POINTARRAY));
                        pa = parse_gml_data(xb->children, hasz, root_srid);
 
                        if (pa->npoints != 4
                                || (!*hasz && !ptarray_is_closed_2d(pa))
                                ||  (*hasz && !ptarray_is_closed_3d(pa)))
                                gml_lwpgerror("invalid GML representation", 46);
 
                        if (srs.reverse_axis) pa = ptarray_flip_coordinates(pa);
                }
        }
 
        /* Exterior Ring is mandatory */
        if (pa == NULL) gml_lwpgerror("invalid GML representation", 47);
 
        if (srs.srid != *root_srid && *root_srid != SRID_UNKNOWN)
                gml_reproject_pa(pa, srs.srid, *root_srid);
 
        geom = (LWGEOM *) lwtriangle_construct(*root_srid, NULL, pa);
 
        return geom;
}
 
 
static LWGEOM* parse_gml_patch(xmlNodePtr xnode, bool *hasz, int *root_srid)
{
        xmlChar *interpolation=NULL;
        POINTARRAY **ppa=NULL;
        LWGEOM *geom=NULL;
        xmlNodePtr xa, xb;
        int i, ring=0;
        gmlSrs srs;
 
        /* PolygonPatch */
        if (!is_gml_element(xnode, "PolygonPatch"))
                gml_lwpgerror("invalid GML representation", 48);
 
        /* GML SF is restricted to planar interpolation  */
        interpolation = gmlGetProp(xnode, "interpolation");
        if (interpolation != NULL)
        {
                if (strcmp((char *) interpolation, "planar"))
                        gml_lwpgerror("invalid GML representation", 48);
                xmlFree(interpolation);
        }
 
        parse_gml_srs(xnode, &srs);
 
        /* PolygonPatch/exterior */
        for (xa = xnode->children ; xa != NULL ; xa = xa->next)
        {
                if (!is_gml_namespace(xa, false)) continue;
                if (!is_gml_element(xa, "exterior")) continue;
 
                /* PolygonPatch/exterior/LinearRing */
                for (xb = xa->children ; xb != NULL ; xb = xb->next)
                {
                        if (xb->type != XML_ELEMENT_NODE) continue;
                        if (!is_gml_namespace(xb, false)) continue;
                        if (!is_gml_element(xb, "LinearRing")) continue;
 
                        ppa = (POINTARRAY**) lwalloc(sizeof(POINTARRAY*));
                        ppa[0] = parse_gml_data(xb->children, hasz, root_srid);
 
                        if (ppa[0]->npoints < 4
                                || (!*hasz && !ptarray_is_closed_2d(ppa[0]))
                                ||  (*hasz && !ptarray_is_closed_3d(ppa[0])))
                                gml_lwpgerror("invalid GML representation", 48);
 
                        if (srs.reverse_axis)
                                ppa[0] = ptarray_flip_coordinates(ppa[0]);
                }
        }
 
        /* Interior but no Exterior ! */
        if ( ! ppa )
                gml_lwpgerror("invalid GML representation", 48);
 
        /* PolygonPatch/interior */
        for (ring=1, xa = xnode->children ; xa != NULL ; xa = xa->next)
        {
                if (xa->type != XML_ELEMENT_NODE) continue;
                if (!is_gml_namespace(xa, false)) continue;
                if (!is_gml_element(xa, "interior")) continue;
 
                /* PolygonPatch/interior/LinearRing */
                for (xb = xa->children ; xb != NULL ; xb = xb->next)
                {
                        if (xb->type != XML_ELEMENT_NODE) continue;
                        if (!is_gml_element(xb, "LinearRing")) continue;
 
                        ppa = (POINTARRAY**) lwrealloc((POINTARRAY *) ppa,
                                                       sizeof(POINTARRAY*) * (ring + 1));
                        ppa[ring] = parse_gml_data(xb->children, hasz, root_srid);
 
                        if (ppa[ring]->npoints < 4
                                || (!*hasz && !ptarray_is_closed_2d(ppa[ring]))
                                || ( *hasz && !ptarray_is_closed_3d(ppa[ring])))
                                gml_lwpgerror("invalid GML representation", 49);
 
                        if (srs.reverse_axis)
                                ppa[ring] = ptarray_flip_coordinates(ppa[ring]);
 
                        ring++;
                }
        }
 
        /* Exterior Ring is mandatory */
        if (ppa == NULL || ppa[0] == NULL) gml_lwpgerror("invalid GML representation", 50);
 
        if (srs.srid != *root_srid && *root_srid != SRID_UNKNOWN)
        {
                for (i=0 ; i < ring ; i++)
                        gml_reproject_pa(ppa[i], srs.srid, *root_srid);
        }
        geom = (LWGEOM *) lwpoly_construct(*root_srid, NULL, ring, ppa);
 
        return geom;
}
 
 
static LWGEOM* parse_gml_surface(xmlNodePtr xnode, bool *hasz, int *root_srid)
{
        xmlNodePtr xa;
        int patch;
        LWGEOM *geom=NULL;
        bool found=false;
 
        if (is_xlink(xnode)) xnode = get_xlink_node(xnode);
 
        /* Looking for gml:patches */
        for (xa = xnode->children ; xa != NULL ; xa = xa->next)
        {
                if (xa->type != XML_ELEMENT_NODE) continue;
                if (!is_gml_namespace(xa, false)) continue;
                if (is_gml_element(xa, "patches"))
                {
                        found = true;
                        break;
                }
        }
        if (!found) gml_lwpgerror("invalid GML representation", 51);
 
        /* Processing gml:PolygonPatch */
        for (patch=0, xa = xa->children ; xa != NULL ; xa = xa->next)
        {
                if (xa->type != XML_ELEMENT_NODE) continue;
                if (!is_gml_namespace(xa, false)) continue;
                if (!is_gml_element(xa, "PolygonPatch")) continue;
                patch++;
 
                /* SQL/MM define ST_CurvePolygon as a single patch only,
                   cf ISO 13249-3:2009 -> 4.2.9 (p27) */
                if (patch > 1) gml_lwpgerror("invalid GML representation", 52);
 
                geom = parse_gml_patch(xa, hasz, root_srid);
        }
 
        if (!patch) gml_lwpgerror("invalid GML representation", 53);
 
        return geom;
}
 
 
static LWGEOM* parse_gml_tin(xmlNodePtr xnode, bool *hasz, int *root_srid)
{
        gmlSrs srs;
        xmlNodePtr xa;
        LWGEOM *geom=NULL;
        bool found=false;
 
        if (is_xlink(xnode)) xnode = get_xlink_node(xnode);
 
        parse_gml_srs(xnode, &srs);
        if (*root_srid == SRID_UNKNOWN && srs.srid != SRID_UNKNOWN)
                *root_srid = srs.srid;
 
        geom = (LWGEOM *)lwcollection_construct_empty(TINTYPE, *root_srid, 1, 0);
 
        if (xnode->children == NULL)
                return geom;
 
        /* Looking for gml:patches or gml:trianglePatches */
        for (xa = xnode->children ; xa != NULL ; xa = xa->next)
        {
                if (xa->type != XML_ELEMENT_NODE) continue;
                if (!is_gml_namespace(xa, false)) continue;
                if (is_gml_element(xa, "patches") ||
                    is_gml_element(xa, "trianglePatches"))
                {
                        found = true;
                        break;
                }
        }
        if (!found) return geom; /* empty one */
 
        /* Processing each gml:Triangle */
        for (xa = xa->children ; xa != NULL ; xa = xa->next)
        {
                if (xa->type != XML_ELEMENT_NODE) continue;
                if (!is_gml_namespace(xa, false)) continue;
                if (!is_gml_element(xa, "Triangle")) continue;
 
                if (xa->children != NULL)
                        geom = (LWGEOM*) lwtin_add_lwtriangle((LWTIN *) geom,
                               (LWTRIANGLE *) parse_gml_triangle(xa, hasz, root_srid));
        }
 
        return geom;
}
 
 
static LWGEOM* parse_gml_mpoint(xmlNodePtr xnode, bool *hasz, int *root_srid)
{
        gmlSrs srs;
        xmlNodePtr xa, xb;
        LWGEOM *geom = NULL;
 
        if (is_xlink(xnode)) xnode = get_xlink_node(xnode);
 
        parse_gml_srs(xnode, &srs);
        if (*root_srid == SRID_UNKNOWN && srs.srid != SRID_UNKNOWN)
                *root_srid = srs.srid;
 
        geom = (LWGEOM *)lwcollection_construct_empty(MULTIPOINTTYPE, *root_srid, 1, 0);
 
        if (xnode->children == NULL)
                return geom;
 
        for (xa = xnode->children ; xa != NULL ; xa = xa->next)
        {
                /* MultiPoint/pointMember */
                if (xa->type != XML_ELEMENT_NODE) continue;
                if (!is_gml_namespace(xa, false)) continue;
                if (is_gml_element(xa, "pointMembers"))
                {
                        for (xb = xa->children ; xb != NULL ; xb = xb->next)
                        {
                                geom = (LWGEOM*)lwmpoint_add_lwpoint(
                                    (LWMPOINT*)geom,
                                    (LWPOINT*)parse_gml(xb, hasz, root_srid));
                        }
                }
                else if (is_gml_element(xa, "pointMember"))
                {
                        if (xa->children != NULL)
                                geom = (LWGEOM*)lwmpoint_add_lwpoint((LWMPOINT*)geom,
                                                                 (LWPOINT*)parse_gml(xa->children, hasz, root_srid));
                }
        }
 
        return geom;
}
 
 
static LWGEOM* parse_gml_mline(xmlNodePtr xnode, bool *hasz, int *root_srid)
{
        gmlSrs srs;
        xmlNodePtr xa;
        LWGEOM *geom = NULL;
 
        if (is_xlink(xnode)) xnode = get_xlink_node(xnode);
 
        parse_gml_srs(xnode, &srs);
        if (*root_srid == SRID_UNKNOWN && srs.srid != SRID_UNKNOWN)
                *root_srid = srs.srid;
 
        geom = (LWGEOM *)lwcollection_construct_empty(MULTILINETYPE, *root_srid, 1, 0);
 
        if (xnode->children == NULL)
                return geom;
 
        for (xa = xnode->children ; xa != NULL ; xa = xa->next)
        {
                /* MultiLineString/lineStringMember */
                if (xa->type != XML_ELEMENT_NODE) continue;
                if (!is_gml_namespace(xa, false)) continue;
                if (!is_gml_element(xa, "lineStringMember")) continue;
                if (xa->children != NULL)
                        geom = (LWGEOM*)lwmline_add_lwline((LWMLINE*)geom,
                                                           (LWLINE*)parse_gml(xa->children, hasz, root_srid));
        }
 
        return geom;
}
 
 
static LWGEOM* parse_gml_mcurve(xmlNodePtr xnode, bool *hasz, int *root_srid)
{
        gmlSrs srs;
        xmlNodePtr xa, xb;
        LWGEOM *geom = NULL;
 
        if (is_xlink(xnode)) xnode = get_xlink_node(xnode);
 
        parse_gml_srs(xnode, &srs);
        if (*root_srid == SRID_UNKNOWN && srs.srid != SRID_UNKNOWN)
                *root_srid = srs.srid;
 
        geom = (LWGEOM *)lwcollection_construct_empty(MULTILINETYPE, *root_srid, 1, 0);
 
        if (xnode->children == NULL)
                return geom;
 
        for (xa = xnode->children ; xa != NULL ; xa = xa->next)
        {
 
                /* MultiCurve/curveMember */
                if (xa->type != XML_ELEMENT_NODE) continue;
                if (!is_gml_namespace(xa, false)) continue;
                if (is_gml_element(xa, "curveMembers"))
                {
                        for (xb = xa->children ; xb != NULL ; xb = xb->next)
                        {
                                if (xb != NULL)
                                        geom = (LWGEOM*)lwmline_add_lwline((LWMLINE*)geom,
                                                                       (LWLINE*)parse_gml(xb, hasz, root_srid));
                        }
                }
                else if (is_gml_element(xa, "curveMember"))
                {
                        if (xa->children != NULL)
                                geom = (LWGEOM*)lwmline_add_lwline((LWMLINE*)geom,
                                                                   (LWLINE*)parse_gml(xa->children, hasz, root_srid));
                }
        }
 
        return geom;
}
 
 
static LWGEOM* parse_gml_mpoly(xmlNodePtr xnode, bool *hasz, int *root_srid)
{
        gmlSrs srs;
        xmlNodePtr xa;
        LWGEOM *geom = NULL;
 
        if (is_xlink(xnode)) xnode = get_xlink_node(xnode);
 
        parse_gml_srs(xnode, &srs);
        if (*root_srid == SRID_UNKNOWN && srs.srid != SRID_UNKNOWN)
                *root_srid = srs.srid;
 
        geom = (LWGEOM *)lwcollection_construct_empty(MULTIPOLYGONTYPE, *root_srid, 1, 0);
 
        if (xnode->children == NULL)
                return geom;
 
        for (xa = xnode->children ; xa != NULL ; xa = xa->next)
        {
                /* MultiPolygon/polygonMember */
                if (xa->type != XML_ELEMENT_NODE) continue;
                if (!is_gml_namespace(xa, false)) continue;
                if (!is_gml_element(xa, "polygonMember")) continue;
                if (xa->children != NULL)
                        geom = (LWGEOM*)lwmpoly_add_lwpoly((LWMPOLY*)geom,
                                                           (LWPOLY*)parse_gml(xa->children, hasz, root_srid));
        }
 
        return geom;
}
 
 
static LWGEOM* parse_gml_msurface(xmlNodePtr xnode, bool *hasz, int *root_srid)
{
        gmlSrs srs;
        xmlNodePtr xa, xb;
        LWGEOM *geom = NULL;
 
        if (is_xlink(xnode)) xnode = get_xlink_node(xnode);
 
        parse_gml_srs(xnode, &srs);
        if (*root_srid == SRID_UNKNOWN && srs.srid != SRID_UNKNOWN)
                *root_srid = srs.srid;
 
        geom = (LWGEOM *)lwcollection_construct_empty(MULTIPOLYGONTYPE, *root_srid, 1, 0);
 
        if (xnode->children == NULL)
                return geom;
 
        for (xa = xnode->children ; xa != NULL ; xa = xa->next)
        {
                /* MultiSurface/surfaceMember */
                if (xa->type != XML_ELEMENT_NODE) continue;
                if (!is_gml_namespace(xa, false)) continue;
                if (is_gml_element(xa, "surfaceMembers"))
                {
                        for (xb = xa->children ; xb != NULL ; xb = xb->next)
                        {
                                if (xb != NULL)
                                        geom = (LWGEOM*)lwmpoly_add_lwpoly((LWMPOLY*)geom,
                                                                       (LWPOLY*)parse_gml(xb, hasz, root_srid));
                        }
                }
                else if (is_gml_element(xa, "surfaceMember"))
                {
                        if (xa->children != NULL)
                                geom = (LWGEOM*)lwmpoly_add_lwpoly((LWMPOLY*)geom,
                                                                   (LWPOLY*)parse_gml(xa->children, hasz, root_srid));
                }
        }
 
        return geom;
}
 
 
static LWGEOM* parse_gml_psurface(xmlNodePtr xnode, bool *hasz, int *root_srid)
{
        gmlSrs srs;
        xmlNodePtr xa;
        bool found = false;
        LWGEOM *geom = NULL;
 
        if (is_xlink(xnode)) xnode = get_xlink_node(xnode);
 
        parse_gml_srs(xnode, &srs);
        if (*root_srid == SRID_UNKNOWN && srs.srid != SRID_UNKNOWN)
                *root_srid = srs.srid;
 
        geom = (LWGEOM *)lwcollection_construct_empty(POLYHEDRALSURFACETYPE, *root_srid, 1, 0);
 
        if (xnode->children == NULL)
                return geom;
 
        /* Looking for gml:polygonPatches */
        for (xa = xnode->children ; xa != NULL ; xa = xa->next)
        {
                if (xa->type != XML_ELEMENT_NODE) continue;
                if (!is_gml_namespace(xa, false)) continue;
                if (is_gml_element(xa, "polygonPatches"))
                {
                        found = true;
                        break;
                }
        }
        if (!found) return geom;
 
        for (xa = xa->children ; xa != NULL ; xa = xa->next)
        {
                /* PolyhedralSurface/polygonPatches/PolygonPatch */
                if (xa->type != XML_ELEMENT_NODE) continue;
                if (!is_gml_namespace(xa, false)) continue;
                if (!is_gml_element(xa, "PolygonPatch")) continue;
 
                geom = (LWGEOM*)lwpsurface_add_lwpoly((LWPSURFACE*)geom,
                                                      (LWPOLY*)parse_gml_patch(xa, hasz, root_srid));
        }
 
        return geom;
}
 
 
static LWGEOM* parse_gml_coll(xmlNodePtr xnode, bool *hasz, int *root_srid)
{
        gmlSrs srs;
        xmlNodePtr xa;
        LWGEOM *geom = NULL;
 
        if (is_xlink(xnode)) xnode = get_xlink_node(xnode);
 
        parse_gml_srs(xnode, &srs);
        if (*root_srid == SRID_UNKNOWN && srs.srid != SRID_UNKNOWN)
                *root_srid = srs.srid;
 
        geom = (LWGEOM *)lwcollection_construct_empty(COLLECTIONTYPE, *root_srid, 1, 0);
 
        if (xnode->children == NULL)
                return geom;
 
        for (xa = xnode->children ; xa != NULL ; xa = xa->next)
        {
                if (xa->type != XML_ELEMENT_NODE) continue;
                if (!is_gml_namespace(xa, false)) continue;
 
                /*
                 * In GML 2.1.2 pointMember, lineStringMember and
                 * polygonMember are parts of geometryMember
                 * substitution group
                 */
                if (       is_gml_element(xa, "pointMember")
                        || is_gml_element(xa, "lineStringMember")
                        || is_gml_element(xa, "polygonMember")
                        || is_gml_element(xa, "geometryMember"))
                {
                        if (xa->children == NULL) break;
                        geom = (LWGEOM*)lwcollection_add_lwgeom((LWCOLLECTION *)geom,
                                                                parse_gml(xa->children, hasz, root_srid));
                }
        }
 
        return geom;
}
 
static LWGEOM *
lwgeom_from_gml(const char *xml, int xml_size)
{
        xmlDocPtr xmldoc;
        xmlNodePtr xmlroot=NULL;
        LWGEOM *lwgeom = NULL;
        bool hasz=true;
        int root_srid=SRID_UNKNOWN;
 
        /* Begin to Parse XML doc */
        xmlInitParser();
 
        xmldoc = xmlReadMemory(xml, xml_size, NULL, NULL, 0);
        if (!xmldoc)
        {
                xmlCleanupParser();
                gml_lwpgerror("invalid GML representation", 1);
                return NULL;
        }
 
        xmlroot = xmlDocGetRootElement(xmldoc);
        if (!xmlroot)
        {
                xmlFreeDoc(xmldoc);
                xmlCleanupParser();
                gml_lwpgerror("invalid GML representation", 1);
                return NULL;
        }
 
        lwgeom = parse_gml(xmlroot, &hasz, &root_srid);
 
        xmlFreeDoc(xmldoc);
        xmlCleanupParser();
        /* shouldn't we be releasing xmldoc too here ? */
 
        if ( root_srid != SRID_UNKNOWN )
                lwgeom->srid = root_srid;
 
        /* GML geometries could be either 2 or 3D and can be nested mixed.
         * Missing Z dimension is even tolerated inside some GML coords
         *
         * So we deal with 3D in all structures allocation, and flag hasz
         * to false if we met once a missing Z dimension
         * In this case, we force recursive 2D.
         */
        if (!hasz)
        {
                LWGEOM *tmp = lwgeom_force_2d(lwgeom);
                lwgeom_free(lwgeom);
                lwgeom = tmp;
        }
 
        return lwgeom;
}
 
 
static LWGEOM* parse_gml(xmlNodePtr xnode, bool *hasz, int *root_srid)
{
        xmlNodePtr xa = xnode;
        gmlSrs srs;
 
        /* Scroll forward to the root node */
        while (xa != NULL &&
              (xa->type != XML_ELEMENT_NODE || !is_gml_namespace(xa, false)))
        {
                xa = xa->next;
        }
 
        if (xa == NULL) gml_lwpgerror("invalid GML representation", 55);
 
        parse_gml_srs(xa, &srs);
        if (*root_srid == SRID_UNKNOWN && srs.srid != SRID_UNKNOWN)
        {
                *root_srid = srs.srid;
        }
 
        if (is_gml_element(xa, "Point"))
                return parse_gml_point(xa, hasz, root_srid);
 
        if (is_gml_element(xa, "LineString"))
                return parse_gml_line(xa, hasz, root_srid);
 
        if (is_gml_element(xa, "Curve"))
                return parse_gml_curve(xa, hasz, root_srid);
 
        if (is_gml_element(xa, "LinearRing"))
                return parse_gml_linearring(xa, hasz, root_srid);
 
        if (is_gml_element(xa, "Polygon"))
                return parse_gml_polygon(xa, hasz, root_srid);
 
        if (is_gml_element(xa, "Triangle"))
                return parse_gml_triangle(xa, hasz, root_srid);
 
        if (is_gml_element(xa, "Surface"))
                return parse_gml_surface(xa, hasz, root_srid);
 
        if (is_gml_element(xa, "MultiPoint"))
                return parse_gml_mpoint(xa, hasz, root_srid);
 
        if (is_gml_element(xa, "MultiLineString"))
                return parse_gml_mline(xa, hasz, root_srid);
 
        if (is_gml_element(xa, "MultiCurve"))
                return parse_gml_mcurve(xa, hasz, root_srid);
 
        if (is_gml_element(xa, "MultiPolygon"))
                return parse_gml_mpoly(xa, hasz, root_srid);
 
        if (is_gml_element(xa, "MultiSurface"))
                return parse_gml_msurface(xa, hasz, root_srid);
 
        if (is_gml_element(xa, "PolyhedralSurface"))
                return parse_gml_psurface(xa, hasz, root_srid);
 
        if (is_gml_element(xa, "Tin") ||
            is_gml_element(xa, "TriangulatedSurface"))
                return parse_gml_tin(xa, hasz, root_srid);
 
        if (is_gml_element(xa, "MultiGeometry"))
                return parse_gml_coll(xa, hasz, root_srid);
 
        gml_lwpgerror("invalid GML representation", 56);
        return NULL; /* Never reach */
}