PostGIS  3.4.0dev-r@@SVN_REVISION@@

◆ rt_raster_to_gdal_mem()

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 to check for pixels with value
count: number of elements in bandNums and exclude_nodata_values
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
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 1935 of file rt_raster.c.

1942  {
1943  GDALDriverH drv = NULL;
1944  GDALDatasetH ds = NULL;
1945  double gt[6] = {0.0};
1946  CPLErr cplerr;
1947  GDALDataType gdal_pt = GDT_Unknown;
1948  GDALRasterBandH band;
1949  void *pVoid;
1950  char *pszDataPointer;
1951  char szGDALOption[50];
1952  char *apszOptions[4];
1953  double nodata = 0.0;
1954  int allocBandNums = 0;
1955  int allocNodataValues = 0;
1956 
1957  int i;
1958  uint32_t numBands;
1959  uint32_t width = 0;
1960  uint32_t height = 0;
1961  rt_band rtband = NULL;
1962  rt_pixtype pt = PT_END;
1963 
1964  assert(NULL != raster);
1965  assert(NULL != rtn_drv);
1966  assert(NULL != destroy_rtn_drv);
1967 
1968  *destroy_rtn_drv = 0;
1969 
1970  /* store raster in GDAL MEM raster */
1971  if (!rt_util_gdal_driver_registered("MEM")) {
1972  RASTER_DEBUG(4, "Registering MEM driver");
1973  GDALRegister_MEM();
1974  *destroy_rtn_drv = 1;
1975  }
1976  drv = GDALGetDriverByName("MEM");
1977  if (NULL == drv) {
1978  rterror("rt_raster_to_gdal_mem: Could not load the MEM GDAL driver");
1979  return 0;
1980  }
1981  *rtn_drv = drv;
1982 
1983  /* unload driver from GDAL driver manager */
1984  if (*destroy_rtn_drv) {
1985  RASTER_DEBUG(4, "Deregistering MEM driver");
1986  GDALDeregisterDriver(drv);
1987  }
1988 
1989  width = rt_raster_get_width(raster);
1990  height = rt_raster_get_height(raster);
1991  ds = GDALCreate(
1992  drv, "",
1993  width, height,
1994  0, GDT_Byte, NULL
1995  );
1996  if (NULL == ds) {
1997  rterror("rt_raster_to_gdal_mem: Could not create a GDALDataset to convert into");
1998  return 0;
1999  }
2000 
2001  /* add geotransform */
2002  rt_raster_get_geotransform_matrix(raster, gt);
2003  cplerr = GDALSetGeoTransform(ds, gt);
2004  if (cplerr != CE_None) {
2005  rterror("rt_raster_to_gdal_mem: Could not set geotransformation");
2006  GDALClose(ds);
2007  return 0;
2008  }
2009 
2010  /* set spatial reference */
2011  if (NULL != srs && strlen(srs)) {
2012  char *_srs = rt_util_gdal_convert_sr(srs, 0);
2013  if (_srs == NULL) {
2014  rterror("rt_raster_to_gdal_mem: Could not convert srs to GDAL accepted format");
2015  GDALClose(ds);
2016  return 0;
2017  }
2018 
2019  cplerr = GDALSetProjection(ds, _srs);
2020  CPLFree(_srs);
2021  if (cplerr != CE_None) {
2022  rterror("rt_raster_to_gdal_mem: Could not set projection");
2023  GDALClose(ds);
2024  return 0;
2025  }
2026  RASTER_DEBUGF(3, "Projection set to: %s", GDALGetProjectionRef(ds));
2027  }
2028 
2029  /* process bandNums */
2030  numBands = rt_raster_get_num_bands(raster);
2031  if (NULL != bandNums && count > 0) {
2032  for (i = 0; i < count; i++) {
2033  if (bandNums[i] >= numBands) {
2034  rterror("rt_raster_to_gdal_mem: The band index %d is invalid", bandNums[i]);
2035  GDALClose(ds);
2036  return 0;
2037  }
2038  }
2039  }
2040  else {
2041  count = numBands;
2042  bandNums = (uint32_t *) rtalloc(sizeof(uint32_t) * count);
2043  if (NULL == bandNums) {
2044  rterror("rt_raster_to_gdal_mem: Could not allocate memory for band indices");
2045  GDALClose(ds);
2046  return 0;
2047  }
2048  allocBandNums = 1;
2049  for (i = 0; i < count; i++) bandNums[i] = i;
2050  }
2051 
2052  /* process exclude_nodata_values */
2053  if (NULL == excludeNodataValues) {
2054  excludeNodataValues = (int *) rtalloc(sizeof(int) * count);
2055  if (NULL == excludeNodataValues) {
2056  rterror("rt_raster_to_gdal_mem: Could not allocate memory for NODATA flags");
2057  GDALClose(ds);
2058  return 0;
2059  }
2060  allocNodataValues = 1;
2061  for (i = 0; i < count; i++) excludeNodataValues[i] = 1;
2062  }
2063 
2064  /* add band(s) */
2065  for (i = 0; i < count; i++) {
2066  rtband = rt_raster_get_band(raster, bandNums[i]);
2067  if (NULL == rtband) {
2068  rterror("rt_raster_to_gdal_mem: Could not get requested band index %d", bandNums[i]);
2069  if (allocBandNums) rtdealloc(bandNums);
2070  if (allocNodataValues) rtdealloc(excludeNodataValues);
2071  GDALClose(ds);
2072  return 0;
2073  }
2074 
2075  pt = rt_band_get_pixtype(rtband);
2076  gdal_pt = rt_util_pixtype_to_gdal_datatype(pt);
2077  if (gdal_pt == GDT_Unknown)
2078  rtwarn("rt_raster_to_gdal_mem: Unknown pixel type for band");
2079 
2080  /*
2081  For all pixel types other than PT_8BSI, set pointer to start of data
2082  */
2083  if (pt != PT_8BSI) {
2084  pVoid = rt_band_get_data(rtband);
2085  RASTER_DEBUGF(4, "Band data is at pos %p", pVoid);
2086 
2087  pszDataPointer = (char *) rtalloc(20 * sizeof (char));
2088  sprintf(pszDataPointer, "%p", pVoid);
2089  RASTER_DEBUGF(4, "rt_raster_to_gdal_mem: szDatapointer is %p",
2090  pszDataPointer);
2091 
2092  if (strncasecmp(pszDataPointer, "0x", 2) == 0)
2093  sprintf(szGDALOption, "DATAPOINTER=%s", pszDataPointer);
2094  else
2095  sprintf(szGDALOption, "DATAPOINTER=0x%s", pszDataPointer);
2096 
2097  RASTER_DEBUG(3, "Storing info for GDAL MEM raster band");
2098 
2099  apszOptions[0] = szGDALOption;
2100  apszOptions[1] = NULL; /* pixel offset, not needed */
2101  apszOptions[2] = NULL; /* line offset, not needed */
2102  apszOptions[3] = NULL;
2103 
2104  /* free */
2105  rtdealloc(pszDataPointer);
2106 
2107  /* add band */
2108  if (GDALAddBand(ds, gdal_pt, apszOptions) == CE_Failure) {
2109  rterror("rt_raster_to_gdal_mem: Could not add GDAL raster band");
2110  if (allocBandNums) rtdealloc(bandNums);
2111  GDALClose(ds);
2112  return 0;
2113  }
2114  }
2115  /*
2116  PT_8BSI is special as GDAL has no equivalent pixel type.
2117  Must convert 8BSI to 16BSI so create basic band
2118  */
2119  else {
2120  /* add band */
2121  if (GDALAddBand(ds, gdal_pt, NULL) == CE_Failure) {
2122  rterror("rt_raster_to_gdal_mem: Could not add GDAL raster band");
2123  if (allocBandNums) rtdealloc(bandNums);
2124  if (allocNodataValues) rtdealloc(excludeNodataValues);
2125  GDALClose(ds);
2126  return 0;
2127  }
2128  }
2129 
2130  /* check band count */
2131  if (GDALGetRasterCount(ds) != i + 1) {
2132  rterror("rt_raster_to_gdal_mem: Error creating GDAL MEM raster band");
2133  if (allocBandNums) rtdealloc(bandNums);
2134  if (allocNodataValues) rtdealloc(excludeNodataValues);
2135  GDALClose(ds);
2136  return 0;
2137  }
2138 
2139  /* get new band */
2140  band = NULL;
2141  band = GDALGetRasterBand(ds, i + 1);
2142  if (NULL == band) {
2143  rterror("rt_raster_to_gdal_mem: Could not get GDAL band for additional processing");
2144  if (allocBandNums) rtdealloc(bandNums);
2145  if (allocNodataValues) rtdealloc(excludeNodataValues);
2146  GDALClose(ds);
2147  return 0;
2148  }
2149 
2150  /* PT_8BSI requires manual setting of pixels */
2151  if (pt == PT_8BSI) {
2152  uint32_t nXBlocks, nYBlocks;
2153  int nXBlockSize, nYBlockSize;
2154  uint32_t iXBlock, iYBlock;
2155  int nXValid, nYValid;
2156  int iX, iY;
2157  int iXMax, iYMax;
2158 
2159  int x, y, z;
2160  uint32_t valueslen = 0;
2161  int16_t *values = NULL;
2162  double value = 0.;
2163 
2164  /* this makes use of GDAL's "natural" blocks */
2165  GDALGetBlockSize(band, &nXBlockSize, &nYBlockSize);
2166  nXBlocks = (width + nXBlockSize - 1) / nXBlockSize;
2167  nYBlocks = (height + nYBlockSize - 1) / nYBlockSize;
2168  RASTER_DEBUGF(4, "(nXBlockSize, nYBlockSize) = (%d, %d)", nXBlockSize, nYBlockSize);
2169  RASTER_DEBUGF(4, "(nXBlocks, nYBlocks) = (%d, %d)", nXBlocks, nYBlocks);
2170 
2171  /* length is for the desired pixel type */
2172  valueslen = rt_pixtype_size(PT_16BSI) * nXBlockSize * nYBlockSize;
2173  values = rtalloc(valueslen);
2174  if (NULL == values) {
2175  rterror("rt_raster_to_gdal_mem: Could not allocate memory for GDAL band pixel values");
2176  if (allocBandNums) rtdealloc(bandNums);
2177  if (allocNodataValues) rtdealloc(excludeNodataValues);
2178  GDALClose(ds);
2179  return 0;
2180  }
2181 
2182  for (iYBlock = 0; iYBlock < nYBlocks; iYBlock++) {
2183  for (iXBlock = 0; iXBlock < nXBlocks; iXBlock++) {
2184  memset(values, 0, valueslen);
2185 
2186  x = iXBlock * nXBlockSize;
2187  y = iYBlock * nYBlockSize;
2188  RASTER_DEBUGF(4, "(iXBlock, iYBlock) = (%d, %d)", iXBlock, iYBlock);
2189  RASTER_DEBUGF(4, "(x, y) = (%d, %d)", x, y);
2190 
2191  /* valid block width */
2192  if ((iXBlock + 1) * nXBlockSize > width)
2193  nXValid = width - (iXBlock * nXBlockSize);
2194  else
2195  nXValid = nXBlockSize;
2196 
2197  /* valid block height */
2198  if ((iYBlock + 1) * nYBlockSize > height)
2199  nYValid = height - (iYBlock * nYBlockSize);
2200  else
2201  nYValid = nYBlockSize;
2202 
2203  RASTER_DEBUGF(4, "(nXValid, nYValid) = (%d, %d)", nXValid, nYValid);
2204 
2205  /* convert 8BSI values to 16BSI */
2206  z = 0;
2207  iYMax = y + nYValid;
2208  iXMax = x + nXValid;
2209  for (iY = y; iY < iYMax; iY++) {
2210  for (iX = x; iX < iXMax; iX++) {
2211  if (rt_band_get_pixel(rtband, iX, iY, &value, NULL) != ES_NONE) {
2212  rterror("rt_raster_to_gdal_mem: Could not get pixel value to convert from 8BSI to 16BSI");
2213  rtdealloc(values);
2214  if (allocBandNums) rtdealloc(bandNums);
2215  if (allocNodataValues) rtdealloc(excludeNodataValues);
2216  GDALClose(ds);
2217  return 0;
2218  }
2219 
2220  values[z++] = rt_util_clamp_to_16BSI(value);
2221  }
2222  }
2223 
2224  /* burn values */
2225  if (GDALRasterIO(
2226  band, GF_Write,
2227  x, y,
2228  nXValid, nYValid,
2229  values, nXValid, nYValid,
2230  gdal_pt,
2231  0, 0
2232  ) != CE_None) {
2233  rterror("rt_raster_to_gdal_mem: Could not write converted 8BSI to 16BSI values to GDAL band");
2234  rtdealloc(values);
2235  if (allocBandNums) rtdealloc(bandNums);
2236  if (allocNodataValues) rtdealloc(excludeNodataValues);
2237  GDALClose(ds);
2238  return 0;
2239  }
2240  }
2241  }
2242 
2243  rtdealloc(values);
2244  }
2245 
2246  /* Add nodata value for band */
2247  if (rt_band_get_hasnodata_flag(rtband) != FALSE && excludeNodataValues[i]) {
2248  rt_band_get_nodata(rtband, &nodata);
2249  if (GDALSetRasterNoDataValue(band, nodata) != CE_None)
2250  rtwarn("rt_raster_to_gdal_mem: Could not set nodata value for band");
2251  RASTER_DEBUGF(3, "nodata value set to %f", GDALGetRasterNoDataValue(band, NULL));
2252  }
2253 
2254 #if POSTGIS_DEBUG_LEVEL > 3
2255  {
2256  GDALRasterBandH _grb = NULL;
2257  double _min;
2258  double _max;
2259  double _mean;
2260  double _stddev;
2261 
2262  _grb = GDALGetRasterBand(ds, i + 1);
2263  GDALComputeRasterStatistics(_grb, FALSE, &_min, &_max, &_mean, &_stddev, NULL, NULL);
2264  RASTER_DEBUGF(4, "GDAL Band %d stats: %f, %f, %f, %f", i + 1, _min, _max, _mean, _stddev);
2265  }
2266 #endif
2267 
2268  }
2269 
2270  /* necessary??? */
2271  GDALFlushCache(ds);
2272 
2273  if (allocBandNums) rtdealloc(bandNums);
2274  if (allocNodataValues) rtdealloc(excludeNodataValues);
2275 
2276  return ds;
2277 }
FALSE
#define FALSE
Definition: dbfopen.c:72
rterror
void rterror(const char *fmt,...)
Wrappers used for reporting errors and info.
Definition: rt_context.c:219
rtalloc
void * rtalloc(size_t size)
Wrappers used for managing memory.
Definition: rt_context.c:191
RASTER_DEBUG
#define RASTER_DEBUG(level, msg)
Definition: librtcore.h:302
RASTER_DEBUGF
#define RASTER_DEBUGF(level, msg,...)
Definition: librtcore.h:306
rt_util_gdal_convert_sr
char * rt_util_gdal_convert_sr(const char *srs, int proj4)
Definition: rt_util.c:218
rt_band_get_hasnodata_flag
int rt_band_get_hasnodata_flag(rt_band band)
Get hasnodata flag value.
Definition: rt_band.c:674
rt_band_get_pixel
rt_errorstate rt_band_get_pixel(rt_band band, int x, int y, double *value, int *nodata)
Get pixel value.
Definition: rt_band.c:1376
rt_pixtype
rt_pixtype
Definition: librtcore.h:187
PT_END
@ PT_END
Definition: librtcore.h:199
PT_8BSI
@ PT_8BSI
Definition: librtcore.h:191
PT_16BSI
@ PT_16BSI
Definition: librtcore.h:193
rt_util_pixtype_to_gdal_datatype
GDALDataType rt_util_pixtype_to_gdal_datatype(rt_pixtype pt)
Convert rt_pixtype to GDALDataType.
Definition: rt_util.c:124
rtwarn
void rtwarn(const char *fmt,...)
Definition: rt_context.c:244
ES_NONE
@ ES_NONE
Definition: librtcore.h:182
rt_util_clamp_to_16BSI
int16_t rt_util_clamp_to_16BSI(double value)
Definition: rt_util.c:60
rt_util_gdal_driver_registered
int rt_util_gdal_driver_registered(const char *drv)
Definition: rt_util.c:361
rt_band_get_nodata
rt_errorstate rt_band_get_nodata(rt_band band, double *nodata)
Get NODATA value.
Definition: rt_band.c:1887
rt_band_get_pixtype
rt_pixtype rt_band_get_pixtype(rt_band band)
Return pixeltype of this band.
Definition: rt_band.c:631
rtdealloc
void rtdealloc(void *mem)
Definition: rt_context.c:206
rt_band_get_data
void * rt_band_get_data(rt_band band)
Get pointer to raster band data.
Definition: rt_band.c:400
rt_pixtype_size
int rt_pixtype_size(rt_pixtype pixtype)
Return size in bytes of a value in the given pixtype.
Definition: rt_pixel.c:39
genraster.value
int value
Definition: genraster.py:62
genraster.count
int count
Definition: genraster.py:57
ovdump.band
band
Definition: ovdump.py:58
pixval.ds
ds
Definition: pixval.py:67
rtrowdump.raster
raster
Be careful!! Zeros function's input parameter can be a (height x width) array, not (width x height): ...
Definition: rtrowdump.py:121
window.gt
gt
Definition: window.py:78
rt_raster_get_num_bands
uint16_t rt_raster_get_num_bands(rt_raster raster)
Definition: rt_raster.c:376
rt_raster_get_height
uint16_t rt_raster_get_height(rt_raster raster)
Definition: rt_raster.c:133
rt_raster_get_band
rt_band rt_raster_get_band(rt_raster raster, int n)
Return Nth band, or NULL if unavailable.
Definition: rt_raster.c:385
rt_raster_get_width
uint16_t rt_raster_get_width(rt_raster raster)
Definition: rt_raster.c:125
rt_raster_get_geotransform_matrix
void rt_raster_get_geotransform_matrix(rt_raster raster, double *gt)
Get 6-element array of raster geotransform matrix.
Definition: rt_raster.c:710
rt_band_t
Definition: librtcore.h:2466

References ovdump::band, genraster::count, pixval::ds, ES_NONE, FALSE, window::gt, PT_16BSI, PT_8BSI, PT_END, rtrowdump::raster, 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, pixval::x, and pixval::y.

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

Here is the call graph for this function:
Here is the caller graph for this function: