PostGIS  2.5.0dev-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 1830 of file rt_raster.c.

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

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

1837  {
1838  GDALDriverH drv = NULL;
1839  GDALDatasetH ds = NULL;
1840  double gt[6] = {0.0};
1841  CPLErr cplerr;
1842  GDALDataType gdal_pt = GDT_Unknown;
1843  GDALRasterBandH band;
1844  void *pVoid;
1845  char *pszDataPointer;
1846  char szGDALOption[50];
1847  char *apszOptions[4];
1848  double nodata = 0.0;
1849  int allocBandNums = 0;
1850  int allocNodataValues = 0;
1851 
1852  int i;
1853  uint32_t numBands;
1854  uint32_t width = 0;
1855  uint32_t height = 0;
1856  rt_band rtband = NULL;
1857  rt_pixtype pt = PT_END;
1858 
1859  assert(NULL != raster);
1860  assert(NULL != rtn_drv);
1861  assert(NULL != destroy_rtn_drv);
1862 
1863  *destroy_rtn_drv = 0;
1864 
1865  /* store raster in GDAL MEM raster */
1866  if (!rt_util_gdal_driver_registered("MEM")) {
1867  RASTER_DEBUG(4, "Registering MEM driver");
1868  GDALRegister_MEM();
1869  *destroy_rtn_drv = 1;
1870  }
1871  drv = GDALGetDriverByName("MEM");
1872  if (NULL == drv) {
1873  rterror("rt_raster_to_gdal_mem: Could not load the MEM GDAL driver");
1874  return 0;
1875  }
1876  *rtn_drv = drv;
1877 
1878  /* unload driver from GDAL driver manager */
1879  if (*destroy_rtn_drv) {
1880  RASTER_DEBUG(4, "Deregistering MEM driver");
1881  GDALDeregisterDriver(drv);
1882  }
1883 
1884  width = rt_raster_get_width(raster);
1885  height = rt_raster_get_height(raster);
1886  ds = GDALCreate(
1887  drv, "",
1888  width, height,
1889  0, GDT_Byte, NULL
1890  );
1891  if (NULL == ds) {
1892  rterror("rt_raster_to_gdal_mem: Could not create a GDALDataset to convert into");
1893  return 0;
1894  }
1895 
1896  /* add geotransform */
1898  cplerr = GDALSetGeoTransform(ds, gt);
1899  if (cplerr != CE_None) {
1900  rterror("rt_raster_to_gdal_mem: Could not set geotransformation");
1901  GDALClose(ds);
1902  return 0;
1903  }
1904 
1905  /* set spatial reference */
1906  if (NULL != srs && strlen(srs)) {
1907  char *_srs = rt_util_gdal_convert_sr(srs, 0);
1908  if (_srs == NULL) {
1909  rterror("rt_raster_to_gdal_mem: Could not convert srs to GDAL accepted format");
1910  GDALClose(ds);
1911  return 0;
1912  }
1913 
1914  cplerr = GDALSetProjection(ds, _srs);
1915  CPLFree(_srs);
1916  if (cplerr != CE_None) {
1917  rterror("rt_raster_to_gdal_mem: Could not set projection");
1918  GDALClose(ds);
1919  return 0;
1920  }
1921  RASTER_DEBUGF(3, "Projection set to: %s", GDALGetProjectionRef(ds));
1922  }
1923 
1924  /* process bandNums */
1925  numBands = rt_raster_get_num_bands(raster);
1926  if (NULL != bandNums && count > 0) {
1927  for (i = 0; i < count; i++) {
1928  if (bandNums[i] >= numBands) {
1929  rterror("rt_raster_to_gdal_mem: The band index %d is invalid", bandNums[i]);
1930  GDALClose(ds);
1931  return 0;
1932  }
1933  }
1934  }
1935  else {
1936  count = numBands;
1937  bandNums = (uint32_t *) rtalloc(sizeof(uint32_t) * count);
1938  if (NULL == bandNums) {
1939  rterror("rt_raster_to_gdal_mem: Could not allocate memory for band indices");
1940  GDALClose(ds);
1941  return 0;
1942  }
1943  allocBandNums = 1;
1944  for (i = 0; i < count; i++) bandNums[i] = i;
1945  }
1946 
1947  /* process exclude_nodata_values */
1948  if (NULL == excludeNodataValues) {
1949  excludeNodataValues = (int *) rtalloc(sizeof(int) * count);
1950  if (NULL == excludeNodataValues) {
1951  rterror("rt_raster_to_gdal_mem: Could not allocate memory for NODATA flags");
1952  GDALClose(ds);
1953  return 0;
1954  }
1955  allocNodataValues = 1;
1956  for (i = 0; i < count; i++) excludeNodataValues[i] = 1;
1957  }
1958 
1959  /* add band(s) */
1960  for (i = 0; i < count; i++) {
1961  rtband = rt_raster_get_band(raster, bandNums[i]);
1962  if (NULL == rtband) {
1963  rterror("rt_raster_to_gdal_mem: Could not get requested band index %d", bandNums[i]);
1964  if (allocBandNums) rtdealloc(bandNums);
1965  if (allocNodataValues) rtdealloc(excludeNodataValues);
1966  GDALClose(ds);
1967  return 0;
1968  }
1969 
1970  pt = rt_band_get_pixtype(rtband);
1971  gdal_pt = rt_util_pixtype_to_gdal_datatype(pt);
1972  if (gdal_pt == GDT_Unknown)
1973  rtwarn("rt_raster_to_gdal_mem: Unknown pixel type for band");
1974 
1975  /*
1976  For all pixel types other than PT_8BSI, set pointer to start of data
1977  */
1978  if (pt != PT_8BSI) {
1979  pVoid = rt_band_get_data(rtband);
1980  RASTER_DEBUGF(4, "Band data is at pos %p", pVoid);
1981 
1982  pszDataPointer = (char *) rtalloc(20 * sizeof (char));
1983  sprintf(pszDataPointer, "%p", pVoid);
1984  RASTER_DEBUGF(4, "rt_raster_to_gdal_mem: szDatapointer is %p",
1985  pszDataPointer);
1986 
1987  if (strncasecmp(pszDataPointer, "0x", 2) == 0)
1988  sprintf(szGDALOption, "DATAPOINTER=%s", pszDataPointer);
1989  else
1990  sprintf(szGDALOption, "DATAPOINTER=0x%s", pszDataPointer);
1991 
1992  RASTER_DEBUG(3, "Storing info for GDAL MEM raster band");
1993 
1994  apszOptions[0] = szGDALOption;
1995  apszOptions[1] = NULL; /* pixel offset, not needed */
1996  apszOptions[2] = NULL; /* line offset, not needed */
1997  apszOptions[3] = NULL;
1998 
1999  /* free */
2000  rtdealloc(pszDataPointer);
2001 
2002  /* add band */
2003  if (GDALAddBand(ds, gdal_pt, apszOptions) == CE_Failure) {
2004  rterror("rt_raster_to_gdal_mem: Could not add GDAL raster band");
2005  if (allocBandNums) rtdealloc(bandNums);
2006  GDALClose(ds);
2007  return 0;
2008  }
2009  }
2010  /*
2011  PT_8BSI is special as GDAL has no equivalent pixel type.
2012  Must convert 8BSI to 16BSI so create basic band
2013  */
2014  else {
2015  /* add band */
2016  if (GDALAddBand(ds, gdal_pt, NULL) == CE_Failure) {
2017  rterror("rt_raster_to_gdal_mem: Could not add GDAL raster band");
2018  if (allocBandNums) rtdealloc(bandNums);
2019  if (allocNodataValues) rtdealloc(excludeNodataValues);
2020  GDALClose(ds);
2021  return 0;
2022  }
2023  }
2024 
2025  /* check band count */
2026  if (GDALGetRasterCount(ds) != i + 1) {
2027  rterror("rt_raster_to_gdal_mem: Error creating GDAL MEM raster band");
2028  if (allocBandNums) rtdealloc(bandNums);
2029  if (allocNodataValues) rtdealloc(excludeNodataValues);
2030  GDALClose(ds);
2031  return 0;
2032  }
2033 
2034  /* get new band */
2035  band = NULL;
2036  band = GDALGetRasterBand(ds, i + 1);
2037  if (NULL == band) {
2038  rterror("rt_raster_to_gdal_mem: Could not get GDAL band for additional processing");
2039  if (allocBandNums) rtdealloc(bandNums);
2040  if (allocNodataValues) rtdealloc(excludeNodataValues);
2041  GDALClose(ds);
2042  return 0;
2043  }
2044 
2045  /* PT_8BSI requires manual setting of pixels */
2046  if (pt == PT_8BSI) {
2047  uint32_t nXBlocks, nYBlocks;
2048  int nXBlockSize, nYBlockSize;
2049  uint32_t iXBlock, iYBlock;
2050  int nXValid, nYValid;
2051  int iX, iY;
2052  int iXMax, iYMax;
2053 
2054  int x, y, z;
2055  uint32_t valueslen = 0;
2056  int16_t *values = NULL;
2057  double value = 0.;
2058 
2059  /* this makes use of GDAL's "natural" blocks */
2060  GDALGetBlockSize(band, &nXBlockSize, &nYBlockSize);
2061  nXBlocks = (width + nXBlockSize - 1) / nXBlockSize;
2062  nYBlocks = (height + nYBlockSize - 1) / nYBlockSize;
2063  RASTER_DEBUGF(4, "(nXBlockSize, nYBlockSize) = (%d, %d)", nXBlockSize, nYBlockSize);
2064  RASTER_DEBUGF(4, "(nXBlocks, nYBlocks) = (%d, %d)", nXBlocks, nYBlocks);
2065 
2066  /* length is for the desired pixel type */
2067  valueslen = rt_pixtype_size(PT_16BSI) * nXBlockSize * nYBlockSize;
2068  values = rtalloc(valueslen);
2069  if (NULL == values) {
2070  rterror("rt_raster_to_gdal_mem: Could not allocate memory for GDAL band pixel values");
2071  if (allocBandNums) rtdealloc(bandNums);
2072  if (allocNodataValues) rtdealloc(excludeNodataValues);
2073  GDALClose(ds);
2074  return 0;
2075  }
2076 
2077  for (iYBlock = 0; iYBlock < nYBlocks; iYBlock++) {
2078  for (iXBlock = 0; iXBlock < nXBlocks; iXBlock++) {
2079  memset(values, 0, valueslen);
2080 
2081  x = iXBlock * nXBlockSize;
2082  y = iYBlock * nYBlockSize;
2083  RASTER_DEBUGF(4, "(iXBlock, iYBlock) = (%d, %d)", iXBlock, iYBlock);
2084  RASTER_DEBUGF(4, "(x, y) = (%d, %d)", x, y);
2085 
2086  /* valid block width */
2087  if ((iXBlock + 1) * nXBlockSize > width)
2088  nXValid = width - (iXBlock * nXBlockSize);
2089  else
2090  nXValid = nXBlockSize;
2091 
2092  /* valid block height */
2093  if ((iYBlock + 1) * nYBlockSize > height)
2094  nYValid = height - (iYBlock * nYBlockSize);
2095  else
2096  nYValid = nYBlockSize;
2097 
2098  RASTER_DEBUGF(4, "(nXValid, nYValid) = (%d, %d)", nXValid, nYValid);
2099 
2100  /* convert 8BSI values to 16BSI */
2101  z = 0;
2102  iYMax = y + nYValid;
2103  iXMax = x + nXValid;
2104  for (iY = y; iY < iYMax; iY++) {
2105  for (iX = x; iX < iXMax; iX++) {
2106  if (rt_band_get_pixel(rtband, iX, iY, &value, NULL) != ES_NONE) {
2107  rterror("rt_raster_to_gdal_mem: Could not get pixel value to convert from 8BSI to 16BSI");
2108  rtdealloc(values);
2109  if (allocBandNums) rtdealloc(bandNums);
2110  if (allocNodataValues) rtdealloc(excludeNodataValues);
2111  GDALClose(ds);
2112  return 0;
2113  }
2114 
2115  values[z++] = rt_util_clamp_to_16BSI(value);
2116  }
2117  }
2118 
2119  /* burn values */
2120  if (GDALRasterIO(
2121  band, GF_Write,
2122  x, y,
2123  nXValid, nYValid,
2124  values, nXValid, nYValid,
2125  gdal_pt,
2126  0, 0
2127  ) != CE_None) {
2128  rterror("rt_raster_to_gdal_mem: Could not write converted 8BSI to 16BSI values to GDAL band");
2129  rtdealloc(values);
2130  if (allocBandNums) rtdealloc(bandNums);
2131  if (allocNodataValues) rtdealloc(excludeNodataValues);
2132  GDALClose(ds);
2133  return 0;
2134  }
2135  }
2136  }
2137 
2138  rtdealloc(values);
2139  }
2140 
2141  /* Add nodata value for band */
2142  if (rt_band_get_hasnodata_flag(rtband) != FALSE && excludeNodataValues[i]) {
2143  rt_band_get_nodata(rtband, &nodata);
2144  if (GDALSetRasterNoDataValue(band, nodata) != CE_None)
2145  rtwarn("rt_raster_to_gdal_mem: Could not set nodata value for band");
2146  RASTER_DEBUGF(3, "nodata value set to %f", GDALGetRasterNoDataValue(band, NULL));
2147  }
2148 
2149 #if POSTGIS_DEBUG_LEVEL > 3
2150  {
2151  GDALRasterBandH _grb = NULL;
2152  double _min;
2153  double _max;
2154  double _mean;
2155  double _stddev;
2156 
2157  _grb = GDALGetRasterBand(ds, i + 1);
2158  GDALComputeRasterStatistics(_grb, FALSE, &_min, &_max, &_mean, &_stddev, NULL, NULL);
2159  RASTER_DEBUGF(4, "GDAL Band %d stats: %f, %f, %f, %f", i + 1, _min, _max, _mean, _stddev);
2160  }
2161 #endif
2162 
2163  }
2164 
2165  /* necessary??? */
2166  GDALFlushCache(ds);
2167 
2168  if (allocBandNums) rtdealloc(bandNums);
2169  if (allocNodataValues) rtdealloc(excludeNodataValues);
2170 
2171  return ds;
2172 }
tuple gt
Definition: window.py:77
void rterror(const char *fmt,...)
Wrappers used for reporting errors and info.
Definition: rt_context.c:199
tuple band
Definition: ovdump.py:57
void * rtalloc(size_t size)
Wrappers used for managing memory.
Definition: rt_context.c:171
rt_pixtype
Definition: librtcore.h:185
uint16_t rt_raster_get_num_bands(rt_raster raster)
Definition: rt_raster.c:372
rt_errorstate rt_band_get_nodata(rt_band band, double *nodata)
Get NODATA value.
Definition: rt_band.c:1597
unsigned int uint32_t
Definition: uthash.h:78
GDALDataType rt_util_pixtype_to_gdal_datatype(rt_pixtype pt)
Convert rt_pixtype to GDALDataType.
Definition: rt_util.c:118
int rt_util_gdal_driver_registered(const char *drv)
Definition: rt_util.c:353
void rtwarn(const char *fmt,...)
Definition: rt_context.c:224
rt_errorstate rt_band_get_pixel(rt_band band, int x, int y, double *value, int *nodata)
Get pixel value.
Definition: rt_band.c:1088
void rt_raster_get_geotransform_matrix(rt_raster raster, double *gt)
Get 6-element array of raster geotransform matrix.
Definition: rt_raster.c:706
int rt_band_get_hasnodata_flag(rt_band band)
Get hasnodata flag value.
Definition: rt_band.c:541
int count
Definition: genraster.py:56
tuple ds
Definition: pixval.py:66
void * rt_band_get_data(rt_band band)
Get pointer to raster band data.
Definition: rt_band.c:302
#define RASTER_DEBUGF(level, msg,...)
Definition: librtcore.h:299
tuple x
Definition: pixval.py:53
int rt_pixtype_size(rt_pixtype pixtype)
Return size in bytes of a value in the given pixtype.
Definition: rt_pixel.c:39
void rtdealloc(void *mem)
Definition: rt_context.c:186
#define FALSE
Definition: dbfopen.c:168
#define RASTER_DEBUG(level, msg)
Definition: librtcore.h:295
char * rt_util_gdal_convert_sr(const char *srs, int proj4)
Definition: rt_util.c:212
rt_pixtype rt_band_get_pixtype(rt_band band)
Return pixeltype of this band.
Definition: rt_band.c:498
uint16_t rt_raster_get_height(rt_raster raster)
Definition: rt_raster.c:129
uint16_t rt_raster_get_width(rt_raster raster)
Definition: rt_raster.c:121
tuple y
Definition: pixval.py:54
rt_band rt_raster_get_band(rt_raster raster, int n)
Return Nth band, or NULL if unavailable.
Definition: rt_raster.c:381
int16_t rt_util_clamp_to_16BSI(double value)
Definition: rt_util.c:59

Here is the call graph for this function:

Here is the caller graph for this function: