rt_raster_from_gdal_dataset()

rt_raster rt_raster_from_gdal_dataset

(

GDALDatasetH

ds

)

Return a raster from a GDAL dataset.

Parameters

ds	: the GDAL dataset to convert to a raster

Returns

raster or NULL

Definition at line 2291 of file rt_raster.c.

                                             {
        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;
}

References ovdump::band, pixval::ds, ES_NONE, FALSE, window::gt, PT_END, rtpixdump::rast, RASTER_DEBUG, RASTER_DEBUGF, rt_band_set_pixel_line(), rt_pixtype_size(), rt_raster_destroy(), rt_raster_generate_new_band(), rt_raster_get_band(), rt_raster_get_num_bands(), rt_raster_new(), rt_raster_set_geotransform_matrix(), rt_raster_set_srid(), rt_util_gdal_datatype_to_pixtype(), rt_util_gdal_sr_auth_info(), rtalloc(), rtdealloc(), rterror(), pixval::x, and pixval::y.

Referenced by RASTER_fromGDALRaster(), rt_band_load_offline_data(), rt_raster_gdal_rasterize(), rt_raster_gdal_warp(), and test_gdal_to_raster().

Here is the call graph for this function:

Here is the caller graph for this function: