rt_raster_gdal_rasterize()

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
	)

Return a raster of the provided geometry.

Parameters

wkb	: WKB representation of the geometry to convert
wkb_len	: length of the WKB representation of the geometry
srs	: the geometry's coordinate system in OGC WKT
num_bands	: number of bands in the output raster
pixtype	: data type of each band
init	: array of values to initialize each band with
value	: array of values for pixels of geometry
nodata	: array of nodata values for each band
hasnodata	: array flagging the presence of nodata for each band
width	: the number of columns of the raster
height	: the number of rows of the raster
scale_x	: the pixel width of the raster
scale_y	: the pixel height of the raster
ul_xw	: the X value of upper-left corner of the raster
ul_yw	: the Y value of upper-left corner of the raster
grid_xw	: the X value of point on grid to align raster to
grid_yw	: the Y value of point on grid to align raster to
skew_x	: the X skew of the raster
skew_y	: the Y skew of the raster
options	: array of options. only option is "ALL_TOUCHED"

Returns

the raster of the provided geometry or NULL

Definition at line 2502 of file rt_raster.c.

  {
        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 > 10800
 
                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((size_t)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, (size_t)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;
}

References _rti_rasterize_arg_destroy(), _rti_rasterize_arg_init(), _rti_rasterize_arg_t::bandlist, DBL_EQ, ES_NONE, FLT_NEQ, _rti_rasterize_arg_t::hasnodata, rt_raster_t::height, _rti_rasterize_arg_t::init, rt_raster_t::ipX, rt_raster_t::ipY, rt_envelope::MaxX, rt_envelope::MaxY, rt_envelope::MinX, rt_envelope::MinY, _rti_rasterize_arg_t::noband, _rti_rasterize_arg_t::nodata, _rti_rasterize_arg_t::numbands, _rti_rasterize_arg_t::pixtype, PT_8BUI, RASTER_DEBUG, RASTER_DEBUGF, rt_band_destroy(), rt_band_get_hasnodata_flag(), rt_band_get_nodata(), rt_band_get_pixel(), rt_band_get_pixtype(), rt_band_new_inline(), rt_band_set_ownsdata_flag(), rt_band_set_pixel(), rt_pixtype_size(), rt_raster_cell_to_geopoint(), rt_raster_compute_skewed_raster(), rt_raster_destroy(), rt_raster_from_gdal_dataset(), rt_raster_geopoint_to_cell(), rt_raster_get_band(), rt_raster_get_geotransform_matrix(), rt_raster_get_height(), rt_raster_get_width(), rt_raster_new(), rt_raster_replace_band(), rt_raster_set_offsets(), rt_raster_set_scale(), rt_raster_set_skews(), rt_util_from_ogr_envelope(), rt_util_gdal_driver_registered(), rt_util_pixtype_to_gdal_datatype(), rtalloc(), rtdealloc(), rterror(), rtinfo(), _rti_rasterize_arg_t::src_sr, rt_envelope::UpperLeftX, rt_envelope::UpperLeftY, _rti_rasterize_arg_t::value, and rt_raster_t::width.

Referenced by RASTER_asRaster(), RASTER_clip(), RASTER_setPixelValuesGeomval(), and test_gdal_rasterize().

Here is the call graph for this function:

Here is the caller graph for this function: