rt_raster.c

Go to the documentation of this file.

/*
 *
 * WKTRaster - Raster Types for PostGIS
 * http://trac.osgeo.org/postgis/wiki/WKTRaster
 *
 * Copyright (C) 2013 Bborie Park <dustymugs@gmail.com>
 * Copyright (C) 2011-2013 Regents of the University of California
 *   <bkpark@ucdavis.edu>
 * Copyright (C) 2010-2011 Jorge Arevalo <jorge.arevalo@deimos-space.com>
 * Copyright (C) 2010-2011 David Zwarg <dzwarg@azavea.com>
 * Copyright (C) 2009-2011 Pierre Racine <pierre.racine@sbf.ulaval.ca>
 * Copyright (C) 2009-2011 Mateusz Loskot <mateusz@loskot.net>
 * Copyright (C) 2008-2009 Sandro Santilli <strk@kbt.io>
 *
 * 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 "librtcore.h"
#include "librtcore_internal.h"
 
#include <math.h>
 
#ifndef NAN
# define NAN 0.0/0.0
#endif
 
rt_raster
rt_raster_new(uint32_t width, uint32_t height) {
        rt_raster ret = NULL;
 
        ret = (rt_raster) rtalloc(sizeof (struct rt_raster_t));
        if (!ret) {
                rterror("rt_raster_new: Out of virtual memory creating an rt_raster");
                return NULL;
        }
 
        RASTER_DEBUGF(3, "Created rt_raster @ %p", ret);
 
        if (width > 65535 || height > 65535) {
                rterror("rt_raster_new: Dimensions requested exceed the maximum (65535 x 65535) permitted for a raster");
                rt_raster_destroy(ret);
                return NULL;
        }
 
        ret->width = width;
        ret->height = height;
        ret->scaleX = 1;
        ret->scaleY = -1;
        ret->ipX = 0.0;
        ret->ipY = 0.0;
        ret->skewX = 0.0;
        ret->skewY = 0.0;
        ret->srid = SRID_UNKNOWN;
 
        ret->numBands = 0;
        ret->bands = NULL;
 
        return ret;
}
 
void
rt_raster_destroy(rt_raster raster) {
        if (raster == NULL)
                return;
 
        RASTER_DEBUGF(3, "Destroying rt_raster @ %p", raster);
 
        if (raster->bands)
                rtdealloc(raster->bands);
 
        rtdealloc(raster);
}
 
static void
_rt_raster_geotransform_warn_offline_band(rt_raster raster) {
        int numband = 0;
        int i = 0;
        rt_band band = NULL;
 
        if (raster == NULL)
                return;
 
        numband = rt_raster_get_num_bands(raster);
        if (numband < 1)
                return;
 
        for (i = 0; i < numband; i++) {
                band = rt_raster_get_band(raster, i);
                if (NULL == band)
                        continue;
 
                if (!rt_band_is_offline(band))
                        continue;
 
                rtwarn("Changes made to raster geotransform matrix may affect out-db band data. Returned band data may be incorrect");
                break;
        }
}
 
uint16_t
rt_raster_get_width(rt_raster raster) {
 
    assert(NULL != raster);
 
    return raster->width;
}
 
uint16_t
rt_raster_get_height(rt_raster raster) {
 
    assert(NULL != raster);
 
    return raster->height;
}
 
void
rt_raster_set_scale(
        rt_raster raster,
        double scaleX, double scaleY
) {
        assert(NULL != raster);
 
        raster->scaleX = scaleX;
        raster->scaleY = scaleY;
 
        _rt_raster_geotransform_warn_offline_band(raster);
}
 
double
rt_raster_get_x_scale(rt_raster raster) {
 
 
    assert(NULL != raster);
 
    return raster->scaleX;
}
 
double
rt_raster_get_y_scale(rt_raster raster) {
 
 
    assert(NULL != raster);
 
    return raster->scaleY;
}
 
void
rt_raster_set_skews(
        rt_raster raster,
        double skewX, double skewY
) {
        assert(NULL != raster);
 
        raster->skewX = skewX;
        raster->skewY = skewY;
 
        _rt_raster_geotransform_warn_offline_band(raster);
}
 
double
rt_raster_get_x_skew(rt_raster raster) {
 
 
    assert(NULL != raster);
 
    return raster->skewX;
}
 
double
rt_raster_get_y_skew(rt_raster raster) {
 
 
    assert(NULL != raster);
 
    return raster->skewY;
}
 
void
rt_raster_set_offsets(
        rt_raster raster,
        double x, double y
) {
 
        assert(NULL != raster);
 
        raster->ipX = x;
        raster->ipY = y;
 
        _rt_raster_geotransform_warn_offline_band(raster);
}
 
double
rt_raster_get_x_offset(rt_raster raster) {
 
 
    assert(NULL != raster);
 
    return raster->ipX;
}
 
double
rt_raster_get_y_offset(rt_raster raster) {
 
 
    assert(NULL != raster);
 
    return raster->ipY;
}
 
void
rt_raster_get_phys_params(rt_raster rast,
        double *i_mag, double *j_mag, double *theta_i, double *theta_ij)
{
    double o11, o12, o21, o22 ; /* geotransform coefficients */
 
    if (rast == NULL) return ;
    if ( (i_mag==NULL) || (j_mag==NULL) || (theta_i==NULL) || (theta_ij==NULL))
        return ;
 
    /* retrieve coefficients from raster */
    o11 = rt_raster_get_x_scale(rast) ;
    o12 = rt_raster_get_x_skew(rast) ;
    o21 = rt_raster_get_y_skew(rast) ;
    o22 = rt_raster_get_y_scale(rast) ;
 
    rt_raster_calc_phys_params(o11, o12, o21, o22, i_mag, j_mag, theta_i, theta_ij);
}
 
void
rt_raster_calc_phys_params(double xscale, double xskew, double yskew, double yscale,
                           double *i_mag, double *j_mag, double *theta_i, double *theta_ij)
 
{
    double theta_test ;
 
    if ( (i_mag==NULL) || (j_mag==NULL) || (theta_i==NULL) || (theta_ij==NULL))
        return ;
 
    /* pixel size in the i direction */
    *i_mag = sqrt(xscale*xscale + yskew*yskew) ;
 
    /* pixel size in the j direction */
    *j_mag = sqrt(xskew*xskew + yscale*yscale) ;
 
    /* Rotation
     * ========
     * Two steps:
     * 1] calculate the magnitude of the angle between the x axis and
     *     the i basis vector.
     * 2] Calculate the sign of theta_i based on the angle between the y axis
     *     and the i basis vector.
     */
    *theta_i = acos(xscale/(*i_mag)) ;  /* magnitude */
    theta_test = acos(yskew/(*i_mag)) ; /* sign */
    if (theta_test < M_PI_2){
        *theta_i = -(*theta_i) ;
    }
 
 
    /* Angular separation of basis vectors
     * ===================================
     * Two steps:
     * 1] calculate the magnitude of the angle between the j basis vector and
     *     the i basis vector.
     * 2] Calculate the sign of theta_ij based on the angle between the
     *    perpendicular of the i basis vector and the j basis vector.
     */
    *theta_ij = acos(((xscale*xskew) + (yskew*yscale))/((*i_mag)*(*j_mag))) ;
    theta_test = acos( ((-yskew*xskew)+(xscale*yscale)) /
            ((*i_mag)*(*j_mag)));
    if (theta_test > M_PI_2) {
        *theta_ij = -(*theta_ij) ;
    }
}
 
void
rt_raster_set_phys_params(rt_raster rast,double i_mag, double j_mag, double theta_i, double theta_ij)
{
    double o11, o12, o21, o22 ; /* calculated geotransform coefficients */
    int success ;
 
    if (rast == NULL) return ;
 
    success = rt_raster_calc_gt_coeff(i_mag, j_mag, theta_i, theta_ij,
                            &o11, &o12, &o21, &o22) ;
 
    if (success) {
        rt_raster_set_scale(rast, o11, o22) ;
        rt_raster_set_skews(rast, o12, o21) ;
    }
}
 
int
rt_raster_calc_gt_coeff(double i_mag, double j_mag, double theta_i, double theta_ij,
                        double *xscale, double *xskew, double *yskew, double *yscale)
{
    double f ;        /* reflection flag 1.0 or -1.0 */
    double k_i ;      /* shearing coefficient */
    double s_i, s_j ; /* scaling coefficients */
    double cos_theta_i, sin_theta_i ;
 
    if ( (xscale==NULL) || (xskew==NULL) || (yskew==NULL) || (yscale==NULL)) {
        return 0;
    }
 
    if ( (theta_ij == 0.0) || (theta_ij == M_PI)) {
        return 0;
    }
 
    /* Reflection across the i axis */
    f=1.0 ;
    if (theta_ij < 0) {
        f = -1.0;
    }
 
    /* scaling along i axis */
    s_i = i_mag ;
 
    /* shearing parallel to i axis */
    k_i = tan(f*M_PI_2 - theta_ij) ;
 
    /* scaling along j axis */
    s_j = j_mag / (sqrt(k_i*k_i + 1)) ;
 
    /* putting it altogether */
    cos_theta_i = cos(theta_i) ;
    sin_theta_i = sin(theta_i) ;
    *xscale = s_i * cos_theta_i ;
    *xskew  = k_i * s_j * f * cos_theta_i + s_j * f * sin_theta_i ;
    *yskew  = -s_i * sin_theta_i ;
    *yscale = -k_i * s_j * f * sin_theta_i + s_j * f * cos_theta_i ;
    return 1;
}
 
int32_t
rt_raster_get_srid(rt_raster raster) {
        assert(NULL != raster);
 
        return clamp_srid(raster->srid);
}
 
void
rt_raster_set_srid(rt_raster raster, int32_t srid) {
        assert(NULL != raster);
 
        raster->srid = clamp_srid(srid);
 
        _rt_raster_geotransform_warn_offline_band(raster);
}
 
uint16_t
rt_raster_get_num_bands(rt_raster raster) {
 
 
    assert(NULL != raster);
 
    return raster->numBands;
}
 
rt_band
rt_raster_get_band(rt_raster raster, int n) {
        assert(NULL != raster);
 
        if (n >= raster->numBands || n < 0)
                return NULL;
 
        return raster->bands[n];
}
 
/******************************************************************************
* rt_raster_add_band()
******************************************************************************/
 
int
rt_raster_add_band(rt_raster raster, rt_band band, int index) {
    rt_band *oldbands = NULL;
    rt_band oldband = NULL;
    rt_band tmpband = NULL;
    uint16_t i = 0;
 
    assert(NULL != raster);
                assert(NULL != band);
 
    RASTER_DEBUGF(3, "Adding band %p to raster %p", band, raster);
 
    if (band->width != raster->width || band->height != raster->height) {
        rterror("rt_raster_add_band: Can't add a %dx%d band to a %dx%d raster",
                band->width, band->height, raster->width, raster->height);
        return -1;
    }
 
    if (index > raster->numBands)
        index = raster->numBands;
 
    if (index < 0)
        index = 0;
 
    oldbands = raster->bands;
 
    RASTER_DEBUGF(3, "Oldbands at %p", oldbands);
 
    raster->bands = (rt_band*) rtrealloc(raster->bands,
            sizeof (rt_band)*(raster->numBands + 1)
            );
 
    RASTER_DEBUG(3, "Checking bands");
 
    if (NULL == raster->bands) {
        rterror("rt_raster_add_band: Out of virtual memory "
                "reallocating band pointers");
        raster->bands = oldbands;
        return -1;
    }
 
    RASTER_DEBUGF(4, "realloc returned %p", raster->bands);
 
    for (i = 0; i <= raster->numBands; ++i) {
        if (i == index) {
            oldband = raster->bands[i];
            raster->bands[i] = band;
        } else if (i > index) {
            tmpband = raster->bands[i];
            raster->bands[i] = oldband;
            oldband = tmpband;
        }
    }
 
                band->raster = raster;
 
    raster->numBands++;
 
    RASTER_DEBUGF(4, "Raster now has %d bands", raster->numBands);
 
    return index;
}
 
/******************************************************************************
* rt_raster_generate_new_band()
******************************************************************************/
 
int
rt_raster_generate_new_band(
        rt_raster raster, rt_pixtype pixtype,
        double initialvalue, uint32_t hasnodata, double nodatavalue,
        int index
) {
    rt_band band = NULL;
    int width = 0;
    int height = 0;
    int numval = 0;
    int datasize = 0;
    int oldnumbands = 0;
    int numbands = 0;
    void * mem = NULL;
    int32_t checkvalint = 0;
    uint32_t checkvaluint = 0;
    double checkvaldouble = 0;
    float checkvalfloat = 0;
    int i;
 
 
    assert(NULL != raster);
 
    /* Make sure index is in a valid range */
    oldnumbands = rt_raster_get_num_bands(raster);
    if (index < 0)
        index = 0;
    else if (index > oldnumbands + 1)
        index = oldnumbands + 1;
 
    /* Determine size of memory block to allocate and allocate it */
    width = rt_raster_get_width(raster);
    height = rt_raster_get_height(raster);
    numval = width * height;
    datasize = rt_pixtype_size(pixtype) * numval;
 
    mem = (int *)rtalloc(datasize);
    if (!mem) {
        rterror("rt_raster_generate_new_band: Could not allocate memory for band");
        return -1;
    }
 
    if (FLT_EQ(initialvalue, 0.0))
        memset(mem, 0, datasize);
    else {
        switch (pixtype)
        {
            case PT_1BB:
            {
                uint8_t *ptr = mem;
                uint8_t clamped_initval = rt_util_clamp_to_1BB(initialvalue);
                for (i = 0; i < numval; i++)
                    ptr[i] = clamped_initval;
                checkvalint = ptr[0];
                break;
            }
            case PT_2BUI:
            {
                uint8_t *ptr = mem;
                uint8_t clamped_initval = rt_util_clamp_to_2BUI(initialvalue);
                for (i = 0; i < numval; i++)
                    ptr[i] = clamped_initval;
                checkvalint = ptr[0];
                break;
            }
            case PT_4BUI:
            {
                uint8_t *ptr = mem;
                uint8_t clamped_initval = rt_util_clamp_to_4BUI(initialvalue);
                for (i = 0; i < numval; i++)
                    ptr[i] = clamped_initval;
                checkvalint = ptr[0];
                break;
            }
            case PT_8BSI:
            {
                int8_t *ptr = mem;
                int8_t clamped_initval = rt_util_clamp_to_8BSI(initialvalue);
                for (i = 0; i < numval; i++)
                    ptr[i] = clamped_initval;
                checkvalint = ptr[0];
                break;
            }
            case PT_8BUI:
            {
                uint8_t *ptr = mem;
                uint8_t clamped_initval = rt_util_clamp_to_8BUI(initialvalue);
                for (i = 0; i < numval; i++)
                    ptr[i] = clamped_initval;
                checkvalint = ptr[0];
                break;
            }
            case PT_16BSI:
            {
                int16_t *ptr = mem;
                int16_t clamped_initval = rt_util_clamp_to_16BSI(initialvalue);
                for (i = 0; i < numval; i++)
                    ptr[i] = clamped_initval;
                checkvalint = ptr[0];
                break;
            }
            case PT_16BUI:
            {
                uint16_t *ptr = mem;
                uint16_t clamped_initval = rt_util_clamp_to_16BUI(initialvalue);
                for (i = 0; i < numval; i++)
                    ptr[i] = clamped_initval;
                checkvalint = ptr[0];
                break;
            }
            case PT_32BSI:
            {
                int32_t *ptr = mem;
                int32_t clamped_initval = rt_util_clamp_to_32BSI(initialvalue);
                for (i = 0; i < numval; i++)
                    ptr[i] = clamped_initval;
                checkvalint = ptr[0];
                break;
            }
            case PT_32BUI:
            {
                uint32_t *ptr = mem;
                uint32_t clamped_initval = rt_util_clamp_to_32BUI(initialvalue);
                for (i = 0; i < numval; i++)
                    ptr[i] = clamped_initval;
                checkvaluint = ptr[0];
                break;
            }
            case PT_32BF:
            {
                float *ptr = mem;
                float clamped_initval = rt_util_clamp_to_32F(initialvalue);
                for (i = 0; i < numval; i++)
                    ptr[i] = clamped_initval;
                checkvalfloat = ptr[0];
                break;
            }
            case PT_64BF:
            {
                double *ptr = mem;
                for (i = 0; i < numval; i++)
                    ptr[i] = initialvalue;
                checkvaldouble = ptr[0];
                break;
            }
            default:
            {
                rterror("rt_raster_generate_new_band: Unknown pixeltype %d", pixtype);
                rtdealloc(mem);
                return -1;
            }
        }
    }
 
    /* Overflow checking */
    rt_util_dbl_trunc_warning(
                        initialvalue,
                        checkvalint, checkvaluint,
                        checkvalfloat, checkvaldouble,
                        pixtype
                );
 
    band = rt_band_new_inline(width, height, pixtype, hasnodata, nodatavalue, mem);
    if (! band) {
        rterror("rt_raster_generate_new_band: Could not add band to raster. Aborting");
        rtdealloc(mem);
        return -1;
    }
                rt_band_set_ownsdata_flag(band, 1); /* we DO own this data!!! */
    index = rt_raster_add_band(raster, band, index);
    numbands = rt_raster_get_num_bands(raster);
    if (numbands == oldnumbands || index == -1) {
        rterror("rt_raster_generate_new_band: Could not add band to raster. Aborting");
        rt_band_destroy(band);
    }
 
                /* set isnodata if hasnodata = TRUE and initial value = nodatavalue */
                if (hasnodata && FLT_EQ(initialvalue, nodatavalue))
                        rt_band_set_isnodata_flag(band, 1);
 
    return index;
}
 
rt_errorstate rt_raster_get_inverse_geotransform_matrix(
        rt_raster raster,
        double *gt, double *igt
) {
        double _gt[6] = {0};
 
        assert((raster != NULL || gt != NULL));
        assert(igt != NULL);
 
        if (gt == NULL)
                rt_raster_get_geotransform_matrix(raster, _gt);
        else
                memcpy(_gt, gt, sizeof(double) * 6);
 
        if (!GDALInvGeoTransform(_gt, igt)) {
                rterror("rt_raster_get_inverse_geotransform_matrix: Could not compute inverse geotransform matrix");
                return ES_ERROR;
        }
 
        return ES_NONE;
}
 
void
rt_raster_get_geotransform_matrix(rt_raster raster,
        double *gt) {
        assert(NULL != raster);
        assert(NULL != gt);
 
        gt[0] = raster->ipX;
        gt[1] = raster->scaleX;
        gt[2] = raster->skewX;
        gt[3] = raster->ipY;
        gt[4] = raster->skewY;
        gt[5] = raster->scaleY;
}
 
void
rt_raster_set_geotransform_matrix(rt_raster raster,
        double *gt) {
        assert(NULL != raster);
        assert(NULL != gt);
 
        raster->ipX = gt[0];
        raster->scaleX = gt[1];
        raster->skewX = gt[2];
        raster->ipY = gt[3];
        raster->skewY = gt[4];
        raster->scaleY = gt[5];
 
        _rt_raster_geotransform_warn_offline_band(raster);
}
 
rt_errorstate
rt_raster_cell_to_geopoint(
        rt_raster raster,
        double xr, double yr,
        double *xw, double *yw,
        double *gt
) {
        double _gt[6] = {0};
 
        assert(NULL != raster);
        assert(NULL != xw && NULL != yw);
 
        if (NULL != gt)
                memcpy(_gt, gt, sizeof(double) * 6);
 
        /* scale of matrix is not set */
        if (FLT_EQ(_gt[1], 0.0) || FLT_EQ(_gt[5], 0.0))
        {
                rt_raster_get_geotransform_matrix(raster, _gt);
        }
 
        RASTER_DEBUGF(4, "gt = (%f, %f, %f, %f, %f, %f)",
                _gt[0],
                _gt[1],
                _gt[2],
                _gt[3],
                _gt[4],
                _gt[5]
        );
 
        GDALApplyGeoTransform(_gt, xr, yr, xw, yw);
        RASTER_DEBUGF(4, "GDALApplyGeoTransform (c -> g) for (%f, %f) = (%f, %f)",
                xr, yr, *xw, *yw);
 
        return ES_NONE;
}
 
rt_errorstate
rt_raster_geopoint_to_cell(
        rt_raster raster,
        double xw, double yw,
        double *xr, double *yr,
        double *igt
) {
        double rnd = 0;
        rt_errorstate err;
 
        err = rt_raster_geopoint_to_rasterpoint(raster, xw, yw, xr, yr, igt);
        if (err != ES_NONE)
                return err;
 
        rnd = ROUND(*xr, 0);
        if (FLT_EQ(rnd, *xr))
                *xr = rnd;
        else
                *xr = floor(*xr);
 
        rnd = ROUND(*yr, 0);
        if (FLT_EQ(rnd, *yr))
                *yr = rnd;
        else
                *yr = floor(*yr);
 
        RASTER_DEBUGF(4, "Corrected GDALApplyGeoTransform (g -> c) for (%f, %f) = (%f, %f)",
                xw, yw, *xr, *yr);
 
        return ES_NONE;
}
 
rt_errorstate
rt_raster_geopoint_to_rasterpoint(
        rt_raster raster,
        double xw, double yw,
        double *xr, double *yr,
        double *igt
) {
        double _igt[6] = {0};
 
        assert(NULL != raster);
        assert(NULL != xr && NULL != yr);
 
        if (igt != NULL)
                memcpy(_igt, igt, sizeof(double) * 6);
 
        /* matrix is not set */
        if (
                FLT_EQ(_igt[0], 0.) &&
                FLT_EQ(_igt[1], 0.) &&
                FLT_EQ(_igt[2], 0.) &&
                FLT_EQ(_igt[3], 0.) &&
                FLT_EQ(_igt[4], 0.) &&
                FLT_EQ(_igt[5], 0.)
        ) {
                if (rt_raster_get_inverse_geotransform_matrix(raster, NULL, _igt) != ES_NONE) {
                        rterror("rt_raster_geopoint_to_cell: Could not get inverse geotransform matrix");
                        return ES_ERROR;
                }
        }
 
        GDALApplyGeoTransform(_igt, xw, yw, xr, yr);
        RASTER_DEBUGF(4, "GDALApplyGeoTransform (g -> c) for (%f, %f) = (%f, %f)",
                xw, yw, *xr, *yr);
 
        return ES_NONE;
}
 
 
/******************************************************************************
* rt_raster_get_envelope()
******************************************************************************/
 
rt_errorstate
rt_raster_get_envelope(
        rt_raster raster,
        rt_envelope *env
) {
        int i;
        int rtn;
        int set = 0;
        double _r[2] = {0.};
        double _w[2] = {0.};
        double _gt[6] = {0.};
 
        assert(raster != NULL);
        assert(env != NULL);
 
        rt_raster_get_geotransform_matrix(raster, _gt);
 
        for (i = 0; i < 4; i++) {
                switch (i) {
                        case 0:
                                _r[0] = 0;
                                _r[1] = 0;
                                break;
                        case 1:
                                _r[0] = 0;
                                _r[1] = raster->height;
                                break;
                        case 2:
                                _r[0] = raster->width;
                                _r[1] = raster->height;
                                break;
                        case 3:
                                _r[0] = raster->width;
                                _r[1] = 0;
                                break;
                }
 
                rtn = rt_raster_cell_to_geopoint(
                        raster,
                        _r[0], _r[1],
                        &(_w[0]), &(_w[1]),
                        _gt
                );
                if (rtn != ES_NONE) {
                        rterror("rt_raster_get_envelope: Could not compute spatial coordinates for raster pixel");
                        return ES_ERROR;
                }
 
                if (!set) {
                        set = 1;
                        env->MinX = _w[0];
                        env->MaxX = _w[0];
                        env->MinY = _w[1];
                        env->MaxY = _w[1];
                }
                else {
                        if (_w[0] < env->MinX)
                                env->MinX = _w[0];
                        else if (_w[0] > env->MaxX)
                                env->MaxX = _w[0];
 
                        if (_w[1] < env->MinY)
                                env->MinY = _w[1];
                        else if (_w[1] > env->MaxY)
                                env->MaxY = _w[1];
                }
        }
 
        return ES_NONE;
}
 
/******************************************************************************
* rt_raster_compute_skewed_raster()
******************************************************************************/
 
/*
 * Compute skewed extent that covers unskewed extent.
 *
 * @param envelope : unskewed extent of type rt_envelope
 * @param skew : pointer to 2-element array (x, y) of skew
 * @param scale : pointer to 2-element array (x, y) of scale
 * @param tolerance : value between 0 and 1 where the smaller the tolerance
 * results in an extent approaching the "minimum" skewed extent.
 * If value <= 0, tolerance = 0.1. If value > 1, tolerance = 1.
 *
 * @return skewed raster who's extent covers unskewed extent, NULL on error
 */
rt_raster
rt_raster_compute_skewed_raster(
        rt_envelope extent,
        double *skew,
        double *scale,
        double tolerance
) {
        uint32_t run = 0;
        uint32_t max_run = 1;
        double dbl_run = 0;
 
        int rtn;
        int covers = 0;
        rt_raster raster;
        double _gt[6] = {0};
        double _igt[6] = {0};
        int _d[2] = {1, -1};
        int _dlast = 0;
        int _dlastpos = 0;
        double _w[2] = {0};
        double _r[2] = {0};
        double _xy[2] = {0};
        int i;
        int j;
        int x;
        int y;
 
        LWGEOM *geom = NULL;
        GEOSGeometry *sgeom = NULL;
        GEOSGeometry *ngeom = NULL;
 
        if (
                (tolerance < 0.) ||
                FLT_EQ(tolerance, 0.)
        ) {
                tolerance = 0.1;
        }
        else if (tolerance > 1.)
                tolerance = 1;
 
        dbl_run = tolerance;
        while (dbl_run < 10) {
                dbl_run *= 10.;
                max_run *= 10;
        }
 
        /* scale must be provided */
        if (scale == NULL)
                return NULL;
        for (i = 0; i < 2; i++) {
                if (FLT_EQ(scale[i], 0.0))
                {
                        rterror("rt_raster_compute_skewed_raster: Scale cannot be zero");
                        return 0;
                }
 
                if (i < 1)
                        _gt[1] = fabs(scale[i] * tolerance);
                else
                        _gt[5] = fabs(scale[i] * tolerance);
        }
        /* conform scale-y to be negative */
        _gt[5] *= -1;
 
        /* skew not provided or skew is zero, return raster of correct dim and spatial attributes */
        if ((skew == NULL) || (FLT_EQ(skew[0], 0.0) && FLT_EQ(skew[1], 0.0)))
        {
                int _dim[2] = {
                        (int) fmax((fabs(extent.MaxX - extent.MinX) + (fabs(scale[0]) / 2.)) / fabs(scale[0]), 1),
                        (int) fmax((fabs(extent.MaxY - extent.MinY) + (fabs(scale[1]) / 2.)) / fabs(scale[1]), 1)
                };
 
                raster = rt_raster_new(_dim[0], _dim[1]);
                if (raster == NULL) {
                        rterror("rt_raster_compute_skewed_raster: Could not create output raster");
                        return NULL;
                }
 
                rt_raster_set_offsets(raster, extent.MinX, extent.MaxY);
                rt_raster_set_scale(raster, fabs(scale[0]), -1 * fabs(scale[1]));
                rt_raster_set_skews(raster, skew[0], skew[1]);
 
                return raster;
        }
 
        /* direction to shift upper-left corner */
        if (skew[0] > 0.)
                _d[0] = -1;
        if (skew[1] < 0.)
                _d[1] = 1;
 
        /* geotransform */
        _gt[0] = extent.UpperLeftX;
        _gt[2] = skew[0] * tolerance;
        _gt[3] = extent.UpperLeftY;
        _gt[4] = skew[1] * tolerance;
 
        RASTER_DEBUGF(4, "Initial geotransform: %f, %f, %f, %f, %f, %f",
                _gt[0], _gt[1], _gt[2], _gt[3], _gt[4], _gt[5]
        );
        RASTER_DEBUGF(4, "Delta: %d, %d", _d[0], _d[1]);
 
        /* simple raster */
        if ((raster = rt_raster_new(1, 1)) == NULL) {
                rterror("rt_raster_compute_skewed_raster: Out of memory allocating extent raster");
                return NULL;
        }
        rt_raster_set_geotransform_matrix(raster, _gt);
 
        /* get inverse geotransform matrix */
        if (!GDALInvGeoTransform(_gt, _igt)) {
                rterror("rt_raster_compute_skewed_raster: Could not compute inverse geotransform matrix");
                rt_raster_destroy(raster);
                return NULL;
        }
        RASTER_DEBUGF(4, "Inverse geotransform: %f, %f, %f, %f, %f, %f",
                _igt[0], _igt[1], _igt[2], _igt[3], _igt[4], _igt[5]
        );
 
        /* shift along axis */
        for (i = 0; i < 2; i++) {
                covers = 0;
                run = 0;
 
                RASTER_DEBUGF(3, "Shifting along %s axis", i < 1 ? "X" : "Y");
 
                do {
 
                        /* prevent possible infinite loop */
                        if (run > max_run) {
                                rterror("rt_raster_compute_skewed_raster: Could not compute skewed extent due to check preventing infinite loop");
                                rt_raster_destroy(raster);
                                return NULL;
                        }
 
                        /*
                                check the four corners that they are covered along the specific axis
                                pixel column should be >= 0
                        */
                        for (j = 0; j < 4; j++) {
                                switch (j) {
                                        /* upper-left */
                                        case 0:
                                                _xy[0] = extent.MinX;
                                                _xy[1] = extent.MaxY;
                                                break;
                                        /* lower-left */
                                        case 1:
                                                _xy[0] = extent.MinX;
                                                _xy[1] = extent.MinY;
                                                break;
                                        /* lower-right */
                                        case 2:
                                                _xy[0] = extent.MaxX;
                                                _xy[1] = extent.MinY;
                                                break;
                                        /* upper-right */
                                        case 3:
                                                _xy[0] = extent.MaxX;
                                                _xy[1] = extent.MaxY;
                                                break;
                                }
 
                                rtn = rt_raster_geopoint_to_cell(
                                        raster,
                                        _xy[0], _xy[1],
                                        &(_r[0]), &(_r[1]),
                                        _igt
                                );
                                if (rtn != ES_NONE) {
                                        rterror("rt_raster_compute_skewed_raster: Could not compute raster pixel for spatial coordinates");
                                        rt_raster_destroy(raster);
                                        return NULL;
                                }
 
                                RASTER_DEBUGF(4, "Point %d at cell %d x %d", j, (int) _r[0], (int) _r[1]);
 
                                /* raster doesn't cover point */
                                if ((int) _r[i] < 0) {
                                        RASTER_DEBUGF(4, "Point outside of skewed extent: %d", j);
                                        covers = 0;
 
                                        if (_dlastpos != j) {
                                                _dlast = (int) _r[i];
                                                _dlastpos = j;
                                        }
                                        else if ((int) _r[i] < _dlast) {
                                                RASTER_DEBUG(4, "Point going in wrong direction.  Reversing direction");
                                                _d[i] *= -1;
                                                _dlastpos = -1;
                                                run = 0;
                                        }
 
                                        break;
                                }
 
                                covers++;
                        }
 
                        if (!covers) {
                                x = 0;
                                y = 0;
                                if (i < 1)
                                        x = _d[i] * fabs(_r[i]);
                                else
                                        y = _d[i] * fabs(_r[i]);
 
                                rtn = rt_raster_cell_to_geopoint(
                                        raster,
                                        x, y,
                                        &(_w[0]), &(_w[1]),
                                        _gt
                                );
                                if (rtn != ES_NONE) {
                                        rterror("rt_raster_compute_skewed_raster: Could not compute spatial coordinates for raster pixel");
                                        rt_raster_destroy(raster);
                                        return NULL;
                                }
 
                                /* adjust ul */
                                if (i < 1)
                                        _gt[0] = _w[i];
                                else
                                        _gt[3] = _w[i];
                                rt_raster_set_geotransform_matrix(raster, _gt);
                                RASTER_DEBUGF(4, "Shifted geotransform: %f, %f, %f, %f, %f, %f",
                                        _gt[0], _gt[1], _gt[2], _gt[3], _gt[4], _gt[5]
                                );
 
                                /* get inverse geotransform matrix */
                                if (!GDALInvGeoTransform(_gt, _igt)) {
                                        rterror("rt_raster_compute_skewed_raster: Could not compute inverse geotransform matrix");
                                        rt_raster_destroy(raster);
                                        return NULL;
                                }
                                RASTER_DEBUGF(4, "Inverse geotransform: %f, %f, %f, %f, %f, %f",
                                        _igt[0], _igt[1], _igt[2], _igt[3], _igt[4], _igt[5]
                                );
                        }
 
                        run++;
                }
                while (!covers);
        }
 
        /* covers test */
        rtn = rt_raster_geopoint_to_cell(
                raster,
                extent.MaxX, extent.MinY,
                &(_r[0]), &(_r[1]),
                _igt
        );
        if (rtn != ES_NONE) {
                rterror("rt_raster_compute_skewed_raster: Could not compute raster pixel for spatial coordinates");
                rt_raster_destroy(raster);
                return NULL;
        }
 
        RASTER_DEBUGF(4, "geopoint %f x %f at cell %d x %d", extent.MaxX, extent.MinY, (int) _r[0], (int) _r[1]);
 
        raster->width = _r[0];
        raster->height = _r[1];
 
        /* initialize GEOS */
        initGEOS(rtinfo, lwgeom_geos_error);
 
        /* create reference LWPOLY */
        {
                LWPOLY *npoly = rt_util_envelope_to_lwpoly(extent);
                if (npoly == NULL) {
                        rterror("rt_raster_compute_skewed_raster: Could not build extent's geometry for covers test");
                        rt_raster_destroy(raster);
                        return NULL;
                }
 
                ngeom = (GEOSGeometry *) LWGEOM2GEOS(lwpoly_as_lwgeom(npoly), 0);
                lwpoly_free(npoly);
        }
 
        do {
                covers = 0;
 
                /* construct sgeom from raster */
                if ((rt_raster_get_convex_hull(raster, &geom) != ES_NONE) || geom == NULL) {
                        rterror("rt_raster_compute_skewed_raster: Could not build skewed extent's geometry for covers test");
                        GEOSGeom_destroy(ngeom);
                        rt_raster_destroy(raster);
                        return NULL;
                }
 
                sgeom = (GEOSGeometry *) LWGEOM2GEOS(geom, 0);
                lwgeom_free(geom);
 
                covers = GEOSRelatePattern(sgeom, ngeom, "******FF*");
                GEOSGeom_destroy(sgeom);
 
                if (covers == 2) {
                        rterror("rt_raster_compute_skewed_raster: Could not run covers test");
                        GEOSGeom_destroy(ngeom);
                        rt_raster_destroy(raster);
                        return NULL;
                }
 
                if (!covers)
                {
                        raster->width++;
                        raster->height++;
                }
        }
        while (!covers);
 
        RASTER_DEBUGF(4, "Skewed extent does cover normal extent with dimensions %d x %d", raster->width, raster->height);
 
        raster->width = (int) ((((double) raster->width) * fabs(_gt[1]) + fabs(scale[0] / 2.)) / fabs(scale[0]));
        raster->height = (int) ((((double) raster->height) * fabs(_gt[5]) + fabs(scale[1] / 2.)) / fabs(scale[1]));
        _gt[1] = fabs(scale[0]);
        _gt[5] = -1 * fabs(scale[1]);
        _gt[2] = skew[0];
        _gt[4] = skew[1];
        rt_raster_set_geotransform_matrix(raster, _gt);
 
        /* minimize width/height */
        for (i = 0; i < 2; i++) {
                covers = 1;
                do {
                        if (i < 1)
                                raster->width--;
                        else
                                raster->height--;
 
                        /* construct sgeom from raster */
                        if ((rt_raster_get_convex_hull(raster, &geom) != ES_NONE) || geom == NULL) {
                                rterror("rt_raster_compute_skewed_raster: Could not build skewed extent's geometry for minimizing dimensions");
                                GEOSGeom_destroy(ngeom);
                                rt_raster_destroy(raster);
                                return NULL;
                        }
 
                        sgeom = (GEOSGeometry *) LWGEOM2GEOS(geom, 0);
                        lwgeom_free(geom);
 
                        covers = GEOSRelatePattern(sgeom, ngeom, "******FF*");
                        GEOSGeom_destroy(sgeom);
 
                        if (covers == 2) {
                                rterror("rt_raster_compute_skewed_raster: Could not run covers test for minimizing dimensions");
                                GEOSGeom_destroy(ngeom);
                                rt_raster_destroy(raster);
                                return NULL;
                        }
                } while (covers);
 
                if (i < 1)
                        raster->width++;
                else
                        raster->height++;
 
        }
 
        GEOSGeom_destroy(ngeom);
 
        return raster;
}
 
int
rt_raster_is_empty(rt_raster raster) {
        return (!raster || raster->height == 0 || raster->width == 0);
}
 
int
rt_raster_has_band(rt_raster raster, int nband) {
        return !(NULL == raster || nband >= raster->numBands || nband < 0);
}
 
/******************************************************************************
* rt_raster_copy_band()
******************************************************************************/
 
int
rt_raster_copy_band(
        rt_raster torast, rt_raster fromrast,
        int fromindex, int toindex
) {
        rt_band srcband = NULL;
        rt_band dstband = NULL;
 
        assert(NULL != torast);
        assert(NULL != fromrast);
 
        /* Check raster dimensions */
        if (torast->height != fromrast->height || torast->width != fromrast->width) {
                rtwarn("rt_raster_copy_band: Attempting to add a band with different width or height");
                return -1;
        }
 
        /* Check bands limits */
        if (fromrast->numBands < 1) {
                rtwarn("rt_raster_copy_band: Second raster has no band");
                return -1;
        }
        else if (fromindex < 0) {
                rtwarn("rt_raster_copy_band: Band index for second raster < 0. Defaulted to 0");
                fromindex = 0;
        }
        else if (fromindex >= fromrast->numBands) {
                rtwarn("rt_raster_copy_band: Band index for second raster > number of bands, truncated from %u to %u", fromindex, fromrast->numBands - 1);
                fromindex = fromrast->numBands - 1;
        }
 
        if (toindex < 0) {
                rtwarn("rt_raster_copy_band: Band index for first raster < 0. Defaulted to 0");
                toindex = 0;
        }
        else if (toindex > torast->numBands) {
                rtwarn("rt_raster_copy_band: Band index for first raster > number of bands, truncated from %u to %u", toindex, torast->numBands);
                toindex = torast->numBands;
        }
 
        /* Get band from source raster */
        srcband = rt_raster_get_band(fromrast, fromindex);
 
        /* duplicate band */
        dstband = rt_band_duplicate(srcband);
 
        /* Add band to the second raster */
        return rt_raster_add_band(torast, dstband, toindex);
}
 
/******************************************************************************
* rt_raster_from_band()
******************************************************************************/
 
rt_raster
rt_raster_from_band(rt_raster raster, uint32_t *bandNums, int count) {
        rt_raster rast = NULL;
        int i = 0;
        int j = 0;
        int idx;
        int32_t flag;
        double gt[6] = {0.};
 
        assert(NULL != raster);
        assert(NULL != bandNums);
 
        RASTER_DEBUGF(3, "rt_raster_from_band: source raster has %d bands",
                rt_raster_get_num_bands(raster));
 
        /* create new raster */
        rast = rt_raster_new(raster->width, raster->height);
        if (NULL == rast) {
                rterror("rt_raster_from_band: Out of memory allocating new raster");
                return NULL;
        }
 
        /* copy raster attributes */
        rt_raster_get_geotransform_matrix(raster, gt);
        rt_raster_set_geotransform_matrix(rast, gt);
 
        /* srid */
        rt_raster_set_srid(rast, raster->srid);
 
        /* copy bands */
        for (i = 0; i < count; i++) {
                idx = bandNums[i];
                flag = rt_raster_copy_band(rast, raster, idx, i);
 
                if (flag < 0) {
                        rterror("rt_raster_from_band: Could not copy band");
                        for (j = 0; j < i; j++) rt_band_destroy(rast->bands[j]);
                        rt_raster_destroy(rast);
                        return NULL;
                }
 
                RASTER_DEBUGF(3, "rt_raster_from_band: band created at index %d",
                        flag);
        }
 
        RASTER_DEBUGF(3, "rt_raster_from_band: new raster has %d bands",
                rt_raster_get_num_bands(rast));
        return rast;
}
 
/******************************************************************************
* rt_raster_replace_band()
******************************************************************************/
 
rt_band
rt_raster_replace_band(rt_raster raster, rt_band band, int index) {
        rt_band oldband = NULL;
        assert(NULL != raster);
        assert(NULL != band);
 
        if (band->width != raster->width || band->height != raster->height) {
                rterror("rt_raster_replace_band: Band does not match raster's dimensions: %dx%d band to %dx%d raster",
                        band->width, band->height, raster->width, raster->height);
                return 0;
        }
 
        if (index >= raster->numBands || index < 0) {
                rterror("rt_raster_replace_band: Band index is not valid");
                return 0;
        }
 
        oldband = rt_raster_get_band(raster, index);
        RASTER_DEBUGF(3, "rt_raster_replace_band: old band at %p", oldband);
        RASTER_DEBUGF(3, "rt_raster_replace_band: new band at %p", band);
 
        raster->bands[index] = band;
        RASTER_DEBUGF(3, "rt_raster_replace_band: new band at %p", raster->bands[index]);
 
        band->raster = raster;
        oldband->raster = NULL;
 
        return oldband;
}
 
/******************************************************************************
* rt_raster_clone()
******************************************************************************/
 
rt_raster
rt_raster_clone(rt_raster raster, uint8_t deep) {
        rt_raster rtn = NULL;
        double gt[6] = {0};
 
        assert(NULL != raster);
 
        if (deep) {
                int numband = rt_raster_get_num_bands(raster);
                uint32_t *nband = NULL;
                int i = 0;
 
                nband = rtalloc(sizeof(uint32_t) * numband);
                if (nband == NULL) {
                        rterror("rt_raster_clone: Could not allocate memory for deep clone");
                        return NULL;
                }
                for (i = 0; i < numband; i++)
                        nband[i] = i;
 
                rtn = rt_raster_from_band(raster, nband, numband);
                rtdealloc(nband);
 
                return rtn;
        }
 
        rtn = rt_raster_new(
                rt_raster_get_width(raster),
                rt_raster_get_height(raster)
        );
        if (rtn == NULL) {
                rterror("rt_raster_clone: Could not create cloned raster");
                return NULL;
        }
 
        rt_raster_get_geotransform_matrix(raster, gt);
        rt_raster_set_geotransform_matrix(rtn, gt);
        rt_raster_set_srid(rtn, rt_raster_get_srid(raster));
 
        return rtn;
}
 
/******************************************************************************
* rt_raster_copy_to_geometry()
******************************************************************************/
 
rt_errorstate
rt_raster_copy_to_geometry(
        rt_raster raster,
        uint32_t bandnum,
        char dim,
        rt_resample_type resample,
        const LWGEOM *lwgeom_in,
        LWGEOM **lwgeom_out
        )
{
        int has_z = lwgeom_has_z(lwgeom_in);
        int has_m = lwgeom_has_m(lwgeom_in);
        LWGEOM *lwgeom;
        LWPOINTITERATOR* it;
        POINT4D p;
        double igt[6] = {0};
        rt_errorstate err;
        rt_band band = NULL;
        double nodatavalue = 0.0;
 
        /* Get the band reference and read the nodatavalue */
        band = rt_raster_get_band(raster, bandnum);
        if (!band) {
                rterror("unable to read requested band");
                return ES_ERROR;
        }
        rt_band_get_nodata(band, &nodatavalue);
 
        /* Fluff up geometry to have space for our new dimension */
        if (dim == 'z') {
                if (has_z)
                        lwgeom = lwgeom_clone(lwgeom_in);
                else if (has_m)
                        lwgeom = lwgeom_force_4d(lwgeom_in, nodatavalue, nodatavalue);
                else
                        lwgeom = lwgeom_force_3dz(lwgeom_in, nodatavalue);
        }
        else if (dim == 'm') {
                if (has_m)
                        lwgeom = lwgeom_clone(lwgeom_in);
                if (has_z)
                        lwgeom = lwgeom_force_4d(lwgeom_in, nodatavalue, nodatavalue);
                else
                        lwgeom = lwgeom_force_3dm(lwgeom_in, nodatavalue);
        }
        else {
                rterror("unknown value for dim");
                return ES_ERROR;
        }
 
        /* Read every point in the geometry */
        it = lwpointiterator_create_rw(lwgeom);
        while (lwpointiterator_has_next(it))
        {
                int nodata;
                double xr, yr, value;
                lwpointiterator_peek(it, &p);
 
                /* Convert X/Y world coordinates into raster coordinates */
                err = rt_raster_geopoint_to_rasterpoint(raster, p.x, p.y, &xr, &yr, igt);
                if (err != ES_NONE) continue;
 
                /* Read the raster value for this point */
                err = rt_band_get_pixel_resample(
                        band,
                        xr, yr,
                        resample,
                        &value, &nodata
                );
 
                if (err != ES_NONE) {
                        value = NAN;
                }
 
                /* Copy in the raster value */
                if (dim == 'z')
                        p.z = value;
                if (dim == 'm')
                        p.m = value;
 
                lwpointiterator_modify_next(it, &p);
        }
        lwpointiterator_destroy(it);
 
        if (lwgeom_out)
                *lwgeom_out = lwgeom;
        return ES_NONE;
}
 
 
/******************************************************************************
* rt_raster_to_gdal()
******************************************************************************/
 
uint8_t*
rt_raster_to_gdal(
        rt_raster raster, const char *srs,
        char *format, char **options, uint64_t *gdalsize
) {
        const char *cc;
        const char *vio;
 
        GDALDriverH src_drv = NULL;
        int destroy_src_drv = 0;
        GDALDatasetH src_ds = NULL;
 
        vsi_l_offset rtn_lenvsi;
        uint64_t rtn_len = 0;
 
        GDALDriverH rtn_drv = NULL;
        GDALDatasetH rtn_ds = NULL;
        uint8_t *rtn = NULL;
 
        assert(NULL != raster);
        assert(NULL != gdalsize);
 
        /* any supported format is possible */
        rt_util_gdal_register_all(0);
        RASTER_DEBUG(3, "loaded all supported GDAL formats");
 
        /* output format not specified */
        if (format == NULL || !strlen(format))
                format = "GTiff";
        RASTER_DEBUGF(3, "output format is %s", format);
 
        /* load raster into a GDAL MEM raster */
        src_ds = rt_raster_to_gdal_mem(raster, srs, NULL, NULL, 0, &src_drv, &destroy_src_drv);
        if (NULL == src_ds) {
                rterror("rt_raster_to_gdal: Could not convert raster to GDAL MEM format");
                return 0;
        }
 
        /* load driver */
        rtn_drv = GDALGetDriverByName(format);
        if (NULL == rtn_drv) {
                rterror("rt_raster_to_gdal: Could not load the output GDAL driver");
                GDALClose(src_ds);
                if (destroy_src_drv) GDALDestroyDriver(src_drv);
                return 0;
        }
        RASTER_DEBUG(3, "Output driver loaded");
 
        /* CreateCopy support */
        cc = GDALGetMetadataItem(rtn_drv, GDAL_DCAP_CREATECOPY, NULL);
        /* VirtualIO support */
        vio = GDALGetMetadataItem(rtn_drv, GDAL_DCAP_VIRTUALIO, NULL);
 
        if (cc == NULL || vio == NULL) {
                rterror("rt_raster_to_gdal: Output GDAL driver does not support CreateCopy and/or VirtualIO");
                GDALClose(src_ds);
                if (destroy_src_drv) GDALDestroyDriver(src_drv);
                return 0;
        }
 
        /* convert GDAL MEM raster to output format */
        RASTER_DEBUG(3, "Copying GDAL MEM raster to memory file in output format");
        rtn_ds = GDALCreateCopy(
                rtn_drv,
                "/vsimem/out.dat", /* should be fine assuming this is in a process */
                src_ds,
                FALSE, /* should copy be strictly equivelent? */
                options, /* format options */
                NULL, /* progress function */
                NULL /* progress data */
        );
 
        /* close source dataset */
        GDALClose(src_ds);
        if (destroy_src_drv) GDALDestroyDriver(src_drv);
        RASTER_DEBUG(3, "Closed GDAL MEM raster");
 
        if (NULL == rtn_ds) {
                rterror("rt_raster_to_gdal: Could not create the output GDAL dataset");
                return 0;
        }
 
        RASTER_DEBUGF(4, "dataset SRS: %s", GDALGetProjectionRef(rtn_ds));
 
        /* close dataset, this also flushes any pending writes */
        GDALClose(rtn_ds);
        RASTER_DEBUG(3, "Closed GDAL output raster");
 
        RASTER_DEBUG(3, "Done copying GDAL MEM raster to memory file in output format");
 
        /* from memory file to buffer */
        RASTER_DEBUG(3, "Copying GDAL memory file to buffer");
        rtn = VSIGetMemFileBuffer("/vsimem/out.dat", &rtn_lenvsi, TRUE);
        RASTER_DEBUG(3, "Done copying GDAL memory file to buffer");
        if (NULL == rtn) {
                rterror("rt_raster_to_gdal: Could not create the output GDAL raster");
                return 0;
        }
 
        rtn_len = (uint64_t) rtn_lenvsi;
        *gdalsize = rtn_len;
 
        return rtn;
}
 
/******************************************************************************
* rt_raster_gdal_drivers()
******************************************************************************/
 
rt_gdaldriver
rt_raster_gdal_drivers(uint32_t *drv_count, uint8_t can_write)
{
        assert(drv_count != NULL);
        uint32_t output_driver = 0;
        rt_util_gdal_register_all(0);
        uint32_t count = (uint32_t)GDALGetDriverCount();
 
        rt_gdaldriver rtn = (rt_gdaldriver)rtalloc(count * sizeof(struct rt_gdaldriver_t));
        if (!rtn)
        {
                rterror("rt_raster_gdal_drivers: Could not allocate memory for gdaldriver structure");
                *drv_count = output_driver;
                return NULL;
        }
 
        for (uint32_t i = 0; i < count; i++)
        {
                GDALDriverH *drv = GDALGetDriver(i);
 
#ifdef GDAL_DCAP_RASTER
                /* Starting with GDAL 2.0, vector drivers can also be returned */
                /* Only keep raster drivers */
                const char *is_raster;
                is_raster = GDALGetMetadataItem(drv, GDAL_DCAP_RASTER, NULL);
                if (!is_raster || !EQUAL(is_raster, "YES"))
                        continue;
#endif
 
                /* CreateCopy support */
                const char *cc = GDALGetMetadataItem(drv, GDAL_DCAP_CREATECOPY, NULL);
                if (can_write && !cc)
                        continue;
 
                /* VirtualIO support */
                const char *vio = GDALGetMetadataItem(drv, GDAL_DCAP_VIRTUALIO, NULL);
                if (can_write && !vio)
                        continue;
 
                /* we can always read what GDAL can load */
                rtn[output_driver].can_read = 1;
                /* we require CreateCopy and VirtualIO support to write to GDAL */
                rtn[output_driver].can_write = (cc != NULL && vio != NULL);
 
                /* index of driver */
                rtn[output_driver].idx = i;
 
                /* short name */
                const char *txt = GDALGetDriverShortName(drv);
                size_t txt_len = strlen(txt);
                txt_len = (txt_len + 1) * sizeof(char);
                rtn[output_driver].short_name = (char *)rtalloc(txt_len);
                memcpy(rtn[output_driver].short_name, txt, txt_len);
 
                /* long name */
                txt = GDALGetDriverLongName(drv);
                txt_len = strlen(txt);
                txt_len = (txt_len + 1) * sizeof(char);
                rtn[output_driver].long_name = (char *)rtalloc(txt_len);
                memcpy(rtn[output_driver].long_name, txt, txt_len);
 
                /* creation options */
                txt = GDALGetDriverCreationOptionList(drv);
                txt_len = strlen(txt);
                txt_len = (txt_len + 1) * sizeof(char);
                rtn[output_driver].create_options = (char *)rtalloc(txt_len);
                memcpy(rtn[output_driver].create_options, txt, txt_len);
 
                output_driver++;
        }
 
        /* free unused memory */
        rtn = rtrealloc(rtn, output_driver * sizeof(struct rt_gdaldriver_t));
        *drv_count = output_driver;
 
        return rtn;
}
 
/******************************************************************************
* rt_raster_to_gdal_mem()
******************************************************************************/
 
GDALDatasetH
rt_raster_to_gdal_mem(
        rt_raster raster,
        const char *srs,
        uint32_t *bandNums,
        int *excludeNodataValues,
        int count,
        GDALDriverH *rtn_drv, int *destroy_rtn_drv
) {
        GDALDriverH drv = NULL;
        GDALDatasetH ds = NULL;
        double gt[6] = {0.0};
        CPLErr cplerr;
        GDALDataType gdal_pt = GDT_Unknown;
        GDALRasterBandH band;
        void *pVoid;
        char *pszDataPointer;
        char szGDALOption[50];
        char *apszOptions[4];
        double nodata = 0.0;
        int allocBandNums = 0;
        int allocNodataValues = 0;
 
        int i;
        uint32_t numBands;
        uint32_t width = 0;
        uint32_t height = 0;
        rt_band rtband = NULL;
        rt_pixtype pt = PT_END;
 
        assert(NULL != raster);
        assert(NULL != rtn_drv);
        assert(NULL != destroy_rtn_drv);
 
        *destroy_rtn_drv = 0;
 
        /* store raster in GDAL MEM raster */
        if (!rt_util_gdal_driver_registered("MEM")) {
                RASTER_DEBUG(4, "Registering MEM driver");
                GDALRegister_MEM();
                *destroy_rtn_drv = 1;
        }
        drv = GDALGetDriverByName("MEM");
        if (NULL == drv) {
                rterror("rt_raster_to_gdal_mem: Could not load the MEM GDAL driver");
                return 0;
        }
        *rtn_drv = drv;
 
        /* unload driver from GDAL driver manager */
        if (*destroy_rtn_drv) {
                RASTER_DEBUG(4, "Deregistering MEM driver");
                GDALDeregisterDriver(drv);
        }
 
        width = rt_raster_get_width(raster);
        height = rt_raster_get_height(raster);
        ds = GDALCreate(
                drv, "",
                width, height,
                0, GDT_Byte, NULL
        );
        if (NULL == ds) {
                rterror("rt_raster_to_gdal_mem: Could not create a GDALDataset to convert into");
                return 0;
        }
 
        /* add geotransform */
        rt_raster_get_geotransform_matrix(raster, gt);
        cplerr = GDALSetGeoTransform(ds, gt);
        if (cplerr != CE_None) {
                rterror("rt_raster_to_gdal_mem: Could not set geotransformation");
                GDALClose(ds);
                return 0;
        }
 
        /* set spatial reference */
        if (NULL != srs && strlen(srs)) {
                char *_srs = rt_util_gdal_convert_sr(srs, 0);
                if (_srs == NULL) {
                        rterror("rt_raster_to_gdal_mem: Could not convert srs to GDAL accepted format");
                        GDALClose(ds);
                        return 0;
                }
 
                cplerr = GDALSetProjection(ds, _srs);
                CPLFree(_srs);
                if (cplerr != CE_None) {
                        rterror("rt_raster_to_gdal_mem: Could not set projection");
                        GDALClose(ds);
                        return 0;
                }
                RASTER_DEBUGF(3, "Projection set to: %s", GDALGetProjectionRef(ds));
        }
 
        /* process bandNums */
        numBands = rt_raster_get_num_bands(raster);
        if (NULL != bandNums && count > 0) {
                for (i = 0; i < count; i++) {
                        if (bandNums[i] >= numBands) {
                                rterror("rt_raster_to_gdal_mem: The band index %d is invalid", bandNums[i]);
                                GDALClose(ds);
                                return 0;
                        }
                }
        }
        else {
                count = numBands;
                bandNums = (uint32_t *) rtalloc(sizeof(uint32_t) * count);
                if (NULL == bandNums) {
                        rterror("rt_raster_to_gdal_mem: Could not allocate memory for band indices");
                        GDALClose(ds);
                        return 0;
                }
                allocBandNums = 1;
                for (i = 0; i < count; i++) bandNums[i] = i;
        }
 
        /* process exclude_nodata_values */
        if (NULL == excludeNodataValues) {
                excludeNodataValues = (int *) rtalloc(sizeof(int) * count);
                if (NULL == excludeNodataValues) {
                        rterror("rt_raster_to_gdal_mem: Could not allocate memory for NODATA flags");
                        GDALClose(ds);
                        return 0;
                }
                allocNodataValues = 1;
                for (i = 0; i < count; i++) excludeNodataValues[i] = 1;
        }
 
        /* add band(s) */
        for (i = 0; i < count; i++) {
                rtband = rt_raster_get_band(raster, bandNums[i]);
                if (NULL == rtband) {
                        rterror("rt_raster_to_gdal_mem: Could not get requested band index %d", bandNums[i]);
                        if (allocBandNums) rtdealloc(bandNums);
                        if (allocNodataValues) rtdealloc(excludeNodataValues);
                        GDALClose(ds);
                        return 0;
                }
 
                pt = rt_band_get_pixtype(rtband);
                gdal_pt = rt_util_pixtype_to_gdal_datatype(pt);
                if (gdal_pt == GDT_Unknown)
                        rtwarn("rt_raster_to_gdal_mem: Unknown pixel type for band");
 
                /*
                        For all pixel types other than PT_8BSI, set pointer to start of data
                */
                if (pt != PT_8BSI) {
                        pVoid = rt_band_get_data(rtband);
                        RASTER_DEBUGF(4, "Band data is at pos %p", pVoid);
 
                        pszDataPointer = (char *) rtalloc(20 * sizeof (char));
                        sprintf(pszDataPointer, "%p", pVoid);
                        RASTER_DEBUGF(4, "rt_raster_to_gdal_mem: szDatapointer is %p",
                                pszDataPointer);
 
                        if (strncasecmp(pszDataPointer, "0x", 2) == 0)
                                sprintf(szGDALOption, "DATAPOINTER=%s", pszDataPointer);
                        else
                                sprintf(szGDALOption, "DATAPOINTER=0x%s", pszDataPointer);
 
                        RASTER_DEBUG(3, "Storing info for GDAL MEM raster band");
 
                        apszOptions[0] = szGDALOption;
                        apszOptions[1] = NULL; /* pixel offset, not needed */
                        apszOptions[2] = NULL; /* line offset, not needed */
                        apszOptions[3] = NULL;
 
                        /* free */
                        rtdealloc(pszDataPointer);
 
                        /* add band */
                        if (GDALAddBand(ds, gdal_pt, apszOptions) == CE_Failure) {
                                rterror("rt_raster_to_gdal_mem: Could not add GDAL raster band");
                                if (allocBandNums) rtdealloc(bandNums);
                                GDALClose(ds);
                                return 0;
                        }
                }
                /*
                        PT_8BSI is special as GDAL has no equivalent pixel type.
                        Must convert 8BSI to 16BSI so create basic band
                */
                else {
                        /* add band */
                        if (GDALAddBand(ds, gdal_pt, NULL) == CE_Failure) {
                                rterror("rt_raster_to_gdal_mem: Could not add GDAL raster band");
                                if (allocBandNums) rtdealloc(bandNums);
                                if (allocNodataValues) rtdealloc(excludeNodataValues);
                                GDALClose(ds);
                                return 0;
                        }
                }
 
                /* check band count */
                if (GDALGetRasterCount(ds) != i + 1) {
                        rterror("rt_raster_to_gdal_mem: Error creating GDAL MEM raster band");
                        if (allocBandNums) rtdealloc(bandNums);
                        if (allocNodataValues) rtdealloc(excludeNodataValues);
                        GDALClose(ds);
                        return 0;
                }
 
                /* get new band */
                band = NULL;
                band = GDALGetRasterBand(ds, i + 1);
                if (NULL == band) {
                        rterror("rt_raster_to_gdal_mem: Could not get GDAL band for additional processing");
                        if (allocBandNums) rtdealloc(bandNums);
                        if (allocNodataValues) rtdealloc(excludeNodataValues);
                        GDALClose(ds);
                        return 0;
                }
 
                /* PT_8BSI requires manual setting of pixels */
                if (pt == PT_8BSI) {
                        uint32_t nXBlocks, nYBlocks;
                        int nXBlockSize, nYBlockSize;
                        uint32_t iXBlock, iYBlock;
                        int nXValid, nYValid;
                        int iX, iY;
                        int iXMax, iYMax;
 
                        int x, y, z;
                        uint32_t valueslen = 0;
                        int16_t *values = NULL;
                        double value = 0.;
 
                        /* this makes use of GDAL's "natural" blocks */
                        GDALGetBlockSize(band, &nXBlockSize, &nYBlockSize);
                        nXBlocks = (width + nXBlockSize - 1) / nXBlockSize;
                        nYBlocks = (height + nYBlockSize - 1) / nYBlockSize;
                        RASTER_DEBUGF(4, "(nXBlockSize, nYBlockSize) = (%d, %d)", nXBlockSize, nYBlockSize);
                        RASTER_DEBUGF(4, "(nXBlocks, nYBlocks) = (%d, %d)", nXBlocks, nYBlocks);
 
                        /* length is for the desired pixel type */
                        valueslen = rt_pixtype_size(PT_16BSI) * nXBlockSize * nYBlockSize;
                        values = rtalloc(valueslen);
                        if (NULL == values) {
                                rterror("rt_raster_to_gdal_mem: Could not allocate memory for GDAL band pixel values");
                                if (allocBandNums) rtdealloc(bandNums);
                                if (allocNodataValues) rtdealloc(excludeNodataValues);
                                GDALClose(ds);
                                return 0;
                        }
 
                        for (iYBlock = 0; iYBlock < nYBlocks; iYBlock++) {
                                for (iXBlock = 0; iXBlock < nXBlocks; iXBlock++) {
                                        memset(values, 0, valueslen);
 
                                        x = iXBlock * nXBlockSize;
                                        y = iYBlock * nYBlockSize;
                                        RASTER_DEBUGF(4, "(iXBlock, iYBlock) = (%d, %d)", iXBlock, iYBlock);
                                        RASTER_DEBUGF(4, "(x, y) = (%d, %d)", x, y);
 
                                        /* valid block width */
                                        if ((iXBlock + 1) * nXBlockSize > width)
                                                nXValid = width - (iXBlock * nXBlockSize);
                                        else
                                                nXValid = nXBlockSize;
 
                                        /* valid block height */
                                        if ((iYBlock + 1) * nYBlockSize > height)
                                                nYValid = height - (iYBlock * nYBlockSize);
                                        else
                                                nYValid = nYBlockSize;
 
                                        RASTER_DEBUGF(4, "(nXValid, nYValid) = (%d, %d)", nXValid, nYValid);
 
                                        /* convert 8BSI values to 16BSI */
                                        z = 0;
                                        iYMax = y + nYValid;
                                        iXMax = x + nXValid;
                                        for (iY = y; iY < iYMax; iY++)  {
                                                for (iX = x; iX < iXMax; iX++)  {
                                                        if (rt_band_get_pixel(rtband, iX, iY, &value, NULL) != ES_NONE) {
                                                                rterror("rt_raster_to_gdal_mem: Could not get pixel value to convert from 8BSI to 16BSI");
                                                                rtdealloc(values);
                                                                if (allocBandNums) rtdealloc(bandNums);
                                                                if (allocNodataValues) rtdealloc(excludeNodataValues);
                                                                GDALClose(ds);
                                                                return 0;
                                                        }
 
                                                        values[z++] = rt_util_clamp_to_16BSI(value);
                                                }
                                        }
 
                                        /* burn values */
                                        if (GDALRasterIO(
                                                band, GF_Write,
                                                x, y,
                                                nXValid, nYValid,
                                                values, nXValid, nYValid,
                                                gdal_pt,
                                                0, 0
                                        ) != CE_None) {
                                                rterror("rt_raster_to_gdal_mem: Could not write converted 8BSI to 16BSI values to GDAL band");
                                                rtdealloc(values);
                                                if (allocBandNums) rtdealloc(bandNums);
                                                if (allocNodataValues) rtdealloc(excludeNodataValues);
                                                GDALClose(ds);
                                                return 0;
                                        }
                                }
                        }
 
                        rtdealloc(values);
                }
 
                /* Add nodata value for band */
                if (rt_band_get_hasnodata_flag(rtband) != FALSE && excludeNodataValues[i]) {
                        rt_band_get_nodata(rtband, &nodata);
                        if (GDALSetRasterNoDataValue(band, nodata) != CE_None)
                                rtwarn("rt_raster_to_gdal_mem: Could not set nodata value for band");
                        RASTER_DEBUGF(3, "nodata value set to %f", GDALGetRasterNoDataValue(band, NULL));
                }
 
#if POSTGIS_DEBUG_LEVEL > 3
                {
                        GDALRasterBandH _grb = NULL;
                        double _min;
                        double _max;
                        double _mean;
                        double _stddev;
 
                        _grb = GDALGetRasterBand(ds, i + 1);
                        GDALComputeRasterStatistics(_grb, FALSE, &_min, &_max, &_mean, &_stddev, NULL, NULL);
                        RASTER_DEBUGF(4, "GDAL Band %d stats: %f, %f, %f, %f", i + 1, _min, _max, _mean, _stddev);
                }
#endif
 
        }
 
        /* necessary??? */
        GDALFlushCache(ds);
 
        if (allocBandNums) rtdealloc(bandNums);
        if (allocNodataValues) rtdealloc(excludeNodataValues);
 
        return ds;
}
 
/******************************************************************************
* rt_raster_from_gdal_dataset()
******************************************************************************/
 
rt_raster
rt_raster_from_gdal_dataset(GDALDatasetH ds) {
        rt_raster rast = NULL;
        double gt[6] = {0};
        CPLErr cplerr;
        uint32_t width = 0;
        uint32_t height = 0;
        uint32_t numBands = 0;
        uint32_t i = 0;
        char *authname = NULL;
        char *authcode = NULL;
 
        GDALRasterBandH gdband = NULL;
        GDALDataType gdpixtype = GDT_Unknown;
        rt_band band;
        int32_t idx;
        rt_pixtype pt = PT_END;
        uint32_t ptlen = 0;
        int hasnodata = 0;
        double nodataval;
 
        int x;
        int y;
 
        uint32_t nXBlocks, nYBlocks;
        int nXBlockSize, nYBlockSize;
        uint32_t iXBlock, iYBlock;
        uint32_t nXValid, nYValid;
        uint32_t iY;
 
        uint8_t *values = NULL;
        uint32_t valueslen = 0;
        uint8_t *ptr = NULL;
 
        assert(NULL != ds);
 
        /* raster size */
        width = GDALGetRasterXSize(ds);
        height = GDALGetRasterYSize(ds);
        RASTER_DEBUGF(3, "Raster dimensions (width x height): %d x %d", width, height);
 
        /* create new raster */
        RASTER_DEBUG(3, "Creating new raster");
        rast = rt_raster_new(width, height);
        if (NULL == rast) {
                rterror("rt_raster_from_gdal_dataset: Out of memory allocating new raster");
                return NULL;
        }
        RASTER_DEBUGF(3, "Created raster dimensions (width x height): %d x %d", rast->width, rast->height);
 
        /* get raster attributes */
        cplerr = GDALGetGeoTransform(ds, gt);
        if (GDALGetGeoTransform(ds, gt) != CE_None) {
                RASTER_DEBUG(4, "Using default geotransform matrix (0, 1, 0, 0, 0, -1)");
                gt[0] = 0;
                gt[1] = 1;
                gt[2] = 0;
                gt[3] = 0;
                gt[4] = 0;
                gt[5] = -1;
        }
 
        /* apply raster attributes */
        rt_raster_set_geotransform_matrix(rast, gt);
 
        RASTER_DEBUGF(3, "Raster geotransform (%f, %f, %f, %f, %f, %f)",
                gt[0], gt[1], gt[2], gt[3], gt[4], gt[5]);
 
        /* srid */
        if (rt_util_gdal_sr_auth_info(ds, &authname, &authcode) == ES_NONE) {
                if (
                        authname != NULL &&
                        strcmp(authname, "EPSG") == 0 &&
                        authcode != NULL
                ) {
                        rt_raster_set_srid(rast, atoi(authcode));
                        RASTER_DEBUGF(3, "New raster's SRID = %d", rast->srid);
                }
 
                if (authname != NULL)
                        rtdealloc(authname);
                if (authcode != NULL)
                        rtdealloc(authcode);
        }
 
        numBands = GDALGetRasterCount(ds);
 
#if POSTGIS_DEBUG_LEVEL > 3
        for (i = 1; i <= numBands; i++) {
                GDALRasterBandH _grb = NULL;
                double _min;
                double _max;
                double _mean;
                double _stddev;
 
                _grb = GDALGetRasterBand(ds, i);
                GDALComputeRasterStatistics(_grb, FALSE, &_min, &_max, &_mean, &_stddev, NULL, NULL);
                RASTER_DEBUGF(4, "GDAL Band %d stats: %f, %f, %f, %f", i, _min, _max, _mean, _stddev);
        }
#endif
 
        /* copy bands */
        for (i = 1; i <= numBands; i++) {
                RASTER_DEBUGF(3, "Processing band %d of %d", i, numBands);
                gdband = NULL;
                gdband = GDALGetRasterBand(ds, i);
                if (NULL == gdband) {
                        rterror("rt_raster_from_gdal_dataset: Could not get GDAL band");
                        rt_raster_destroy(rast);
                        return NULL;
                }
                RASTER_DEBUGF(4, "gdband @ %p", gdband);
 
                /* pixtype */
                gdpixtype = GDALGetRasterDataType(gdband);
                RASTER_DEBUGF(4, "gdpixtype, size = %s, %d", GDALGetDataTypeName(gdpixtype), GDALGetDataTypeSize(gdpixtype) / 8);
                pt = rt_util_gdal_datatype_to_pixtype(gdpixtype);
                if (pt == PT_END) {
                        rterror("rt_raster_from_gdal_dataset: Unknown pixel type for GDAL band");
                        rt_raster_destroy(rast);
                        return NULL;
                }
                ptlen = rt_pixtype_size(pt);
 
                /* size: width and height */
                width = GDALGetRasterBandXSize(gdband);
                height = GDALGetRasterBandYSize(gdband);
                RASTER_DEBUGF(3, "GDAL band dimensions (width x height): %d x %d", width, height);
 
                /* nodata */
                nodataval = GDALGetRasterNoDataValue(gdband, &hasnodata);
                RASTER_DEBUGF(3, "(hasnodata, nodataval) = (%d, %f)", hasnodata, nodataval);
 
                /* create band object */
                idx = rt_raster_generate_new_band(
                        rast, pt,
                        (hasnodata ? nodataval : 0),
                        hasnodata, nodataval, rt_raster_get_num_bands(rast)
                );
                if (idx < 0) {
                        rterror("rt_raster_from_gdal_dataset: Could not allocate memory for raster band");
                        rt_raster_destroy(rast);
                        return NULL;
                }
                band = rt_raster_get_band(rast, idx);
                RASTER_DEBUGF(3, "Created band of dimension (width x height): %d x %d", band->width, band->height);
 
                /* this makes use of GDAL's "natural" blocks */
                GDALGetBlockSize(gdband, &nXBlockSize, &nYBlockSize);
                nXBlocks = (width + nXBlockSize - 1) / nXBlockSize;
                nYBlocks = (height + nYBlockSize - 1) / nYBlockSize;
                RASTER_DEBUGF(4, "(nXBlockSize, nYBlockSize) = (%d, %d)", nXBlockSize, nYBlockSize);
                RASTER_DEBUGF(4, "(nXBlocks, nYBlocks) = (%d, %d)", nXBlocks, nYBlocks);
 
                /* allocate memory for values */
                valueslen = ptlen * nXBlockSize * nYBlockSize;
                values = rtalloc(valueslen);
                if (values == NULL) {
                        rterror("rt_raster_from_gdal_dataset: Could not allocate memory for GDAL band pixel values");
                        rt_raster_destroy(rast);
                        return NULL;
                }
                RASTER_DEBUGF(3, "values @ %p of length = %d", values, valueslen);
 
                for (iYBlock = 0; iYBlock < nYBlocks; iYBlock++) {
                        for (iXBlock = 0; iXBlock < nXBlocks; iXBlock++) {
                                x = iXBlock * nXBlockSize;
                                y = iYBlock * nYBlockSize;
                                RASTER_DEBUGF(4, "(iXBlock, iYBlock) = (%d, %d)", iXBlock, iYBlock);
                                RASTER_DEBUGF(4, "(x, y) = (%d, %d)", x, y);
 
                                memset(values, 0, valueslen);
 
                                /* valid block width */
                                if ((iXBlock + 1) * nXBlockSize > width)
                                        nXValid = width - (iXBlock * nXBlockSize);
                                else
                                        nXValid = nXBlockSize;
 
                                /* valid block height */
                                if ((iYBlock + 1) * nYBlockSize > height)
                                        nYValid = height - (iYBlock * nYBlockSize);
                                else
                                        nYValid = nYBlockSize;
 
                                RASTER_DEBUGF(4, "(nXValid, nYValid) = (%d, %d)", nXValid, nYValid);
 
                                cplerr = GDALRasterIO(
                                        gdband, GF_Read,
                                        x, y,
                                        nXValid, nYValid,
                                        values, nXValid, nYValid,
                                        gdpixtype,
                                        0, 0
                                );
                                if (cplerr != CE_None) {
                                        rterror("rt_raster_from_gdal_dataset: Could not get data from GDAL raster");
                                        rtdealloc(values);
                                        rt_raster_destroy(rast);
                                        return NULL;
                                }
 
                                /* if block width is same as raster width, shortcut */
                                if (nXBlocks == 1 && nYBlockSize > 1 && nXValid == width) {
                                        x = 0;
                                        y = nYBlockSize * iYBlock;
 
                                        RASTER_DEBUGF(4, "Setting set of pixel lines at (%d, %d) for %d pixels", x, y, nXValid * nYValid);
                                        rt_band_set_pixel_line(band, x, y, values, nXValid * nYValid);
                                }
                                else {
                                        ptr = values;
                                        x = nXBlockSize * iXBlock;
                                        for (iY = 0; iY < nYValid; iY++) {
                                                y = iY + (nYBlockSize * iYBlock);
 
                                                RASTER_DEBUGF(4, "Setting pixel line at (%d, %d) for %d pixels", x, y, nXValid);
                                                rt_band_set_pixel_line(band, x, y, ptr, nXValid);
                                                ptr += (nXValid * ptlen);
                                        }
                                }
                        }
                }
 
                /* free memory */
                rtdealloc(values);
        }
 
        return rast;
}
 
/******************************************************************************
* rt_raster_gdal_rasterize()
******************************************************************************/
 
typedef struct _rti_rasterize_arg_t* _rti_rasterize_arg;
struct _rti_rasterize_arg_t {
        uint8_t noband;
 
        uint32_t numbands;
 
        OGRSpatialReferenceH src_sr;
 
        rt_pixtype *pixtype;
        double *init;
        double *nodata;
        uint8_t *hasnodata;
        double *value;
        int *bandlist;
};
 
static _rti_rasterize_arg
_rti_rasterize_arg_init() {
        _rti_rasterize_arg arg = NULL;
 
        arg = rtalloc(sizeof(struct _rti_rasterize_arg_t));
        if (arg == NULL) {
                rterror("_rti_rasterize_arg_init: Could not allocate memory for _rti_rasterize_arg");
                return NULL;
        }
 
        arg->noband = 0;
 
        arg->numbands = 0;
 
        arg->src_sr = NULL;
 
        arg->pixtype = NULL;
        arg->init = NULL;
        arg->nodata = NULL;
        arg->hasnodata = NULL;
        arg->value = NULL;
        arg->bandlist = NULL;
 
        return arg;
}
 
static void
_rti_rasterize_arg_destroy(_rti_rasterize_arg arg) {
        if (arg->noband) {
                if (arg->pixtype != NULL)
                        rtdealloc(arg->pixtype);
                if (arg->init != NULL)
                        rtdealloc(arg->init);
                if (arg->nodata != NULL)
                        rtdealloc(arg->nodata);
                if (arg->hasnodata != NULL)
                        rtdealloc(arg->hasnodata);
                if (arg->value != NULL)
                        rtdealloc(arg->value);
        }
 
        if (arg->bandlist != NULL)
                rtdealloc(arg->bandlist);
 
        if (arg->src_sr != NULL)
                OSRDestroySpatialReference(arg->src_sr);
 
        rtdealloc(arg);
}
 
rt_raster
rt_raster_gdal_rasterize(
        const unsigned char *wkb, uint32_t wkb_len,
        const char *srs,
        uint32_t num_bands, rt_pixtype *pixtype,
        double *init, double *value,
        double *nodata, uint8_t *hasnodata,
        int *width, int *height,
        double *scale_x, double *scale_y,
        double *ul_xw, double *ul_yw,
        double *grid_xw, double *grid_yw,
        double *skew_x, double *skew_y,
        char **options
) {
        rt_raster rast = NULL;
        uint32_t i = 0;
        int err = 0;
 
        _rti_rasterize_arg arg = NULL;
 
        int _dim[2] = {0};
        double _scale[2] = {0};
        double _skew[2] = {0};
 
        OGRErr ogrerr;
        OGRGeometryH src_geom;
        OGREnvelope src_env;
        rt_envelope extent;
        OGRwkbGeometryType wkbtype = wkbUnknown;
 
        int ul_user = 0;
 
        CPLErr cplerr;
        double _gt[6] = {0};
        GDALDriverH _drv = NULL;
        int unload_drv = 0;
        GDALDatasetH _ds = NULL;
        GDALRasterBandH _band = NULL;
 
        uint16_t _width = 0;
        uint16_t _height = 0;
 
        RASTER_DEBUG(3, "starting");
 
        assert(NULL != wkb);
        assert(0 != wkb_len);
 
        /* internal variables */
        arg = _rti_rasterize_arg_init();
        if (arg == NULL) {
                rterror("rt_raster_gdal_rasterize: Could not initialize internal variables");
                return NULL;
        }
 
        /* no bands, raster is a mask */
        if (num_bands < 1) {
                arg->noband = 1;
                arg->numbands = 1;
 
                arg->pixtype = (rt_pixtype *) rtalloc(sizeof(rt_pixtype));
                arg->pixtype[0] = PT_8BUI;
 
                arg->init = (double *) rtalloc(sizeof(double));
                arg->init[0] = 0;
 
                arg->nodata = (double *) rtalloc(sizeof(double));
                arg->nodata[0] = 0;
 
                arg->hasnodata = (uint8_t *) rtalloc(sizeof(uint8_t));
                arg->hasnodata[0] = 1;
 
                arg->value = (double *) rtalloc(sizeof(double));
                arg->value[0] = 1;
        }
        else {
                arg->noband = 0;
                arg->numbands = num_bands;
 
                arg->pixtype = pixtype;
                arg->init = init;
                arg->nodata = nodata;
                arg->hasnodata = hasnodata;
                arg->value = value;
        }
 
        /* OGR spatial reference */
        if (NULL != srs && strlen(srs)) {
                arg->src_sr = OSRNewSpatialReference(NULL);
                if (OSRSetFromUserInput(arg->src_sr, srs) != OGRERR_NONE) {
                        rterror("rt_raster_gdal_rasterize: Could not create OSR spatial reference using the provided srs: %s", srs);
                        _rti_rasterize_arg_destroy(arg);
                        return NULL;
                }
        }
 
        /* convert WKB to OGR Geometry */
        ogrerr = OGR_G_CreateFromWkb((unsigned char *) wkb, arg->src_sr, &src_geom, wkb_len);
        if (ogrerr != OGRERR_NONE) {
                rterror("rt_raster_gdal_rasterize: Could not create OGR Geometry from WKB");
 
                _rti_rasterize_arg_destroy(arg);
                /* OGRCleanupAll(); */
 
                return NULL;
        }
 
        /* OGR Geometry is empty */
        if (OGR_G_IsEmpty(src_geom)) {
                rtinfo("Geometry provided is empty. Returning empty raster");
 
                OGR_G_DestroyGeometry(src_geom);
                _rti_rasterize_arg_destroy(arg);
                /* OGRCleanupAll(); */
 
                return rt_raster_new(0, 0);
        }
 
        /* get envelope */
        OGR_G_GetEnvelope(src_geom, &src_env);
        rt_util_from_ogr_envelope(src_env, &extent);
 
        RASTER_DEBUGF(3, "Suggested raster envelope: %f, %f, %f, %f",
                extent.MinX, extent.MinY, extent.MaxX, extent.MaxY);
 
        /* user-defined scale */
        if (
                (NULL != scale_x) &&
                (NULL != scale_y) &&
                (FLT_NEQ(*scale_x, 0.0)) &&
                (FLT_NEQ(*scale_y, 0.0))
        ) {
                /* for now, force scale to be in left-right, top-down orientation */
                _scale[0] = fabs(*scale_x);
                _scale[1] = fabs(*scale_y);
        }
        /* user-defined width/height */
        else if ((NULL != width) && (NULL != height) && (*width != 0) && (*height != 0))
        {
                _dim[0] = abs(*width);
                _dim[1] = abs(*height);
 
                if (FLT_NEQ(extent.MaxX, extent.MinX))
                        _scale[0] = fabs((extent.MaxX - extent.MinX) / _dim[0]);
                else
                        _scale[0] = 1.;
 
                if (FLT_NEQ(extent.MaxY, extent.MinY))
                        _scale[1] = fabs((extent.MaxY - extent.MinY) / _dim[1]);
                else
                        _scale[1] = 1.;
        }
        else {
                rterror("rt_raster_gdal_rasterize: Values must be provided for width and height or X and Y of scale");
 
                OGR_G_DestroyGeometry(src_geom);
                _rti_rasterize_arg_destroy(arg);
                /* OGRCleanupAll(); */
 
                return NULL;
        }
        RASTER_DEBUGF(3, "scale (x, y) = %f, %f", _scale[0], -1 * _scale[1]);
        RASTER_DEBUGF(3, "dim (x, y) = %d, %d", _dim[0], _dim[1]);
 
        /* user-defined skew */
        if (NULL != skew_x) {
                _skew[0] = *skew_x;
 
                /*
                        negative scale-x affects skew
                        for now, force skew to be in left-right, top-down orientation
                */
                if (
                        NULL != scale_x &&
                        *scale_x < 0.
                ) {
                        _skew[0] *= -1;
                }
        }
        if (NULL != skew_y) {
                _skew[1] = *skew_y;
 
                /*
                        positive scale-y affects skew
                        for now, force skew to be in left-right, top-down orientation
                */
                if (
                        NULL != scale_y &&
                        *scale_y > 0.
                ) {
                        _skew[1] *= -1;
                }
        }
 
        /*
                if geometry is a point, a linestring or set of either and bounds not set,
                increase extent by a pixel to avoid missing points on border
 
                a whole pixel is used instead of half-pixel due to backward
                compatibility with GDAL 1.6, 1.7 and 1.8.  1.9+ works fine with half-pixel.
        */
        wkbtype = wkbFlatten(OGR_G_GetGeometryType(src_geom));
        if ((
                        (wkbtype == wkbPoint) ||
                        (wkbtype == wkbMultiPoint) ||
                        (wkbtype == wkbLineString) ||
                        (wkbtype == wkbMultiLineString)
                ) &&
                _dim[0] == 0 &&
                _dim[1] == 0
        ) {
 
#if POSTGIS_GDAL_VERSION > 18
 
                RASTER_DEBUG(3, "Adjusting extent for GDAL > 1.8 by half the scale on X-axis");
                extent.MinX -= (_scale[0] / 2.);
                extent.MaxX += (_scale[0] / 2.);
 
                RASTER_DEBUG(3, "Adjusting extent for GDAL > 1.8 by half the scale on Y-axis");
                extent.MinY -= (_scale[1] / 2.);
                extent.MaxY += (_scale[1] / 2.);
 
#else
 
                RASTER_DEBUG(3, "Adjusting extent for GDAL <= 1.8 by the scale on X-axis");
                extent.MinX -= _scale[0];
                extent.MaxX += _scale[0];
 
                RASTER_DEBUG(3, "Adjusting extent for GDAL <= 1.8 by the scale on Y-axis");
                extent.MinY -= _scale[1];
                extent.MaxY += _scale[1];
 
#endif
 
 
                RASTER_DEBUGF(3, "Adjusted extent: %f, %f, %f, %f",
                        extent.MinX, extent.MinY, extent.MaxX, extent.MaxY);
 
                extent.UpperLeftX = extent.MinX;
                extent.UpperLeftY = extent.MaxY;
        }
 
        /* reprocess extent if skewed */
        if (FLT_NEQ(_skew[0], 0.0) || FLT_NEQ(_skew[1], 0.0))
        {
                rt_raster skewedrast;
 
                RASTER_DEBUG(3, "Computing skewed extent's envelope");
 
                skewedrast = rt_raster_compute_skewed_raster(
                        extent,
                        _skew,
                        _scale,
                        0.01
                );
                if (skewedrast == NULL) {
                        rterror("rt_raster_gdal_rasterize: Could not compute skewed raster");
 
                        OGR_G_DestroyGeometry(src_geom);
                        _rti_rasterize_arg_destroy(arg);
                        /* OGRCleanupAll(); */
 
                        return NULL;
                }
 
                _dim[0] = skewedrast->width;
                _dim[1] = skewedrast->height;
 
                extent.UpperLeftX = skewedrast->ipX;
                extent.UpperLeftY = skewedrast->ipY;
 
                rt_raster_destroy(skewedrast);
        }
 
        /* raster dimensions */
        if (!_dim[0])
                _dim[0] = (int) fmax((fabs(extent.MaxX - extent.MinX) + (_scale[0] / 2.)) / _scale[0], 1);
        if (!_dim[1])
                _dim[1] = (int) fmax((fabs(extent.MaxY - extent.MinY) + (_scale[1] / 2.)) / _scale[1], 1);
 
        /* temporary raster */
        rast = rt_raster_new(_dim[0], _dim[1]);
        if (rast == NULL) {
                rterror("rt_raster_gdal_rasterize: Out of memory allocating temporary raster");
 
                OGR_G_DestroyGeometry(src_geom);
                _rti_rasterize_arg_destroy(arg);
                /* OGRCleanupAll(); */
 
                return NULL;
        }
 
        /* set raster's spatial attributes */
        rt_raster_set_offsets(rast, extent.UpperLeftX, extent.UpperLeftY);
        rt_raster_set_scale(rast, _scale[0], -1 * _scale[1]);
        rt_raster_set_skews(rast, _skew[0], _skew[1]);
 
        rt_raster_get_geotransform_matrix(rast, _gt);
        RASTER_DEBUGF(3, "Temp raster's geotransform: %f, %f, %f, %f, %f, %f",
                _gt[0], _gt[1], _gt[2], _gt[3], _gt[4], _gt[5]);
        RASTER_DEBUGF(3, "Temp raster's dimensions (width x height): %d x %d",
                _dim[0], _dim[1]);
 
        /* user-specified upper-left corner */
        if (
                NULL != ul_xw &&
                NULL != ul_yw
        ) {
                ul_user = 1;
 
                RASTER_DEBUGF(4, "Using user-specified upper-left corner: %f, %f", *ul_xw, *ul_yw);
 
                /* set upper-left corner */
                rt_raster_set_offsets(rast, *ul_xw, *ul_yw);
                extent.UpperLeftX = *ul_xw;
                extent.UpperLeftY = *ul_yw;
        }
        else if (
                ((NULL != ul_xw) && (NULL == ul_yw)) ||
                ((NULL == ul_xw) && (NULL != ul_yw))
        ) {
                rterror("rt_raster_gdal_rasterize: Both X and Y upper-left corner values must be provided");
 
                rt_raster_destroy(rast);
                OGR_G_DestroyGeometry(src_geom);
                _rti_rasterize_arg_destroy(arg);
                /* OGRCleanupAll(); */
 
                return NULL;
        }
 
        /* alignment only considered if upper-left corner not provided */
        if (
                !ul_user && (
                        (NULL != grid_xw) || (NULL != grid_yw)
                )
        ) {
 
                if (
                        ((NULL != grid_xw) && (NULL == grid_yw)) ||
                        ((NULL == grid_xw) && (NULL != grid_yw))
                ) {
                        rterror("rt_raster_gdal_rasterize: Both X and Y alignment values must be provided");
 
                        rt_raster_destroy(rast);
                        OGR_G_DestroyGeometry(src_geom);
                        _rti_rasterize_arg_destroy(arg);
                        /* OGRCleanupAll(); */
 
                        return NULL;
                }
 
                RASTER_DEBUGF(4, "Aligning extent to user-specified grid: %f, %f", *grid_xw, *grid_yw);
 
                do {
                        double _r[2] = {0};
                        double _w[2] = {0};
 
                        /* raster is already aligned */
                        if (DBL_EQ(*grid_xw, extent.UpperLeftX) && DBL_EQ(*grid_yw, extent.UpperLeftY)) {
                                RASTER_DEBUG(3, "Skipping raster alignment as it is already aligned to grid");
                                break;
                        }
 
                        extent.UpperLeftX = rast->ipX;
                        extent.UpperLeftY = rast->ipY;
                        rt_raster_set_offsets(rast, *grid_xw, *grid_yw);
 
                        /* process upper-left corner */
                        if (rt_raster_geopoint_to_cell(
                                rast,
                                extent.UpperLeftX, extent.UpperLeftY,
                                &(_r[0]), &(_r[1]),
                                NULL
                        ) != ES_NONE) {
                                rterror("rt_raster_gdal_rasterize: Could not compute raster pixel for spatial coordinates");
 
                                rt_raster_destroy(rast);
                                OGR_G_DestroyGeometry(src_geom);
                                _rti_rasterize_arg_destroy(arg);
                                /* OGRCleanupAll(); */
 
                                return NULL;
                        }
 
                        if (rt_raster_cell_to_geopoint(
                                rast,
                                _r[0], _r[1],
                                &(_w[0]), &(_w[1]),
                                NULL
                        ) != ES_NONE) {
                                rterror("rt_raster_gdal_rasterize: Could not compute spatial coordinates for raster pixel");
 
                                rt_raster_destroy(rast);
                                OGR_G_DestroyGeometry(src_geom);
                                _rti_rasterize_arg_destroy(arg);
                                /* OGRCleanupAll(); */
 
                                return NULL;
                        }
 
                        /* shift occurred */
                        if (FLT_NEQ(_w[0], extent.UpperLeftX)) {
                                if (NULL == width)
                                        rast->width++;
                                else if (NULL == scale_x) {
                                        double _c[2] = {0};
 
                                        rt_raster_set_offsets(rast, extent.UpperLeftX, extent.UpperLeftY);
 
                                        /* get upper-right corner */
                                        if (rt_raster_cell_to_geopoint(
                                                rast,
                                                rast->width, 0,
                                                &(_c[0]), &(_c[1]),
                                                NULL
                                        ) != ES_NONE) {
                                                rterror("rt_raster_gdal_rasterize: Could not compute spatial coordinates for raster pixel");
 
                                                rt_raster_destroy(rast);
                                                OGR_G_DestroyGeometry(src_geom);
                                                _rti_rasterize_arg_destroy(arg);
                                                /* OGRCleanupAll(); */
 
                                                return NULL;
                                        }
 
                                        rast->scaleX = fabs((_c[0] - _w[0]) / ((double) rast->width));
                                }
                        }
                        if (FLT_NEQ(_w[1], extent.UpperLeftY)) {
                                if (NULL == height)
                                        rast->height++;
                                else if (NULL == scale_y) {
                                        double _c[2] = {0};
 
                                        rt_raster_set_offsets(rast, extent.UpperLeftX, extent.UpperLeftY);
 
                                        /* get upper-right corner */
                                        if (rt_raster_cell_to_geopoint(
                                                rast,
                                                0, rast->height,
                                                &(_c[0]), &(_c[1]),
                                                NULL
                                        ) != ES_NONE) {
                                                rterror("rt_raster_gdal_rasterize: Could not compute spatial coordinates for raster pixel");
 
                                                rt_raster_destroy(rast);
                                                OGR_G_DestroyGeometry(src_geom);
                                                _rti_rasterize_arg_destroy(arg);
                                                /* OGRCleanupAll(); */
 
                                                return NULL;
                                        }
 
                                        rast->scaleY = -1 * fabs((_c[1] - _w[1]) / ((double) rast->height));
                                }
                        }
 
                        rt_raster_set_offsets(rast, _w[0], _w[1]);
                }
                while (0);
        }
 
        /*
                after this point, rt_envelope extent is no longer used
        */
 
        /* get key attributes from rast */
        _dim[0] = rast->width;
        _dim[1] = rast->height;
        rt_raster_get_geotransform_matrix(rast, _gt);
 
        /* scale-x is negative or scale-y is positive */
        if ((
                (NULL != scale_x) && (*scale_x < 0.)
        ) || (
                (NULL != scale_y) && (*scale_y > 0)
        )) {
                double _w[2] = {0};
 
                /* negative scale-x */
                if (
                        (NULL != scale_x) &&
                        (*scale_x < 0.)
                ) {
                        RASTER_DEBUG(3, "Processing negative scale-x");
 
                        if (rt_raster_cell_to_geopoint(
                                rast,
                                _dim[0], 0,
                                &(_w[0]), &(_w[1]),
                                NULL
                        ) != ES_NONE) {
                                rterror("rt_raster_gdal_rasterize: Could not compute spatial coordinates for raster pixel");
 
                                rt_raster_destroy(rast);
                                OGR_G_DestroyGeometry(src_geom);
                                _rti_rasterize_arg_destroy(arg);
                                /* OGRCleanupAll(); */
 
                                return NULL;
                        }
 
                        _gt[0] = _w[0];
                        _gt[1] = *scale_x;
 
                        /* check for skew */
                        if (NULL != skew_x && FLT_NEQ(*skew_x, 0.0))
                                _gt[2] = *skew_x;
                }
                /* positive scale-y */
                if (
                        (NULL != scale_y) &&
                        (*scale_y > 0)
                ) {
                        RASTER_DEBUG(3, "Processing positive scale-y");
 
                        if (rt_raster_cell_to_geopoint(
                                rast,
                                0, _dim[1],
                                &(_w[0]), &(_w[1]),
                                NULL
                        ) != ES_NONE) {
                                rterror("rt_raster_gdal_rasterize: Could not compute spatial coordinates for raster pixel");
 
                                rt_raster_destroy(rast);
                                OGR_G_DestroyGeometry(src_geom);
                                _rti_rasterize_arg_destroy(arg);
                                /* OGRCleanupAll(); */
 
                                return NULL;
                        }
 
                        _gt[3] = _w[1];
                        _gt[5] = *scale_y;
 
                        /* check for skew */
                        if (NULL != skew_y && FLT_NEQ(*skew_y, 0.0))
                                _gt[4] = *skew_y;
                }
        }
 
        rt_raster_destroy(rast);
        rast = NULL;
 
        RASTER_DEBUGF(3, "Applied geotransform: %f, %f, %f, %f, %f, %f",
                _gt[0], _gt[1], _gt[2], _gt[3], _gt[4], _gt[5]);
        RASTER_DEBUGF(3, "Raster dimensions (width x height): %d x %d",
                _dim[0], _dim[1]);
 
        /* load GDAL mem */
        if (!rt_util_gdal_driver_registered("MEM")) {
                RASTER_DEBUG(4, "Registering MEM driver");
                GDALRegister_MEM();
                unload_drv = 1;
        }
        _drv = GDALGetDriverByName("MEM");
        if (NULL == _drv) {
                rterror("rt_raster_gdal_rasterize: Could not load the MEM GDAL driver");
 
                OGR_G_DestroyGeometry(src_geom);
                _rti_rasterize_arg_destroy(arg);
                /* OGRCleanupAll(); */
 
                return NULL;
        }
 
        /* unload driver from GDAL driver manager */
        if (unload_drv) {
                RASTER_DEBUG(4, "Deregistering MEM driver");
                GDALDeregisterDriver(_drv);
        }
 
        _ds = GDALCreate(_drv, "", _dim[0], _dim[1], 0, GDT_Byte, NULL);
        if (NULL == _ds) {
                rterror("rt_raster_gdal_rasterize: Could not create a GDALDataset to rasterize the geometry into");
 
                OGR_G_DestroyGeometry(src_geom);
                _rti_rasterize_arg_destroy(arg);
                /* OGRCleanupAll(); */
                if (unload_drv) GDALDestroyDriver(_drv);
 
                return NULL;
        }
 
        /* set geotransform */
        cplerr = GDALSetGeoTransform(_ds, _gt);
        if (cplerr != CE_None) {
                rterror("rt_raster_gdal_rasterize: Could not set geotransform on GDALDataset");
 
                OGR_G_DestroyGeometry(src_geom);
                _rti_rasterize_arg_destroy(arg);
                /* OGRCleanupAll(); */
 
                GDALClose(_ds);
                if (unload_drv) GDALDestroyDriver(_drv);
 
                return NULL;
        }
 
        /* set SRS */
        if (NULL != arg->src_sr) {
                char *_srs = NULL;
                OSRExportToWkt(arg->src_sr, &_srs);
 
                cplerr = GDALSetProjection(_ds, _srs);
                CPLFree(_srs);
                if (cplerr != CE_None) {
                        rterror("rt_raster_gdal_rasterize: Could not set projection on GDALDataset");
 
                        OGR_G_DestroyGeometry(src_geom);
                        _rti_rasterize_arg_destroy(arg);
                        /* OGRCleanupAll(); */
 
                        GDALClose(_ds);
                if (unload_drv) GDALDestroyDriver(_drv);
 
                        return NULL;
                }
        }
 
        /* set bands */
        for (i = 0; i < arg->numbands; i++) {
                err = 0;
 
                do {
                        /* add band */
                        cplerr = GDALAddBand(_ds, rt_util_pixtype_to_gdal_datatype(arg->pixtype[i]), NULL);
                        if (cplerr != CE_None) {
                                rterror("rt_raster_gdal_rasterize: Could not add band to GDALDataset");
                                err = 1;
                                break;
                        }
 
                        _band = GDALGetRasterBand(_ds, i + 1);
                        if (NULL == _band) {
                                rterror("rt_raster_gdal_rasterize: Could not get band %d from GDALDataset", i + 1);
                                err = 1;
                                break;
                        }
 
                        /* nodata value */
                        if (arg->hasnodata[i]) {
                                RASTER_DEBUGF(4, "Setting NODATA value of band %d to %f", i, arg->nodata[i]);
                                cplerr = GDALSetRasterNoDataValue(_band, arg->nodata[i]);
                                if (cplerr != CE_None) {
                                        rterror("rt_raster_gdal_rasterize: Could not set nodata value");
                                        err = 1;
                                        break;
                                }
                                RASTER_DEBUGF(4, "NODATA value set to %f", GDALGetRasterNoDataValue(_band, NULL));
                        }
 
                        /* initial value */
                        cplerr = GDALFillRaster(_band, arg->init[i], 0);
                        if (cplerr != CE_None) {
                                rterror("rt_raster_gdal_rasterize: Could not set initial value");
                                err = 1;
                                break;
                        }
                }
                while (0);
 
                if (err) {
 
                        OGR_G_DestroyGeometry(src_geom);
                        _rti_rasterize_arg_destroy(arg);
 
                        /* OGRCleanupAll(); */
 
                        GDALClose(_ds);
                        if (unload_drv) GDALDestroyDriver(_drv);
 
                        return NULL;
                }
        }
 
        arg->bandlist = (int *) rtalloc(sizeof(int) * arg->numbands);
        for (i = 0; i < arg->numbands; i++) arg->bandlist[i] = i + 1;
 
        /* burn geometry */
        cplerr = GDALRasterizeGeometries(
                _ds,
                arg->numbands, arg->bandlist,
                1, &src_geom,
                NULL, NULL,
                arg->value,
                options,
                NULL, NULL
        );
        if (cplerr != CE_None) {
                rterror("rt_raster_gdal_rasterize: Could not rasterize geometry");
 
                OGR_G_DestroyGeometry(src_geom);
                _rti_rasterize_arg_destroy(arg);
                /* OGRCleanupAll(); */
 
                GDALClose(_ds);
                if (unload_drv) GDALDestroyDriver(_drv);
 
                return NULL;
        }
 
        /* convert gdal dataset to raster */
        GDALFlushCache(_ds);
        RASTER_DEBUG(3, "Converting GDAL dataset to raster");
        rast = rt_raster_from_gdal_dataset(_ds);
 
        OGR_G_DestroyGeometry(src_geom);
        /* OGRCleanupAll(); */
 
        GDALClose(_ds);
        if (unload_drv) GDALDestroyDriver(_drv);
 
        if (NULL == rast) {
                rterror("rt_raster_gdal_rasterize: Could not rasterize geometry");
                return NULL;
        }
 
        /* width, height */
        _width = rt_raster_get_width(rast);
        _height = rt_raster_get_height(rast);
 
        /* check each band for pixtype */
        for (i = 0; i < arg->numbands; i++) {
                uint8_t *data = NULL;
                rt_band band = NULL;
                rt_band oldband = NULL;
 
                double val = 0;
                int nodata = 0;
                int hasnodata = 0;
                double nodataval = 0;
                int x = 0;
                int y = 0;
 
                oldband = rt_raster_get_band(rast, i);
                if (oldband == NULL) {
                        rterror("rt_raster_gdal_rasterize: Could not get band %d of output raster", i);
                        _rti_rasterize_arg_destroy(arg);
                        rt_raster_destroy(rast);
                        return NULL;
                }
 
                /* band is of user-specified type */
                if (rt_band_get_pixtype(oldband) == arg->pixtype[i])
                        continue;
 
                /* hasnodata, nodataval */
                hasnodata = rt_band_get_hasnodata_flag(oldband);
                if (hasnodata)
                        rt_band_get_nodata(oldband, &nodataval);
 
                /* allocate data */
                data = rtalloc(rt_pixtype_size(arg->pixtype[i]) * _width * _height);
                if (data == NULL) {
                        rterror("rt_raster_gdal_rasterize: Could not allocate memory for band data");
                        _rti_rasterize_arg_destroy(arg);
                        rt_raster_destroy(rast);
                        return NULL;
                }
                memset(data, 0, rt_pixtype_size(arg->pixtype[i]) * _width * _height);
 
                /* create new band of correct type */
                band = rt_band_new_inline(
                        _width, _height,
                        arg->pixtype[i],
                        hasnodata, nodataval,
                        data
                );
                if (band == NULL) {
                        rterror("rt_raster_gdal_rasterize: Could not create band");
                        rtdealloc(data);
                        _rti_rasterize_arg_destroy(arg);
                        rt_raster_destroy(rast);
                        return NULL;
                }
 
                /* give ownership of data to band */
                rt_band_set_ownsdata_flag(band, 1);
 
                /* copy pixel by pixel */
                for (x = 0; x < _width; x++) {
                        for (y = 0; y < _height; y++) {
                                err = rt_band_get_pixel(oldband, x, y, &val, &nodata);
                                if (err != ES_NONE) {
                                        rterror("rt_raster_gdal_rasterize: Could not get pixel value");
                                        _rti_rasterize_arg_destroy(arg);
                                        rt_raster_destroy(rast);
                                        rt_band_destroy(band);
                                        return NULL;
                                }
 
                                if (nodata)
                                        val = nodataval;
 
                                err = rt_band_set_pixel(band, x, y, val, NULL);
                                if (err != ES_NONE) {
                                        rterror("rt_raster_gdal_rasterize: Could not set pixel value");
                                        _rti_rasterize_arg_destroy(arg);
                                        rt_raster_destroy(rast);
                                        rt_band_destroy(band);
                                        return NULL;
                                }
                        }
                }
 
                /* replace band */
                oldband = rt_raster_replace_band(rast, band, i);
                if (oldband == NULL) {
                        rterror("rt_raster_gdal_rasterize: Could not replace band %d of output raster", i);
                        _rti_rasterize_arg_destroy(arg);
                        rt_raster_destroy(rast);
                        rt_band_destroy(band);
                        return NULL;
                }
 
                /* free oldband */
                rt_band_destroy(oldband);
        }
 
        _rti_rasterize_arg_destroy(arg);
 
        RASTER_DEBUG(3, "done");
 
        return rast;
}
 
/******************************************************************************
* rt_raster_from_two_rasters()
******************************************************************************/
 
/*
 * Return raster of computed extent specified extenttype applied
 * on two input rasters.  The raster returned should be freed by
 * the caller
 *
 * @param rast1 : the first raster
 * @param rast2 : the second raster
 * @param extenttype : type of extent for the output raster
 * @param *rtnraster : raster of computed extent
 * @param *offset : 4-element array indicating the X,Y offsets
 * for each raster. 0,1 for rast1 X,Y. 2,3 for rast2 X,Y.
 *
 * @return ES_NONE if success, ES_ERROR if error
 */
rt_errorstate
rt_raster_from_two_rasters(
        rt_raster rast1, rt_raster rast2,
        rt_extenttype extenttype,
        rt_raster *rtnraster, double *offset
) {
        int i;
 
        rt_raster _rast[2] = {rast1, rast2};
        double _offset[2][4] = {{0.}};
        uint16_t _dim[2][2] = {{0}};
 
        rt_raster raster = NULL;
        int aligned = 0;
        int dim[2] = {0};
        double gt[6] = {0.};
 
        assert(NULL != rast1);
        assert(NULL != rast2);
        assert(NULL != rtnraster);
 
        /* set rtnraster to NULL */
        *rtnraster = NULL;
 
        /* rasters must be aligned */
        if (rt_raster_same_alignment(rast1, rast2, &aligned, NULL) != ES_NONE) {
                rterror("rt_raster_from_two_rasters: Could not test for alignment on the two rasters");
                return ES_ERROR;
        }
        if (!aligned) {
                rterror("rt_raster_from_two_rasters: The two rasters provided do not have the same alignment");
                return ES_ERROR;
        }
 
        /* dimensions */
        _dim[0][0] = rast1->width;
        _dim[0][1] = rast1->height;
        _dim[1][0] = rast2->width;
        _dim[1][1] = rast2->height;
 
        /* get raster offsets */
        if (rt_raster_geopoint_to_cell(
                _rast[1],
                _rast[0]->ipX, _rast[0]->ipY,
                &(_offset[1][0]), &(_offset[1][1]),
                NULL
        ) != ES_NONE) {
                rterror("rt_raster_from_two_rasters: Could not compute offsets of the second raster relative to the first raster");
                return ES_ERROR;
        }
        _offset[1][0] = -1 * _offset[1][0];
        _offset[1][1] = -1 * _offset[1][1];
        _offset[1][2] = _offset[1][0] + _dim[1][0] - 1;
        _offset[1][3] = _offset[1][1] + _dim[1][1] - 1;
 
        i = -1;
        switch (extenttype) {
                case ET_FIRST:
                        i = 0;
                        _offset[0][0] = 0.;
                        _offset[0][1] = 0.;
                        /* FALLTHROUGH */
                case ET_LAST:
                case ET_SECOND:
                        if (i < 0) {
                                i = 1;
                                _offset[0][0] = -1 * _offset[1][0];
                                _offset[0][1] = -1 * _offset[1][1];
                                _offset[1][0] = 0.;
                                _offset[1][1] = 0.;
                        }
 
                        dim[0] = _dim[i][0];
                        dim[1] = _dim[i][1];
                        raster = rt_raster_new(
                                dim[0],
                                dim[1]
                        );
                        if (raster == NULL) {
                                rterror("rt_raster_from_two_rasters: Could not create output raster");
                                return ES_ERROR;
                        }
                        rt_raster_set_srid(raster, _rast[i]->srid);
                        rt_raster_get_geotransform_matrix(_rast[i], gt);
                        rt_raster_set_geotransform_matrix(raster, gt);
                        break;
                case ET_UNION: {
                        double off[4] = {0};
 
                        rt_raster_get_geotransform_matrix(_rast[0], gt);
                        RASTER_DEBUGF(4, "gt = (%f, %f, %f, %f, %f, %f)",
                                gt[0],
                                gt[1],
                                gt[2],
                                gt[3],
                                gt[4],
                                gt[5]
                        );
 
                        /* new raster upper left offset */
                        off[0] = 0;
                        if (_offset[1][0] < 0)
                                off[0] = _offset[1][0];
                        off[1] = 0;
                        if (_offset[1][1] < 0)
                                off[1] = _offset[1][1];
 
                        /* new raster lower right offset */
                        off[2] = _dim[0][0] - 1;
                        if ((int) _offset[1][2] >= _dim[0][0])
                                off[2] = _offset[1][2];
                        off[3] = _dim[0][1] - 1;
                        if ((int) _offset[1][3] >= _dim[0][1])
                                off[3] = _offset[1][3];
 
                        /* upper left corner */
                        if (rt_raster_cell_to_geopoint(
                                _rast[0],
                                off[0], off[1],
                                &(gt[0]), &(gt[3]),
                                NULL
                        ) != ES_NONE) {
                                rterror("rt_raster_from_two_rasters: Could not get spatial coordinates of upper-left pixel of output raster");
                                return ES_ERROR;
                        }
 
                        dim[0] = off[2] - off[0] + 1;
                        dim[1] = off[3] - off[1] + 1;
                        RASTER_DEBUGF(4, "off = (%f, %f, %f, %f)",
                                off[0],
                                off[1],
                                off[2],
                                off[3]
                        );
                        RASTER_DEBUGF(4, "dim = (%d, %d)", dim[0], dim[1]);
 
                        raster = rt_raster_new(
                                dim[0],
                                dim[1]
                        );
                        if (raster == NULL) {
                                rterror("rt_raster_from_two_rasters: Could not create output raster");
                                return ES_ERROR;
                        }
                        rt_raster_set_srid(raster, _rast[0]->srid);
                        rt_raster_set_geotransform_matrix(raster, gt);
                        RASTER_DEBUGF(4, "gt = (%f, %f, %f, %f, %f, %f)",
                                gt[0],
                                gt[1],
                                gt[2],
                                gt[3],
                                gt[4],
                                gt[5]
                        );
 
                        /* get offsets */
                        if (rt_raster_geopoint_to_cell(
                                _rast[0],
                                gt[0], gt[3],
                                &(_offset[0][0]), &(_offset[0][1]),
                                NULL
                        ) != ES_NONE) {
                                rterror("rt_raster_from_two_rasters: Could not get offsets of the FIRST raster relative to the output raster");
                                rt_raster_destroy(raster);
                                return ES_ERROR;
                        }
                        _offset[0][0] *= -1;
                        _offset[0][1] *= -1;
 
                        if (rt_raster_geopoint_to_cell(
                                _rast[1],
                                gt[0], gt[3],
                                &(_offset[1][0]), &(_offset[1][1]),
                                NULL
                        ) != ES_NONE) {
                                rterror("rt_raster_from_two_rasters: Could not get offsets of the SECOND raster relative to the output raster");
                                rt_raster_destroy(raster);
                                return ES_ERROR;
                        }
                        _offset[1][0] *= -1;
                        _offset[1][1] *= -1;
                        break;
                }
                case ET_INTERSECTION: {
                        double off[4] = {0};
 
                        /* no intersection */
                        if (
                                (_offset[1][2] < 0 || _offset[1][0] > (_dim[0][0] - 1)) ||
                                (_offset[1][3] < 0 || _offset[1][1] > (_dim[0][1] - 1))
                        ) {
                                RASTER_DEBUG(3, "The two rasters provided have no intersection.  Returning no band raster");
 
                                raster = rt_raster_new(0, 0);
                                if (raster == NULL) {
                                        rterror("rt_raster_from_two_rasters: Could not create output raster");
                                        return ES_ERROR;
                                }
                                rt_raster_set_srid(raster, _rast[0]->srid);
                                rt_raster_set_scale(raster, 0, 0);
 
                                /* set offsets if provided */
                                if (NULL != offset) {
                                        for (i = 0; i < 4; i++)
                                                offset[i] = _offset[i / 2][i % 2];
                                }
 
                                *rtnraster = raster;
                                return ES_NONE;
                        }
 
                        if (_offset[1][0] > 0)
                                off[0] = _offset[1][0];
                        if (_offset[1][1] > 0)
                                off[1] = _offset[1][1];
 
                        off[2] = _dim[0][0] - 1;
                        if (_offset[1][2] < _dim[0][0])
                                off[2] = _offset[1][2];
                        off[3] = _dim[0][1] - 1;
                        if (_offset[1][3] < _dim[0][1])
                                off[3] = _offset[1][3];
 
                        dim[0] = off[2] - off[0] + 1;
                        dim[1] = off[3] - off[1] + 1;
                        raster = rt_raster_new(
                                dim[0],
                                dim[1]
                        );
                        if (raster == NULL) {
                                rterror("rt_raster_from_two_rasters: Could not create output raster");
                                return ES_ERROR;
                        }
                        rt_raster_set_srid(raster, _rast[0]->srid);
 
                        /* get upper-left corner */
                        rt_raster_get_geotransform_matrix(_rast[0], gt);
                        if (rt_raster_cell_to_geopoint(
                                _rast[0],
                                off[0], off[1],
                                &(gt[0]), &(gt[3]),
                                gt
                        ) != ES_NONE) {
                                rterror("rt_raster_from_two_rasters: Could not get spatial coordinates of upper-left pixel of output raster");
                                rt_raster_destroy(raster);
                                return ES_ERROR;
                        }
 
                        rt_raster_set_geotransform_matrix(raster, gt);
 
                        /* get offsets */
                        if (rt_raster_geopoint_to_cell(
                                _rast[0],
                                gt[0], gt[3],
                                &(_offset[0][0]), &(_offset[0][1]),
                                NULL
                        ) != ES_NONE) {
                                rterror("rt_raster_from_two_rasters: Could not get pixel coordinates to compute the offsets of the FIRST raster relative to the output raster");
                                rt_raster_destroy(raster);
                                return ES_ERROR;
                        }
                        _offset[0][0] *= -1;
                        _offset[0][1] *= -1;
 
                        if (rt_raster_geopoint_to_cell(
                                _rast[1],
                                gt[0], gt[3],
                                &(_offset[1][0]), &(_offset[1][1]),
                                NULL
                        ) != ES_NONE) {
                                rterror("rt_raster_from_two_rasters: Could not get pixel coordinates to compute the offsets of the SECOND raster relative to the output raster");
                                rt_raster_destroy(raster);
                                return ES_ERROR;
                        }
                        _offset[1][0] *= -1;
                        _offset[1][1] *= -1;
                        break;
                }
                case ET_CUSTOM:
                        rterror("rt_raster_from_two_rasters: Extent type ET_CUSTOM is not supported by this function");
                        break;
        }
 
        /* set offsets if provided */
        if (NULL != offset) {
                for (i = 0; i < 4; i++)
                        offset[i] = _offset[i / 2][i % 2];
        }
 
        *rtnraster = raster;
        return ES_NONE;
}