PostGIS  2.1.10dev-r@@SVN_REVISION@@
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 
)

Return GDAL dataset using GDAL MEM driver from raster.

Parameters
raster: raster to convert to GDAL MEM
srs: the raster's coordinate system in OGC WKT
bandNums: array of band numbers to extract from raster and include in the GDAL dataset (0 based)
excludeNodataValues: array of zero, nonzero where if non-zero, ignore nodata values for the band
count: number of elements in bandNums
rtn_drv: is set to the GDAL driver object
destroy_rtn_drv: if non-zero, caller must destroy the MEM driver
Returns
GDAL dataset using GDAL MEM driver

Definition at line 8989 of file rt_api.c.

References ovdump::band, genraster::count, pixval::ds, ES_NONE, FALSE, window::gt, rt_raster_serialized_t::height, rt_raster_serialized_t::numBands, PT_16BSI, PT_8BSI, PT_END, RASTER_DEBUG, RASTER_DEBUGF, rt_band_get_data(), rt_band_get_hasnodata_flag(), rt_band_get_nodata(), rt_band_get_pixel(), rt_band_get_pixtype(), rt_pixtype_size(), rt_raster_get_band(), rt_raster_get_geotransform_matrix(), rt_raster_get_height(), rt_raster_get_num_bands(), rt_raster_get_width(), rt_util_clamp_to_16BSI(), rt_util_gdal_convert_sr(), rt_util_gdal_driver_registered(), rt_util_pixtype_to_gdal_datatype(), rtalloc(), rtdealloc(), rterror(), rtwarn(), genraster::value, rt_raster_serialized_t::width, pixval::x, and pixval::y.

Referenced by rt_raster_gdal_polygonize(), rt_raster_gdal_warp(), rt_raster_to_gdal(), and test_gdal_to_raster().

8996  {
8997  GDALDriverH drv = NULL;
8998  GDALDatasetH ds = NULL;
8999  double gt[6] = {0.0};
9000  CPLErr cplerr;
9001  GDALDataType gdal_pt = GDT_Unknown;
9002  GDALRasterBandH band;
9003  void *pVoid;
9004  char *pszDataPointer;
9005  char szGDALOption[50];
9006  char *apszOptions[4];
9007  double nodata = 0.0;
9008  int allocBandNums = 0;
9009  int allocNodataValues = 0;
9010 
9011  int i;
9012  int numBands;
9013  uint32_t width = 0;
9014  uint32_t height = 0;
9015  rt_band rtband = NULL;
9016  rt_pixtype pt = PT_END;
9017 
9018  assert(NULL != raster);
9019  assert(NULL != rtn_drv);
9020  assert(NULL != destroy_rtn_drv);
9021 
9022  /* store raster in GDAL MEM raster */
9023  if (!rt_util_gdal_driver_registered("MEM")) {
9024  RASTER_DEBUG(4, "Registering MEM driver");
9025  GDALRegister_MEM();
9026  *destroy_rtn_drv = 1;
9027  }
9028  drv = GDALGetDriverByName("MEM");
9029  if (NULL == drv) {
9030  rterror("rt_raster_to_gdal_mem: Could not load the MEM GDAL driver");
9031  return 0;
9032  }
9033  *rtn_drv = drv;
9034 
9035  /* unload driver from GDAL driver manager */
9036  if (*destroy_rtn_drv) {
9037  RASTER_DEBUG(4, "Deregistering MEM driver");
9038  GDALDeregisterDriver(drv);
9039  }
9040 
9041  width = rt_raster_get_width(raster);
9042  height = rt_raster_get_height(raster);
9043  ds = GDALCreate(
9044  drv, "",
9045  width, height,
9046  0, GDT_Byte, NULL
9047  );
9048  if (NULL == ds) {
9049  rterror("rt_raster_to_gdal_mem: Could not create a GDALDataset to convert into");
9050  return 0;
9051  }
9052 
9053  /* add geotransform */
9055  cplerr = GDALSetGeoTransform(ds, gt);
9056  if (cplerr != CE_None) {
9057  rterror("rt_raster_to_gdal_mem: Could not set geotransformation");
9058  GDALClose(ds);
9059  return 0;
9060  }
9061 
9062  /* set spatial reference */
9063  if (NULL != srs && strlen(srs)) {
9064  char *_srs = rt_util_gdal_convert_sr(srs, 0);
9065  if (_srs == NULL) {
9066  rterror("rt_raster_to_gdal_mem: Could not convert srs to GDAL accepted format");
9067  GDALClose(ds);
9068  return 0;
9069  }
9070 
9071  cplerr = GDALSetProjection(ds, _srs);
9072  CPLFree(_srs);
9073  if (cplerr != CE_None) {
9074  rterror("rt_raster_to_gdal_mem: Could not set projection");
9075  GDALClose(ds);
9076  return 0;
9077  }
9078  RASTER_DEBUGF(3, "Projection set to: %s", GDALGetProjectionRef(ds));
9079  }
9080 
9081  /* process bandNums */
9082  numBands = rt_raster_get_num_bands(raster);
9083  if (NULL != bandNums && count > 0) {
9084  for (i = 0; i < count; i++) {
9085  if (bandNums[i] >= numBands) {
9086  rterror("rt_raster_to_gdal_mem: The band index %d is invalid", bandNums[i]);
9087  GDALClose(ds);
9088  return 0;
9089  }
9090  }
9091  }
9092  else {
9093  count = numBands;
9094  bandNums = (uint32_t *) rtalloc(sizeof(uint32_t) * count);
9095  if (NULL == bandNums) {
9096  rterror("rt_raster_to_gdal_mem: Could not allocate memory for band indices");
9097  GDALClose(ds);
9098  return 0;
9099  }
9100  allocBandNums = 1;
9101  for (i = 0; i < count; i++) bandNums[i] = i;
9102  }
9103 
9104  /* process exclude_nodata_values */
9105  if (NULL == excludeNodataValues) {
9106  excludeNodataValues = (int *) rtalloc(sizeof(int) * count);
9107  if (NULL == excludeNodataValues) {
9108  rterror("rt_raster_to_gdal_mem: Could not allocate memory for NODATA flags");
9109  GDALClose(ds);
9110  return 0;
9111  }
9112  allocNodataValues = 1;
9113  for (i = 0; i < count; i++) excludeNodataValues[i] = 1;
9114  }
9115 
9116  /* add band(s) */
9117  for (i = 0; i < count; i++) {
9118  rtband = rt_raster_get_band(raster, bandNums[i]);
9119  if (NULL == rtband) {
9120  rterror("rt_raster_to_gdal_mem: Could not get requested band index %d", bandNums[i]);
9121  if (allocBandNums) rtdealloc(bandNums);
9122  if (allocNodataValues) rtdealloc(excludeNodataValues);
9123  GDALClose(ds);
9124  return 0;
9125  }
9126 
9127  pt = rt_band_get_pixtype(rtband);
9128  gdal_pt = rt_util_pixtype_to_gdal_datatype(pt);
9129  if (gdal_pt == GDT_Unknown)
9130  rtwarn("rt_raster_to_gdal_mem: Unknown pixel type for band");
9131 
9132  /*
9133  For all pixel types other than PT_8BSI, set pointer to start of data
9134  */
9135  if (pt != PT_8BSI) {
9136  pVoid = rt_band_get_data(rtband);
9137  RASTER_DEBUGF(4, "Band data is at pos %p", pVoid);
9138 
9139  pszDataPointer = (char *) rtalloc(20 * sizeof (char));
9140  sprintf(pszDataPointer, "%p", pVoid);
9141  RASTER_DEBUGF(4, "rt_raster_to_gdal_mem: szDatapointer is %p",
9142  pszDataPointer);
9143 
9144  if (strnicmp(pszDataPointer, "0x", 2) == 0)
9145  sprintf(szGDALOption, "DATAPOINTER=%s", pszDataPointer);
9146  else
9147  sprintf(szGDALOption, "DATAPOINTER=0x%s", pszDataPointer);
9148 
9149  RASTER_DEBUG(3, "Storing info for GDAL MEM raster band");
9150 
9151  apszOptions[0] = szGDALOption;
9152  apszOptions[1] = NULL; /* pixel offset, not needed */
9153  apszOptions[2] = NULL; /* line offset, not needed */
9154  apszOptions[3] = NULL;
9155 
9156  /* free */
9157  rtdealloc(pszDataPointer);
9158 
9159  /* add band */
9160  if (GDALAddBand(ds, gdal_pt, apszOptions) == CE_Failure) {
9161  rterror("rt_raster_to_gdal_mem: Could not add GDAL raster band");
9162  if (allocBandNums) rtdealloc(bandNums);
9163  GDALClose(ds);
9164  return 0;
9165  }
9166  }
9167  /*
9168  PT_8BSI is special as GDAL has no equivalent pixel type.
9169  Must convert 8BSI to 16BSI so create basic band
9170  */
9171  else {
9172  /* add band */
9173  if (GDALAddBand(ds, gdal_pt, NULL) == CE_Failure) {
9174  rterror("rt_raster_to_gdal_mem: Could not add GDAL raster band");
9175  if (allocBandNums) rtdealloc(bandNums);
9176  if (allocNodataValues) rtdealloc(excludeNodataValues);
9177  GDALClose(ds);
9178  return 0;
9179  }
9180  }
9181 
9182  /* check band count */
9183  if (GDALGetRasterCount(ds) != i + 1) {
9184  rterror("rt_raster_to_gdal_mem: Error creating GDAL MEM raster band");
9185  if (allocBandNums) rtdealloc(bandNums);
9186  if (allocNodataValues) rtdealloc(excludeNodataValues);
9187  GDALClose(ds);
9188  return 0;
9189  }
9190 
9191  /* get new band */
9192  band = NULL;
9193  band = GDALGetRasterBand(ds, i + 1);
9194  if (NULL == band) {
9195  rterror("rt_raster_to_gdal_mem: Could not get GDAL band for additional processing");
9196  if (allocBandNums) rtdealloc(bandNums);
9197  if (allocNodataValues) rtdealloc(excludeNodataValues);
9198  GDALClose(ds);
9199  return 0;
9200  }
9201 
9202  /* PT_8BSI requires manual setting of pixels */
9203  if (pt == PT_8BSI) {
9204  int nXBlocks, nYBlocks;
9205  int nXBlockSize, nYBlockSize;
9206  int iXBlock, iYBlock;
9207  int nXValid, nYValid;
9208  int iX, iY;
9209  int iXMax, iYMax;
9210 
9211  int x, y, z;
9212  uint32_t valueslen = 0;
9213  int16_t *values = NULL;
9214  double value = 0.;
9215 
9216  /* this makes use of GDAL's "natural" blocks */
9217  GDALGetBlockSize(band, &nXBlockSize, &nYBlockSize);
9218  nXBlocks = (width + nXBlockSize - 1) / nXBlockSize;
9219  nYBlocks = (height + nYBlockSize - 1) / nYBlockSize;
9220  RASTER_DEBUGF(4, "(nXBlockSize, nYBlockSize) = (%d, %d)", nXBlockSize, nYBlockSize);
9221  RASTER_DEBUGF(4, "(nXBlocks, nYBlocks) = (%d, %d)", nXBlocks, nYBlocks);
9222 
9223  /* length is for the desired pixel type */
9224  valueslen = rt_pixtype_size(PT_16BSI) * nXBlockSize * nYBlockSize;
9225  values = rtalloc(valueslen);
9226  if (NULL == values) {
9227  rterror("rt_raster_to_gdal_mem: Could not allocate memory for GDAL band pixel values");
9228  if (allocBandNums) rtdealloc(bandNums);
9229  if (allocNodataValues) rtdealloc(excludeNodataValues);
9230  GDALClose(ds);
9231  return 0;
9232  }
9233 
9234  for (iYBlock = 0; iYBlock < nYBlocks; iYBlock++) {
9235  for (iXBlock = 0; iXBlock < nXBlocks; iXBlock++) {
9236  memset(values, 0, valueslen);
9237 
9238  x = iXBlock * nXBlockSize;
9239  y = iYBlock * nYBlockSize;
9240  RASTER_DEBUGF(4, "(iXBlock, iYBlock) = (%d, %d)", iXBlock, iYBlock);
9241  RASTER_DEBUGF(4, "(x, y) = (%d, %d)", x, y);
9242 
9243  /* valid block width */
9244  if ((iXBlock + 1) * nXBlockSize > width)
9245  nXValid = width - (iXBlock * nXBlockSize);
9246  else
9247  nXValid = nXBlockSize;
9248 
9249  /* valid block height */
9250  if ((iYBlock + 1) * nYBlockSize > height)
9251  nYValid = height - (iYBlock * nYBlockSize);
9252  else
9253  nYValid = nYBlockSize;
9254 
9255  RASTER_DEBUGF(4, "(nXValid, nYValid) = (%d, %d)", nXValid, nYValid);
9256 
9257  /* convert 8BSI values to 16BSI */
9258  z = 0;
9259  iYMax = y + nYValid;
9260  iXMax = x + nXValid;
9261  for (iY = y; iY < iYMax; iY++) {
9262  for (iX = x; iX < iXMax; iX++) {
9263  if (rt_band_get_pixel(rtband, iX, iY, &value, NULL) != ES_NONE) {
9264  rterror("rt_raster_to_gdal_mem: Could not get pixel value to convert from 8BSI to 16BSI");
9265  rtdealloc(values);
9266  if (allocBandNums) rtdealloc(bandNums);
9267  if (allocNodataValues) rtdealloc(excludeNodataValues);
9268  GDALClose(ds);
9269  return 0;
9270  }
9271 
9272  values[z++] = rt_util_clamp_to_16BSI(value);
9273  }
9274  }
9275 
9276  /* burn values */
9277  if (GDALRasterIO(
9278  band, GF_Write,
9279  x, y,
9280  nXValid, nYValid,
9281  values, nXValid, nYValid,
9282  gdal_pt,
9283  0, 0
9284  ) != CE_None) {
9285  rterror("rt_raster_to_gdal_mem: Could not write converted 8BSI to 16BSI values to GDAL band");
9286  rtdealloc(values);
9287  if (allocBandNums) rtdealloc(bandNums);
9288  if (allocNodataValues) rtdealloc(excludeNodataValues);
9289  GDALClose(ds);
9290  return 0;
9291  }
9292  }
9293  }
9294 
9295  rtdealloc(values);
9296  }
9297 
9298  /* Add nodata value for band */
9299  if (rt_band_get_hasnodata_flag(rtband) != FALSE && excludeNodataValues[i]) {
9300  rt_band_get_nodata(rtband, &nodata);
9301  if (GDALSetRasterNoDataValue(band, nodata) != CE_None)
9302  rtwarn("rt_raster_to_gdal_mem: Could not set nodata value for band");
9303  RASTER_DEBUGF(3, "nodata value set to %f", GDALGetRasterNoDataValue(band, NULL));
9304  }
9305 
9306 #if POSTGIS_DEBUG_LEVEL > 3
9307  {
9308  GDALRasterBandH _grb = NULL;
9309  double _min;
9310  double _max;
9311  double _mean;
9312  double _stddev;
9313 
9314  _grb = GDALGetRasterBand(ds, i + 1);
9315  GDALComputeRasterStatistics(_grb, FALSE, &_min, &_max, &_mean, &_stddev, NULL, NULL);
9316  RASTER_DEBUGF(4, "GDAL Band %d stats: %f, %f, %f, %f", i + 1, _min, _max, _mean, _stddev);
9317  }
9318 #endif
9319 
9320  }
9321 
9322  /* necessary??? */
9323  GDALFlushCache(ds);
9324 
9325  if (allocBandNums) rtdealloc(bandNums);
9326  if (allocNodataValues) rtdealloc(excludeNodataValues);
9327 
9328  return ds;
9329 }
int rt_raster_get_num_bands(rt_raster raster)
Definition: rt_api.c:5677
void rtdealloc(void *mem)
Definition: rt_api.c:882
void * rt_band_get_data(rt_band band)
Get pointer to raster band data.
Definition: rt_api.c:1710
tuple gt
Definition: window.py:79
Definition: rt_api.h:184
rt_pixtype rt_band_get_pixtype(rt_band band)
Return pixeltype of this band.
Definition: rt_api.c:1900
tuple band
Definition: ovdump.py:57
int rt_pixtype_size(rt_pixtype pixtype)
Return size in bytes of a value in the given pixtype.
Definition: rt_api.c:1097
void rt_raster_get_geotransform_matrix(rt_raster raster, double *gt)
Get 6-element array of raster geotransform matrix.
Definition: rt_api.c:6005
#define RASTER_DEBUG(level, msg)
Definition: rt_api.h:281
int rt_util_gdal_driver_registered(const char *drv)
Definition: rt_api.c:461
void rtwarn(const char *fmt,...)
Definition: rt_api.c:920
rt_pixtype
Definition: rt_api.h:172
#define RASTER_DEBUGF(level, msg,...)
Definition: rt_api.h:285
rt_errorstate rt_band_get_nodata(rt_band band, double *nodata)
Get NODATA value.
Definition: rt_api.c:3058
rt_errorstate rt_band_get_pixel(rt_band band, int x, int y, double *value, int *nodata)
Get pixel value.
Definition: rt_api.c:2549
int count
Definition: genraster.py:57
tuple ds
Definition: pixval.py:66
tuple x
Definition: pixval.py:53
uint16_t rt_raster_get_height(rt_raster raster)
Definition: rt_api.c:5434
void * rtalloc(size_t size)
Raster core memory management functions.
Definition: rt_api.c:867
void rterror(const char *fmt,...)
Raster core error and info handlers.
Definition: rt_api.c:895
rt_band rt_raster_get_band(rt_raster raster, int n)
Return Nth band, or NULL if unavailable.
Definition: rt_api.c:5686
#define FALSE
Definition: dbfopen.c:169
char * rt_util_gdal_convert_sr(const char *srs, int proj4)
Definition: rt_api.c:323
GDALDataType rt_util_pixtype_to_gdal_datatype(rt_pixtype pt)
Convert rt_pixtype to GDALDataType.
Definition: rt_api.c:229
int rt_band_get_hasnodata_flag(rt_band band)
Get hasnodata flag value.
Definition: rt_api.c:2002
tuple y
Definition: pixval.py:54
uint16_t rt_raster_get_width(rt_raster raster)
Definition: rt_api.c:5426
int16_t rt_util_clamp_to_16BSI(double value)
Definition: rt_api.c:114

Here is the call graph for this function:

Here is the caller graph for this function: