PostGIS  2.1.10dev-r@@SVN_REVISION@@
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 9339 of file rt_api.c.

References ovdump::band, ES_NONE, FALSE, window::gt, rt_raster_serialized_t::height, rt_raster_t::height, rt_band_t::height, rt_raster_serialized_t::numBands, 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(), rt_raster_t::srid, rt_raster_serialized_t::width, rt_raster_t::width, rt_band_t::width, pixval::x, and pixval::y.

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

9339  {
9340  rt_raster rast = NULL;
9341  double gt[6] = {0};
9342  CPLErr cplerr;
9343  uint32_t width = 0;
9344  uint32_t height = 0;
9345  uint32_t numBands = 0;
9346  int i = 0;
9347  char *authname = NULL;
9348  char *authcode = NULL;
9349 
9350  GDALRasterBandH gdband = NULL;
9351  GDALDataType gdpixtype = GDT_Unknown;
9352  rt_band band;
9353  int32_t idx;
9354  rt_pixtype pt = PT_END;
9355  uint32_t ptlen = 0;
9356  int hasnodata = 0;
9357  double nodataval;
9358 
9359  int x;
9360  int y;
9361 
9362  int nXBlocks, nYBlocks;
9363  int nXBlockSize, nYBlockSize;
9364  int iXBlock, iYBlock;
9365  int nXValid, nYValid;
9366  int iY;
9367 
9368  uint8_t *values = NULL;
9369  uint32_t valueslen = 0;
9370  uint8_t *ptr = NULL;
9371 
9372  assert(NULL != ds);
9373 
9374  /* raster size */
9375  width = GDALGetRasterXSize(ds);
9376  height = GDALGetRasterYSize(ds);
9377  RASTER_DEBUGF(3, "Raster dimensions (width x height): %d x %d", width, height);
9378 
9379  /* create new raster */
9380  RASTER_DEBUG(3, "Creating new raster");
9381  rast = rt_raster_new(width, height);
9382  if (NULL == rast) {
9383  rterror("rt_raster_from_gdal_dataset: Out of memory allocating new raster");
9384  return NULL;
9385  }
9386  RASTER_DEBUGF(3, "Created raster dimensions (width x height): %d x %d", rast->width, rast->height);
9387 
9388  /* get raster attributes */
9389  cplerr = GDALGetGeoTransform(ds, gt);
9390  if (GDALGetGeoTransform(ds, gt) != CE_None) {
9391  RASTER_DEBUG(4, "Using default geotransform matrix (0, 1, 0, 0, 0, -1)");
9392  gt[0] = 0;
9393  gt[1] = 1;
9394  gt[2] = 0;
9395  gt[3] = 0;
9396  gt[4] = 0;
9397  gt[5] = -1;
9398  }
9399 
9400  /* apply raster attributes */
9402 
9403  RASTER_DEBUGF(3, "Raster geotransform (%f, %f, %f, %f, %f, %f)",
9404  gt[0], gt[1], gt[2], gt[3], gt[4], gt[5]);
9405 
9406  /* srid */
9407  if (rt_util_gdal_sr_auth_info(ds, &authname, &authcode) == ES_NONE) {
9408  if (
9409  authname != NULL &&
9410  strcmp(authname, "EPSG") == 0 &&
9411  authcode != NULL
9412  ) {
9413  rt_raster_set_srid(rast, atoi(authcode));
9414  RASTER_DEBUGF(3, "New raster's SRID = %d", rast->srid);
9415  }
9416 
9417  if (authname != NULL)
9418  rtdealloc(authname);
9419  if (authcode != NULL)
9420  rtdealloc(authcode);
9421  }
9422 
9423  numBands = GDALGetRasterCount(ds);
9424 
9425 #if POSTGIS_DEBUG_LEVEL > 3
9426  for (i = 1; i <= numBands; i++) {
9427  GDALRasterBandH _grb = NULL;
9428  double _min;
9429  double _max;
9430  double _mean;
9431  double _stddev;
9432 
9433  _grb = GDALGetRasterBand(ds, i);
9434  GDALComputeRasterStatistics(_grb, FALSE, &_min, &_max, &_mean, &_stddev, NULL, NULL);
9435  RASTER_DEBUGF(4, "GDAL Band %d stats: %f, %f, %f, %f", i, _min, _max, _mean, _stddev);
9436  }
9437 #endif
9438 
9439  /* copy bands */
9440  for (i = 1; i <= numBands; i++) {
9441  RASTER_DEBUGF(3, "Processing band %d of %d", i, numBands);
9442  gdband = NULL;
9443  gdband = GDALGetRasterBand(ds, i);
9444  if (NULL == gdband) {
9445  rterror("rt_raster_from_gdal_dataset: Could not get GDAL band");
9446  rt_raster_destroy(rast);
9447  return NULL;
9448  }
9449  RASTER_DEBUGF(4, "gdband @ %p", gdband);
9450 
9451  /* pixtype */
9452  gdpixtype = GDALGetRasterDataType(gdband);
9453  RASTER_DEBUGF(4, "gdpixtype, size = %s, %d", GDALGetDataTypeName(gdpixtype), GDALGetDataTypeSize(gdpixtype) / 8);
9454  pt = rt_util_gdal_datatype_to_pixtype(gdpixtype);
9455  if (pt == PT_END) {
9456  rterror("rt_raster_from_gdal_dataset: Unknown pixel type for GDAL band");
9457  rt_raster_destroy(rast);
9458  return NULL;
9459  }
9460  ptlen = rt_pixtype_size(pt);
9461 
9462  /* size: width and height */
9463  width = GDALGetRasterBandXSize(gdband);
9464  height = GDALGetRasterBandYSize(gdband);
9465  RASTER_DEBUGF(3, "GDAL band dimensions (width x height): %d x %d", width, height);
9466 
9467  /* nodata */
9468  nodataval = GDALGetRasterNoDataValue(gdband, &hasnodata);
9469  RASTER_DEBUGF(3, "(hasnodata, nodataval) = (%d, %f)", hasnodata, nodataval);
9470 
9471  /* create band object */
9473  rast, pt,
9474  (hasnodata ? nodataval : 0),
9475  hasnodata, nodataval, rt_raster_get_num_bands(rast)
9476  );
9477  if (idx < 0) {
9478  rterror("rt_raster_from_gdal_dataset: Could not allocate memory for raster band");
9479  rt_raster_destroy(rast);
9480  return NULL;
9481  }
9482  band = rt_raster_get_band(rast, idx);
9483  RASTER_DEBUGF(3, "Created band of dimension (width x height): %d x %d", band->width, band->height);
9484 
9485  /* this makes use of GDAL's "natural" blocks */
9486  GDALGetBlockSize(gdband, &nXBlockSize, &nYBlockSize);
9487  nXBlocks = (width + nXBlockSize - 1) / nXBlockSize;
9488  nYBlocks = (height + nYBlockSize - 1) / nYBlockSize;
9489  RASTER_DEBUGF(4, "(nXBlockSize, nYBlockSize) = (%d, %d)", nXBlockSize, nYBlockSize);
9490  RASTER_DEBUGF(4, "(nXBlocks, nYBlocks) = (%d, %d)", nXBlocks, nYBlocks);
9491 
9492  /* allocate memory for values */
9493  valueslen = ptlen * nXBlockSize * nYBlockSize;
9494  values = rtalloc(valueslen);
9495  if (values == NULL) {
9496  rterror("rt_raster_from_gdal_dataset: Could not allocate memory for GDAL band pixel values");
9497  rt_raster_destroy(rast);
9498  return NULL;
9499  }
9500  RASTER_DEBUGF(3, "values @ %p of length = %d", values, valueslen);
9501 
9502  for (iYBlock = 0; iYBlock < nYBlocks; iYBlock++) {
9503  for (iXBlock = 0; iXBlock < nXBlocks; iXBlock++) {
9504  x = iXBlock * nXBlockSize;
9505  y = iYBlock * nYBlockSize;
9506  RASTER_DEBUGF(4, "(iXBlock, iYBlock) = (%d, %d)", iXBlock, iYBlock);
9507  RASTER_DEBUGF(4, "(x, y) = (%d, %d)", x, y);
9508 
9509  memset(values, 0, valueslen);
9510 
9511  /* valid block width */
9512  if ((iXBlock + 1) * nXBlockSize > width)
9513  nXValid = width - (iXBlock * nXBlockSize);
9514  else
9515  nXValid = nXBlockSize;
9516 
9517  /* valid block height */
9518  if ((iYBlock + 1) * nYBlockSize > height)
9519  nYValid = height - (iYBlock * nYBlockSize);
9520  else
9521  nYValid = nYBlockSize;
9522 
9523  RASTER_DEBUGF(4, "(nXValid, nYValid) = (%d, %d)", nXValid, nYValid);
9524 
9525  cplerr = GDALRasterIO(
9526  gdband, GF_Read,
9527  x, y,
9528  nXValid, nYValid,
9529  values, nXValid, nYValid,
9530  gdpixtype,
9531  0, 0
9532  );
9533  if (cplerr != CE_None) {
9534  rterror("rt_raster_from_gdal_dataset: Could not get data from GDAL raster");
9535  rtdealloc(values);
9536  rt_raster_destroy(rast);
9537  return NULL;
9538  }
9539 
9540  /* if block width is same as raster width, shortcut */
9541  if (nXBlocks == 1 && nYBlockSize > 1 && nXValid == width) {
9542  x = 0;
9543  y = nYBlockSize * iYBlock;
9544 
9545  RASTER_DEBUGF(4, "Setting set of pixel lines at (%d, %d) for %d pixels", x, y, nXValid * nYValid);
9546  rt_band_set_pixel_line(band, x, y, values, nXValid * nYValid);
9547  }
9548  else {
9549  ptr = values;
9550  x = nXBlockSize * iXBlock;
9551  for (iY = 0; iY < nYValid; iY++) {
9552  y = iY + (nYBlockSize * iYBlock);
9553 
9554  RASTER_DEBUGF(4, "Setting pixel line at (%d, %d) for %d pixels", x, y, nXValid);
9555  rt_band_set_pixel_line(band, x, y, ptr, nXValid);
9556  ptr += (nXValid * ptlen);
9557  }
9558  }
9559  }
9560  }
9561 
9562  /* free memory */
9563  rtdealloc(values);
9564  }
9565 
9566  return rast;
9567 }
int rt_raster_get_num_bands(rt_raster raster)
Definition: rt_api.c:5677
int32_t srid
Definition: rt_api.h:2225
void rtdealloc(void *mem)
Definition: rt_api.c:882
void rt_raster_destroy(rt_raster raster)
Release memory associated to a raster.
Definition: rt_api.c:5387
tuple gt
Definition: window.py:79
Definition: rt_api.h:184
tuple band
Definition: ovdump.py:57
tuple rast
Definition: rtpixdump.py:62
uint16_t height
Definition: rt_api.h:2242
rt_errorstate rt_band_set_pixel_line(rt_band band, int x, int y, void *vals, uint32_t len)
Set values of multiple pixels.
Definition: rt_api.c:2181
void rt_raster_set_srid(rt_raster raster, int32_t srid)
Set raster's SRID.
Definition: rt_api.c:5668
uint16_t height
Definition: rt_api.h:2227
int rt_pixtype_size(rt_pixtype pixtype)
Return size in bytes of a value in the given pixtype.
Definition: rt_api.c:1097
#define RASTER_DEBUG(level, msg)
Definition: rt_api.h:281
rt_pixtype rt_util_gdal_datatype_to_pixtype(GDALDataType gdt)
Convert GDALDataType to rt_pixtype.
Definition: rt_api.c:264
rt_pixtype
Definition: rt_api.h:172
#define RASTER_DEBUGF(level, msg,...)
Definition: rt_api.h:285
uint16_t width
Definition: rt_api.h:2226
void rt_raster_set_geotransform_matrix(rt_raster raster, double *gt)
Set raster's geotransform using 6-element array.
Definition: rt_api.c:6026
tuple ds
Definition: pixval.py:66
rt_errorstate rt_util_gdal_sr_auth_info(GDALDatasetH hds, char **authname, char **authcode)
Get auth name and code.
Definition: rt_api.c:381
uint16_t width
Definition: rt_api.h:2241
tuple x
Definition: pixval.py:53
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
rt_raster rt_raster_new(uint32_t width, uint32_t height)
Construct a raster with given dimensions.
Definition: rt_api.c:5353
tuple y
Definition: pixval.py:54
int rt_raster_generate_new_band(rt_raster raster, rt_pixtype pixtype, double initialvalue, uint32_t hasnodata, double nodatavalue, int index)
Generate a new inline band and add it to a raster.
Definition: rt_api.c:5784

Here is the call graph for this function:

Here is the caller graph for this function: