lwgeom_in_marc21.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/>.
 *
 **********************************************************************/
 
#include "postgres.h"
#include "utils/builtins.h"
 
#include <libxml/tree.h>
#include <libxml/parser.h>
#include <string.h>
 
#include "../postgis_config.h"
#include "lwgeom_pg.h"
#include <math.h>
#include "liblwgeom.h"
 
static LWGEOM* parse_marc21(xmlNodePtr xnode);
 
/**********************************************************************
 * Ability to parse geographic data contained in MARC21/XML records
 * to return an LWGEOM or an error message. It returns NULL if the
 * MARC21/XML record is valid but does not contain any geographic
 * data (datafield:034).
 *
 * MARC21/XML version supported: 1.1
 * MARC21/XML Cartographic Mathematical Data Definition:
 *    https://www.loc.gov/marc/bibliographic/bd034.html
 *
 * Copyright (C) 2021 University of Münster (WWU), Germany
 * Written by Jim Jones <jim.jones@uni-muenster.de>
 *
 **********************************************************************/
 
PG_FUNCTION_INFO_V1(ST_GeomFromMARC21);
Datum ST_GeomFromMARC21(PG_FUNCTION_ARGS) {
        GSERIALIZED *geom;
        LWGEOM *lwgeom;
        xmlDocPtr xmldoc;
        text *xml_input;
        int xml_size;
        char *xml;
        xmlNodePtr xmlroot = NULL;
 
        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);
 
        xmlInitParser();
        xmldoc = xmlReadMemory(xml, xml_size, NULL, NULL, 0);
 
        if (!xmldoc || (xmlroot = xmlDocGetRootElement(xmldoc)) == NULL) {
                xmlFreeDoc(xmldoc);
                xmlCleanupParser();
                lwpgerror("invalid MARC21/XML document.");
        }
 
        lwgeom = parse_marc21(xmlroot);
 
        xmlFreeDoc(xmldoc);
        xmlCleanupParser();
 
        if (lwgeom == NULL) {
 
                //lwgeom_free(lwgeom);
                PG_RETURN_NULL();
 
        }
 
        geom = geometry_serialize(lwgeom);
 
        lwgeom_free(lwgeom);
 
        PG_RETURN_POINTER(geom);
}
 
static inline bool
is_xml_element(xmlNodePtr xn, const char *xml_name)
{
        char *colon_pos, *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, xml_name) == 0;
}
 
static int is_literal_valid(const char *literal) {
 
        int num_dec_sep;
        int coord_start;
        int literal_length;
 
        if(literal == NULL) return LW_FALSE;
 
        literal_length = strlen(literal);
 
        POSTGIS_DEBUGF(2, "is_literal_valid called (%s)", literal);
 
        if (literal_length < 3) return LW_FALSE;
 
        coord_start = 0;
        num_dec_sep = 0;
 
        if (literal[0] == 'N' || literal[0] == 'E' || literal[0] == 'S' || literal[0] == 'W' || literal[0] == '+' || literal[0] == '-') {
 
                if (literal_length < 4) {
                        POSTGIS_DEBUGF(3, "  invalid literal length (%d): \"%s\"", literal_length, literal);
                        return LW_FALSE;
                }
 
                coord_start = 1;
        }
 
        for (int j = coord_start; j < literal_length; j++) {
 
                if (!isdigit(literal[j])) {
 
 
                        if (j < 3) {
 
                                POSTGIS_DEBUGF(3,"  invalid character '%c' at the degrees section: \"%s\"",     literal[j], literal);
                                return LW_FALSE;
 
                        }
 
                        if (literal[j] == '.' || literal[j] == ',') {
 
                                num_dec_sep++;
 
                                if (num_dec_sep > 1) return LW_FALSE;
 
                        } else {
                                POSTGIS_DEBUGF(3, "  invalid character '%c' in %d: \"%s\"",     literal[j], j, literal);
                                return LW_FALSE;
 
                        }
 
                }
 
        }
 
        POSTGIS_DEBUGF(2, "=> is_literal_valid returns LW_TRUE for \"%s\"",     literal);
        return LW_TRUE;
 
}
 
static double parse_geo_literal(char *literal) {
 
        char *dgr;
        char *min;
        char *sec;
        size_t literal_length;
 
        char start_character = literal[0];
        int start_literal = 0;
        double result = 0.0;
 
        const size_t numdigits_degrees = 3;
        const size_t numdigits_minutes = 2;
        const size_t numdigits_seconds = 2;
 
        POSTGIS_DEBUGF(2, "parse_geo_literal called (%s)", literal);
        POSTGIS_DEBUGF(2, "  start character: %c", start_character);
 
        literal_length = strlen(literal);
 
        if (!isdigit(start_character)) start_literal = 1;
 
        POSTGIS_DEBUGF(2, "    start_literal=%d", start_literal);
 
        dgr = palloc(sizeof(char)*numdigits_degrees+1);
        snprintf(dgr, numdigits_degrees+1, "%s", &literal[start_literal]);
 
        if (strchr(literal, '.') == NULL && strchr(literal, ',') == NULL) {
 
                POSTGIS_DEBUG(2, "  lat/lon integer coordinates detected");
                POSTGIS_DEBUGF(2, "    parsed degrees (lon/lat): %s", dgr);
 
                /* literal contain at least degrees.
                 * minutes and seconds are optional */
                result = atof(dgr);
 
                /* checks if the literal contains minutes */
                if (literal_length > (start_literal + numdigits_degrees)) {
 
                        min = palloc(sizeof(char)*numdigits_minutes+1);
                        snprintf(min, numdigits_minutes+1, "%s", &literal[start_literal+numdigits_degrees]);
                        POSTGIS_DEBUGF(2, "    parsed minutes (lon/lat): %s", min);
                        result = result + atof(min) / 60;
                        pfree(min);
 
                        /* checks if the literal contains seconds */
                        if (literal_length >= (start_literal + numdigits_degrees + numdigits_minutes)) {
 
                                sec = palloc(sizeof(char)*numdigits_seconds+1);
                                snprintf(sec, numdigits_seconds+1, "%s", &literal[start_literal+numdigits_degrees+numdigits_minutes]);
                                POSTGIS_DEBUGF(2, "    parsed seconds (lon/lat): %s", sec);
 
                                result = result + atof(sec) / 3600;
                                pfree(sec);
 
                        }
 
 
                }
 
 
        } else {
 
                POSTGIS_DEBUG(2, "  decimal coordinates detected");
 
                if (strchr(literal, ',')) {
 
                        /* changes the literal decimal sign from comma to period to avoid problems with atof.
                         * from the docs "In MARC21/XML coordinates, the decimal sign may be either a period or a comma." */
 
                        literal[literal_length-strlen(strchr(literal, ','))]='.';
                        POSTGIS_DEBUGF(2, "  decimal separator changed to '.': %s",literal);
 
                }
 
                /* checks if the literal is encoded in decimal degrees */
                if (literal[start_literal + numdigits_degrees] == '.') {
 
                        char *dec = palloc(sizeof(char)*literal_length+1);
                        snprintf(dec, literal_length+1, "%s", &literal[start_literal]);
                        result = atof(dec);
 
                        POSTGIS_DEBUGF(2, "    parsed decimal degrees: %s", dec);
                        pfree(dec);
 
                /* checks if the literal is encoded in decimal minutes */
                } else if (literal[start_literal + numdigits_degrees + numdigits_minutes] == '.') {
 
                        size_t len_decimal_minutes = literal_length - (start_literal + numdigits_degrees);
 
                        min = palloc(sizeof(char)*len_decimal_minutes+1);
                        snprintf(min, len_decimal_minutes+1, "%s", &literal[start_literal + numdigits_degrees]);
 
                        POSTGIS_DEBUGF(2, "    parsed degrees: %s", dgr);
                        POSTGIS_DEBUGF(2, "    parsed decimal minutes: %s", min);
 
                        result = atof(dgr) + (atof(min) / 60);
 
                        pfree(min);
 
                /* checks if the literal is encoded in decimal seconds */
                } else if (literal[start_literal + numdigits_degrees + numdigits_minutes + numdigits_seconds] == '.') {
 
                        size_t len_decimal_seconds = literal_length - (start_literal + numdigits_degrees + numdigits_minutes);
 
                        min = palloc(sizeof(char)*numdigits_minutes+1);
                        snprintf(min, numdigits_minutes+1, "%s", &literal[start_literal + numdigits_degrees]);
 
                        sec = palloc(sizeof(char)*len_decimal_seconds+1);
                        snprintf(sec, len_decimal_seconds+1, "%s", &literal[start_literal + numdigits_degrees + numdigits_minutes]);
 
                        result = atof(dgr) + (atof(min) / 60) + (atof(sec) / 3600);
 
                        POSTGIS_DEBUGF(2, "    parsed degrees: %s", dgr);
                        POSTGIS_DEBUGF(2, "    parsed minutes: %s", min);
                        POSTGIS_DEBUGF(2, "    parsed decimal seconds: %s", sec);
                        pfree(min);
                        pfree(sec);
 
                }
 
        }
 
        pfree(dgr);
 
        if (start_character == 'S' || start_character == 'W' || start_character == '-') {
 
                POSTGIS_DEBUGF(2, "  switching sign due to start character: '%c'", start_character);
                result = -result;
 
        }
 
        POSTGIS_DEBUGF(2, "=> parse_geo_literal returns: %.*f (in decimal degrees)", literal_length-(3+start_literal), result);
        return result;
}
 
static LWGEOM*
parse_marc21(xmlNodePtr xnode) {
 
        int ngeoms;
        int i;
        xmlNodePtr datafield;
        xmlNodePtr subfield;
        LWGEOM *result;
        LWGEOM **lwgeoms = (LWGEOM**) lwalloc(sizeof(LWGEOM*));
        uint8_t geometry_type;
        uint8_t result_type;
        char *code;
        char *literal;
 
        POSTGIS_DEBUGF(2, "parse_marc21 called: root '<%s>'", xnode->name);
 
        if (!is_xml_element(xnode, "record"))
                lwpgerror("invalid MARC21/XML document. Root element <record> expected but <%s> found.",xnode->name);
 
        result_type = 0;
        ngeoms = 0;
 
        for (datafield = xnode->children; datafield != NULL; datafield = datafield->next) {
 
                char *lw = NULL;
                char *le = NULL;
                char *ln = NULL;
                char *ls = NULL;
 
                if (datafield->type != XML_ELEMENT_NODE) continue;
 
                if (!is_xml_element(datafield, "datafield")     || xmlStrcmp(xmlGetProp(datafield, (xmlChar*) "tag"),(xmlChar*) "034") != 0) continue;
 
                POSTGIS_DEBUG(3, "  datafield found");
 
                for (subfield = datafield->children; subfield != NULL; subfield = subfield->next) {
 
                        if (subfield->type != XML_ELEMENT_NODE) continue;
                        if (!is_xml_element(subfield, "subfield"))
                                continue;
 
                        code = (char*) xmlGetProp(subfield, (xmlChar*) "code");
 
                        if ((strcmp(code, "d") != 0 &&
                             strcmp(code, "e") != 0 &&
                             strcmp(code, "f") != 0 &&
                             strcmp(code, "g")) != 0)
                                continue;
 
                        literal = (char*) xmlNodeGetContent(subfield);
 
                        POSTGIS_DEBUGF(3, "    subfield code '%s': %s", code, literal);
 
                        if (is_literal_valid(literal) == LW_TRUE) {
 
                                if (strcmp(code, "d") == 0) lw = literal;
                                else if (strcmp(code, "e") == 0) le = literal;
                                else if (strcmp(code, "f") == 0) ln = literal;
                                else if (strcmp(code, "g") == 0) ls = literal;
 
                        } else {
 
                                lwpgerror("parse error - invalid literal at 034$%s: \"%s\"", code, literal);
 
                        }
 
                }
 
                xmlFreeNode(subfield);
 
                if (lw && le && ln && ls) {
 
                        double w = parse_geo_literal(lw);
                        double e = parse_geo_literal(le);
                        double n = parse_geo_literal(ln);
                        double s = parse_geo_literal(ls);
                        geometry_type = 0;
 
                        if (ngeoms > 0) lwgeoms = (LWGEOM**)
                                lwrealloc(lwgeoms,      sizeof(LWGEOM*) * (ngeoms + 1));
 
                        if (fabs(w - e) < 0.0000001f && fabs(n - s) < 0.0000001f) {
 
                                lwgeoms[ngeoms] = (LWGEOM*) lwpoint_make2d(SRID_UNKNOWN, w, s);
                                geometry_type = MULTIPOINTTYPE;
 
                        } else {
 
                                lwgeoms[ngeoms] = (LWGEOM*) lwpoly_construct_envelope(SRID_UNKNOWN, w, n, e, s);
                                geometry_type = MULTIPOLYGONTYPE;
 
                        }
 
                        if (ngeoms && result_type != geometry_type) {
                                result_type = COLLECTIONTYPE;
                        } else {
                                result_type = geometry_type;
                        }
 
                        ngeoms++;
 
                } else {
 
                        if (lw || le || ln || ls) {
 
                                lwpgerror("parse error - the Coded Cartographic Mathematical Data (datafield:034) in the given MARC21/XML is incomplete. Coordinates for subfields \"$d\",\"$e\",\"$f\" and \"$g\" are expected.");
                        }
 
                }
 
        }
 
        POSTGIS_DEBUG(5, "  xmlFreeNode(datafield)");
        xmlFreeNode(datafield);
 
        if (ngeoms == 1) {
 
                POSTGIS_DEBUGF(2, "=> parse_marc21 returns single geometry: %s",lwtype_name(lwgeom_get_type(lwgeoms[0])));
                lwgeom_force_clockwise(lwgeoms[0]);
                return lwgeoms[0];
 
        } else if (ngeoms > 1) {
 
                result = (LWGEOM*) lwcollection_construct_empty(result_type,SRID_UNKNOWN, 0, 0);
 
                for (i = 0; i < ngeoms; i++) {
 
                        POSTGIS_DEBUGF(3, "  adding geometry to result set: %s",lwtype_name(lwgeom_get_type(lwgeoms[i])));
                        lwgeom_force_clockwise(lwgeoms[i]);
                        result = (LWGEOM*) lwcollection_add_lwgeom((LWCOLLECTION*) result,lwgeoms[i]);
 
                }
 
                POSTGIS_DEBUGF(2, "=> parse_marc21 returns a collection: %s", lwtype_name(lwgeom_get_type(result)));
                return result;
 
        }
 
        POSTGIS_DEBUG(2, "=> parse_marc21 returns NULL");
        return NULL;
 
}