rt_gdal.c

Go to the documentation of this file.

/*
 *
 * WKTRaster - Raster Types for PostGIS
 * http://trac.osgeo.org/postgis/wiki/WKTRaster
 *
 * Copyright (C) 2021 Paul Ramsey <pramsey@cleverelephant.ca>
 *
 * This program 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.
 *
 * This program 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 this program; if not, write to the Free Software Foundation,
 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 *
 */
 
#include "../../postgis_config.h"
/* #define POSTGIS_DEBUG_LEVEL 4 */
 
#include "librtcore.h"
#include "librtcore_internal.h"
#include "stringbuffer.h"
 
/******************************************************************************
* rt_raster_gdal_contour()
******************************************************************************/
 
typedef struct
{
        struct {
                GDALDatasetH ds;
                GDALDriverH drv;
                int destroy_drv;
        } src;
 
        struct {
                OGRSFDriverH drv;
                OGRDataSourceH ds;
                OGRLayerH lyr;
                int srid;
                OGRwkbGeometryType gtype;
        } dst;
 
} _rti_contour_arg;
 
 
/* ---------------------------------------------------------------- */
/*  GDAL progress callback for interrupt handling */
/* ---------------------------------------------------------------- */
 
static int rt_util_gdal_progress_func(
        double dfComplete,
        const char *pszMessage,
        void *pProgressArg)
{
        (void)dfComplete;
        (void)pszMessage;
 
        if (_lwgeom_interrupt_requested)
        {
                // rtwarn("%s interrupted at %g", (const char*)pProgressArg, dfComplete);
                _lwgeom_interrupt_requested = 0;
                return FALSE;
        }
        else
                return TRUE;
}
 
static void
_rti_contour_arg_init(_rti_contour_arg* arg)
{
        memset(arg, 0, sizeof(_rti_contour_arg));
        return;
}
 
 
static int
_rti_contour_arg_destroy(_rti_contour_arg* arg)
{
        if(arg->src.ds != NULL)
                GDALClose(arg->src.ds);
 
        if (arg->src.drv != NULL && arg->src.destroy_drv) {
                GDALDeregisterDriver(arg->src.drv);
                GDALDestroyDriver(arg->src.drv);
        }
 
        if (arg->dst.ds != NULL)
                OGR_DS_Destroy(arg->dst.ds);
 
        return FALSE;
}
 
 
 
int rt_raster_gdal_contour(
        /* input parameters */
        rt_raster src_raster,
        int src_band,
        int src_srid,
        const char* src_srs,
        double contour_interval,
        double contour_base,
        int fixed_level_count,
        double *fixed_levels,
        int polygonize,
        /* output parameters */
        size_t *ncontours,
        struct rt_contour_t **contours
        )
{
        CPLErr cplerr;
        OGRErr ogrerr;
        GDALRasterBandH hBand;
        int nfeatures = 0, i = 0;
        OGRFeatureH hFeat;
 
        /* For building out options list */
        stringbuffer_t sb;
        char **papszOptList = NULL;
        const char* elev_field = polygonize ? "ELEV_FIELD_MIN" : "ELEV_FIELD";
 
        _rti_contour_arg arg;
        _rti_contour_arg_init(&arg);
 
        /* Load raster into GDAL memory */
        arg.src.ds = rt_raster_to_gdal_mem(src_raster, src_srs, NULL, NULL, 0, &(arg.src.drv), &(arg.src.destroy_drv));
        /* Extract the desired band */
        hBand = GDALGetRasterBand(arg.src.ds, src_band);
 
        /* Set up the OGR destination data store */
        arg.dst.srid = src_srid;
        arg.dst.drv = OGRGetDriverByName("Memory");
        if (!arg.dst.drv)
                return _rti_contour_arg_destroy(&arg);
 
        arg.dst.ds = OGR_Dr_CreateDataSource(arg.dst.drv, "contour_ds", NULL);
        if (!arg.dst.ds)
                return _rti_contour_arg_destroy(&arg);
 
        /* Polygonize is 2.4+ only */
#if POSTGIS_GDAL_VERSION >= 20400
        arg.dst.gtype = polygonize ? wkbPolygon : wkbLineString;
#else
        arg.dst.gtype = wkbLineString;
#endif
 
        /* Layer has geometry, elevation, id */
        arg.dst.lyr = OGR_DS_CreateLayer(arg.dst.ds, "contours", NULL, arg.dst.gtype, NULL);
        if (!arg.dst.lyr)
                return _rti_contour_arg_destroy(&arg);
 
        /* ID field */
        OGRFieldDefnH hFldId = OGR_Fld_Create("id", OFTInteger);
        ogrerr = OGR_L_CreateField(arg.dst.lyr, hFldId, TRUE);
        if (ogrerr != OGRERR_NONE)
                return _rti_contour_arg_destroy(&arg);
 
        /* ELEVATION field */
        OGRFieldDefnH hFldElevation = OGR_Fld_Create("elevation", OFTReal);
        ogrerr = OGR_L_CreateField(arg.dst.lyr, hFldElevation, TRUE);
        if (ogrerr != OGRERR_NONE)
                return _rti_contour_arg_destroy(&arg);
 
        int use_no_data = 0;
        double no_data_value = GDALGetRasterNoDataValue(hBand, &use_no_data);
 
        // LEVEL_INTERVAL=f
        //   The elevation interval between contours generated.
        // LEVEL_BASE=f
        //   The "base" relative to which contour intervals are applied. This is normally zero, but could be different. To generate 10m contours at 5, 15, 25, ... the LEVEL_BASE would be 5.
        // LEVEL_EXP_BASE=f
        //   If greater than 0, contour levels are generated on an exponential scale. Levels will then be generated by LEVEL_EXP_BASE^k where k is a positive integer.
        // FIXED_LEVELS=f[,f]*
        //   The list of fixed contour levels at which contours should be generated. This option has precedence on LEVEL_INTERVAL
        // NODATA=f
        //   The value to use as a "nodata" value. That is, a pixel value which should be ignored in generating contours as if the value of the pixel were not known.
        // ID_FIELD=d
        //   This will be used as a field index to indicate where a unique id should be written for each feature (contour) written.
        // ELEV_FIELD=d
        //   This will be used as a field index to indicate where the elevation value of the contour should be written. Only used in line contouring mode.
        // ELEV_FIELD_MIN=d
        //   This will be used as a field index to indicate where the minimum elevation value of the polygon contour should be written. Only used in polygonal contouring mode.
        // ELEV_FIELD_MAX=d
        //   This will be used as a field index to indicate where the maximum elevation value of the polygon contour should be written. Only used in polygonal contouring mode.
        // POLYGONIZE=YES|NO
 
        /* Options strings list */
        stringbuffer_init(&sb);
 
        if (use_no_data)
                stringbuffer_aprintf(&sb, "NODATA=%g ", no_data_value);
 
        if (fixed_level_count > 0) {
                int i = 0;
                stringbuffer_append(&sb, "FIXED_LEVELS=");
                for (i = 0; i < fixed_level_count; i++) {
                        if (i) stringbuffer_append_char(&sb, ',');
                        stringbuffer_aprintf(&sb, "%g", fixed_levels[i]);
                }
                stringbuffer_append_char(&sb, ' ');
        }
        else {
                stringbuffer_aprintf(&sb, "LEVEL_INTERVAL=%g ", contour_interval);
                stringbuffer_aprintf(&sb, "LEVEL_BASE=%g ", contour_base);
        }
 
        stringbuffer_aprintf(&sb, "ID_FIELD=%d ", 0);
        stringbuffer_aprintf(&sb, "%s=%d ", elev_field, 1);
 
        stringbuffer_aprintf(&sb, "POLYGONIZE=%s ", polygonize ? "YES" : "NO");
 
        papszOptList = CSLTokenizeString(stringbuffer_getstring(&sb));
 
        // CPLSetConfigOption("OGR_GEOMETRY_ACCEPT_UNCLOSED_RING", "NO");
 
        /* Run the contouring routine, filling up the OGR layer */
        cplerr = GDALContourGenerateEx(
                hBand,
                arg.dst.lyr,
                papszOptList,      // Options
                rt_util_gdal_progress_func,  // GDALProgressFunc pfnProgress
                (void*)"GDALContourGenerateEx" // void *pProgressArg
                );
 
        // /* Run the contouring routine, filling up the OGR layer */
        // cplerr = GDALContourGenerate(
        //      hBand,
        //      contour_interval, contour_base,
        //      fixed_level_count, fixed_levels,
        //      use_no_data, no_data_value,
        //      arg.dst.lyr, 0, 1, // OGRLayer, idFieldNum, elevFieldNum
        //      NULL,              // GDALProgressFunc pfnProgress
        //      NULL               // void *pProgressArg
        //      );
 
        if (cplerr >= CE_Failure) {
                return _rti_contour_arg_destroy(&arg); // FALSE
        }
 
        /* Convert the OGR layer into PostGIS geometries */
        nfeatures = OGR_L_GetFeatureCount(arg.dst.lyr, TRUE);
        if (nfeatures < 0)
                return _rti_contour_arg_destroy(&arg);
 
        *contours = rtalloc(sizeof(struct rt_contour_t) * nfeatures);
        OGR_L_ResetReading(arg.dst.lyr);
        while ((hFeat = OGR_L_GetNextFeature(arg.dst.lyr))) {
                size_t szWkb;
                unsigned char *bufWkb;
                struct rt_contour_t contour;
                OGRGeometryH hGeom;
                LWGEOM *geom;
 
                /* Somehow we're still iterating, should not happen. */
                if(i >= nfeatures) break;
 
                /* Store elevation/id */
                contour.id = OGR_F_GetFieldAsInteger(hFeat, 0);
                contour.elevation = OGR_F_GetFieldAsDouble(hFeat, 1);
                /* Convert OGR geometry to LWGEOM via WKB */
                if (!(hGeom = OGR_F_GetGeometryRef(hFeat))) continue;
                szWkb = OGR_G_WkbSize(hGeom);
                bufWkb = rtalloc(szWkb);
                if (OGR_G_ExportToWkb(hGeom, wkbNDR, bufWkb) != OGRERR_NONE) continue;
                /* Reclaim feature and associated geometry memory */
                OGR_F_Destroy(hFeat);
                geom = lwgeom_from_wkb(bufWkb, szWkb, LW_PARSER_CHECK_NONE);
                if (!geom) rterror("%s: invalid wkb", __func__);
                lwgeom_set_srid(geom, arg.dst.srid);
                contour.geom = gserialized_from_lwgeom(geom, NULL);
                lwgeom_free(geom);
                rtdealloc(bufWkb);
                (*contours)[i++] = contour;
        }
 
        /* Return total number of contours saved */
        *ncontours = i;
 
        /* Free all the non-database allocated structures */
        _rti_contour_arg_destroy(&arg);
        stringbuffer_release(&sb);
 
        /* Done */
        return TRUE;
}