PostGIS  3.1.6dev-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 2164 of file rt_raster.c.

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

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: