PostGIS  2.5.2dev-r@@SVN_REVISION@@

◆ 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 2182 of file rt_raster.c.

References ovdump::band, ES_NONE, FALSE, window::gt, rt_raster_t::height, rt_band_t::height, 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_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().

2182  {
2183  rt_raster rast = NULL;
2184  double gt[6] = {0};
2185  CPLErr cplerr;
2186  uint32_t width = 0;
2187  uint32_t height = 0;
2188  uint32_t numBands = 0;
2189  uint32_t i = 0;
2190  char *authname = NULL;
2191  char *authcode = NULL;
2192 
2193  GDALRasterBandH gdband = NULL;
2194  GDALDataType gdpixtype = GDT_Unknown;
2195  rt_band band;
2196  int32_t idx;
2197  rt_pixtype pt = PT_END;
2198  uint32_t ptlen = 0;
2199  int hasnodata = 0;
2200  double nodataval;
2201 
2202  int x;
2203  int y;
2204 
2205  uint32_t nXBlocks, nYBlocks;
2206  int nXBlockSize, nYBlockSize;
2207  uint32_t iXBlock, iYBlock;
2208  uint32_t nXValid, nYValid;
2209  uint32_t iY;
2210 
2211  uint8_t *values = NULL;
2212  uint32_t valueslen = 0;
2213  uint8_t *ptr = NULL;
2214 
2215  assert(NULL != ds);
2216 
2217  /* raster size */
2218  width = GDALGetRasterXSize(ds);
2219  height = GDALGetRasterYSize(ds);
2220  RASTER_DEBUGF(3, "Raster dimensions (width x height): %d x %d", width, height);
2221 
2222  /* create new raster */
2223  RASTER_DEBUG(3, "Creating new raster");
2224  rast = rt_raster_new(width, height);
2225  if (NULL == rast) {
2226  rterror("rt_raster_from_gdal_dataset: Out of memory allocating new raster");
2227  return NULL;
2228  }
2229  RASTER_DEBUGF(3, "Created raster dimensions (width x height): %d x %d", rast->width, rast->height);
2230 
2231  /* get raster attributes */
2232  cplerr = GDALGetGeoTransform(ds, gt);
2233  if (GDALGetGeoTransform(ds, gt) != CE_None) {
2234  RASTER_DEBUG(4, "Using default geotransform matrix (0, 1, 0, 0, 0, -1)");
2235  gt[0] = 0;
2236  gt[1] = 1;
2237  gt[2] = 0;
2238  gt[3] = 0;
2239  gt[4] = 0;
2240  gt[5] = -1;
2241  }
2242 
2243  /* apply raster attributes */
2245 
2246  RASTER_DEBUGF(3, "Raster geotransform (%f, %f, %f, %f, %f, %f)",
2247  gt[0], gt[1], gt[2], gt[3], gt[4], gt[5]);
2248 
2249  /* srid */
2250  if (rt_util_gdal_sr_auth_info(ds, &authname, &authcode) == ES_NONE) {
2251  if (
2252  authname != NULL &&
2253  strcmp(authname, "EPSG") == 0 &&
2254  authcode != NULL
2255  ) {
2256  rt_raster_set_srid(rast, atoi(authcode));
2257  RASTER_DEBUGF(3, "New raster's SRID = %d", rast->srid);
2258  }
2259 
2260  if (authname != NULL)
2261  rtdealloc(authname);
2262  if (authcode != NULL)
2263  rtdealloc(authcode);
2264  }
2265 
2266  numBands = GDALGetRasterCount(ds);
2267 
2268 #if POSTGIS_DEBUG_LEVEL > 3
2269  for (i = 1; i <= numBands; i++) {
2270  GDALRasterBandH _grb = NULL;
2271  double _min;
2272  double _max;
2273  double _mean;
2274  double _stddev;
2275 
2276  _grb = GDALGetRasterBand(ds, i);
2277  GDALComputeRasterStatistics(_grb, FALSE, &_min, &_max, &_mean, &_stddev, NULL, NULL);
2278  RASTER_DEBUGF(4, "GDAL Band %d stats: %f, %f, %f, %f", i, _min, _max, _mean, _stddev);
2279  }
2280 #endif
2281 
2282  /* copy bands */
2283  for (i = 1; i <= numBands; i++) {
2284  RASTER_DEBUGF(3, "Processing band %d of %d", i, numBands);
2285  gdband = NULL;
2286  gdband = GDALGetRasterBand(ds, i);
2287  if (NULL == gdband) {
2288  rterror("rt_raster_from_gdal_dataset: Could not get GDAL band");
2289  rt_raster_destroy(rast);
2290  return NULL;
2291  }
2292  RASTER_DEBUGF(4, "gdband @ %p", gdband);
2293 
2294  /* pixtype */
2295  gdpixtype = GDALGetRasterDataType(gdband);
2296  RASTER_DEBUGF(4, "gdpixtype, size = %s, %d", GDALGetDataTypeName(gdpixtype), GDALGetDataTypeSize(gdpixtype) / 8);
2297  pt = rt_util_gdal_datatype_to_pixtype(gdpixtype);
2298  if (pt == PT_END) {
2299  rterror("rt_raster_from_gdal_dataset: Unknown pixel type for GDAL band");
2300  rt_raster_destroy(rast);
2301  return NULL;
2302  }
2303  ptlen = rt_pixtype_size(pt);
2304 
2305  /* size: width and height */
2306  width = GDALGetRasterBandXSize(gdband);
2307  height = GDALGetRasterBandYSize(gdband);
2308  RASTER_DEBUGF(3, "GDAL band dimensions (width x height): %d x %d", width, height);
2309 
2310  /* nodata */
2311  nodataval = GDALGetRasterNoDataValue(gdband, &hasnodata);
2312  RASTER_DEBUGF(3, "(hasnodata, nodataval) = (%d, %f)", hasnodata, nodataval);
2313 
2314  /* create band object */
2316  rast, pt,
2317  (hasnodata ? nodataval : 0),
2318  hasnodata, nodataval, rt_raster_get_num_bands(rast)
2319  );
2320  if (idx < 0) {
2321  rterror("rt_raster_from_gdal_dataset: Could not allocate memory for raster band");
2322  rt_raster_destroy(rast);
2323  return NULL;
2324  }
2325  band = rt_raster_get_band(rast, idx);
2326  RASTER_DEBUGF(3, "Created band of dimension (width x height): %d x %d", band->width, band->height);
2327 
2328  /* this makes use of GDAL's "natural" blocks */
2329  GDALGetBlockSize(gdband, &nXBlockSize, &nYBlockSize);
2330  nXBlocks = (width + nXBlockSize - 1) / nXBlockSize;
2331  nYBlocks = (height + nYBlockSize - 1) / nYBlockSize;
2332  RASTER_DEBUGF(4, "(nXBlockSize, nYBlockSize) = (%d, %d)", nXBlockSize, nYBlockSize);
2333  RASTER_DEBUGF(4, "(nXBlocks, nYBlocks) = (%d, %d)", nXBlocks, nYBlocks);
2334 
2335  /* allocate memory for values */
2336  valueslen = ptlen * nXBlockSize * nYBlockSize;
2337  values = rtalloc(valueslen);
2338  if (values == NULL) {
2339  rterror("rt_raster_from_gdal_dataset: Could not allocate memory for GDAL band pixel values");
2340  rt_raster_destroy(rast);
2341  return NULL;
2342  }
2343  RASTER_DEBUGF(3, "values @ %p of length = %d", values, valueslen);
2344 
2345  for (iYBlock = 0; iYBlock < nYBlocks; iYBlock++) {
2346  for (iXBlock = 0; iXBlock < nXBlocks; iXBlock++) {
2347  x = iXBlock * nXBlockSize;
2348  y = iYBlock * nYBlockSize;
2349  RASTER_DEBUGF(4, "(iXBlock, iYBlock) = (%d, %d)", iXBlock, iYBlock);
2350  RASTER_DEBUGF(4, "(x, y) = (%d, %d)", x, y);
2351 
2352  memset(values, 0, valueslen);
2353 
2354  /* valid block width */
2355  if ((iXBlock + 1) * nXBlockSize > width)
2356  nXValid = width - (iXBlock * nXBlockSize);
2357  else
2358  nXValid = nXBlockSize;
2359 
2360  /* valid block height */
2361  if ((iYBlock + 1) * nYBlockSize > height)
2362  nYValid = height - (iYBlock * nYBlockSize);
2363  else
2364  nYValid = nYBlockSize;
2365 
2366  RASTER_DEBUGF(4, "(nXValid, nYValid) = (%d, %d)", nXValid, nYValid);
2367 
2368  cplerr = GDALRasterIO(
2369  gdband, GF_Read,
2370  x, y,
2371  nXValid, nYValid,
2372  values, nXValid, nYValid,
2373  gdpixtype,
2374  0, 0
2375  );
2376  if (cplerr != CE_None) {
2377  rterror("rt_raster_from_gdal_dataset: Could not get data from GDAL raster");
2378  rtdealloc(values);
2379  rt_raster_destroy(rast);
2380  return NULL;
2381  }
2382 
2383  /* if block width is same as raster width, shortcut */
2384  if (nXBlocks == 1 && nYBlockSize > 1 && nXValid == width) {
2385  x = 0;
2386  y = nYBlockSize * iYBlock;
2387 
2388  RASTER_DEBUGF(4, "Setting set of pixel lines at (%d, %d) for %d pixels", x, y, nXValid * nYValid);
2389  rt_band_set_pixel_line(band, x, y, values, nXValid * nYValid);
2390  }
2391  else {
2392  ptr = values;
2393  x = nXBlockSize * iXBlock;
2394  for (iY = 0; iY < nYValid; iY++) {
2395  y = iY + (nYBlockSize * iYBlock);
2396 
2397  RASTER_DEBUGF(4, "Setting pixel line at (%d, %d) for %d pixels", x, y, nXValid);
2398  rt_band_set_pixel_line(band, x, y, ptr, nXValid);
2399  ptr += (nXValid * ptlen);
2400  }
2401  }
2402  }
2403  }
2404 
2405  /* free memory */
2406  rtdealloc(values);
2407  }
2408 
2409  return rast;
2410 }
int32_t srid
Definition: librtcore.h:2300
band
Definition: ovdump.py:57
void rterror(const char *fmt,...)
Wrappers used for reporting errors and info.
Definition: rt_context.c:199
void * rtalloc(size_t size)
Wrappers used for managing memory.
Definition: rt_context.c:171
uint16_t height
Definition: librtcore.h:2317
void rt_raster_set_srid(rt_raster raster, int32_t srid)
Set raster&#39;s SRID.
Definition: rt_raster.c:363
gt
Definition: window.py:77
rt_pixtype
Definition: librtcore.h:185
uint16_t rt_raster_get_num_bands(rt_raster raster)
Definition: rt_raster.c:372
void rt_raster_set_geotransform_matrix(rt_raster raster, double *gt)
Set raster&#39;s geotransform using 6-element array.
Definition: rt_raster.c:727
uint16_t height
Definition: librtcore.h:2302
ds
Definition: pixval.py:66
unsigned int uint32_t
Definition: uthash.h:78
uint16_t width
Definition: librtcore.h:2301
rt_pixtype rt_util_gdal_datatype_to_pixtype(GDALDataType gdt)
Convert GDALDataType to rt_pixtype.
Definition: rt_util.c:153
#define RASTER_DEBUGF(level, msg,...)
Definition: librtcore.h:299
uint16_t width
Definition: librtcore.h:2316
int rt_pixtype_size(rt_pixtype pixtype)
Return size in bytes of a value in the given pixtype.
Definition: rt_pixel.c:39
rt_raster rt_raster_new(uint32_t width, uint32_t height)
Construct a raster with given dimensions.
Definition: rt_raster.c:48
void rtdealloc(void *mem)
Definition: rt_context.c:186
#define FALSE
Definition: dbfopen.c:168
#define RASTER_DEBUG(level, msg)
Definition: librtcore.h:295
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_band.c:853
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_raster.c:485
rt_errorstate rt_util_gdal_sr_auth_info(GDALDatasetH hds, char **authname, char **authcode)
Get auth name and code.
Definition: rt_util.c:270
unsigned char uint8_t
Definition: uthash.h:79
void rt_raster_destroy(rt_raster raster)
Release memory associated to a raster.
Definition: rt_raster.c:82
rt_band rt_raster_get_band(rt_raster raster, int n)
Return Nth band, or NULL if unavailable.
Definition: rt_raster.c:381
Here is the call graph for this function:
Here is the caller graph for this function: