PostGIS  2.4.9dev-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 1809 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().

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