PostGIS  2.4.9dev-r@@SVN_REVISION@@

◆ convert_raster()

static int convert_raster ( int  idx,
RTLOADERCFG config,
RASTERINFO info,
STRINGBUFFER tileset,
STRINGBUFFER buffer 
)
static

Definition at line 1563 of file raster2pgsql.c.

References _, append_stringbuffer(), ovdump::band, rasterinfo_t::bandtype, calc_tile_size(), raster_loader_config::copy_statements, rasterinfo_t::dim, FALSE, raster_loader_config::file_column, raster_loader_config::file_column_name, rasterinfo_t::gdalbandtype, window::gt, rasterinfo_t::gt, raster_loader_config::hasnodata, rasterinfo_t::hasnodata, if(), insert_records(), stringbuffer_t::length, MAXTILESIZE, pixval::nband, raster_loader_config::nband, rasterinfo_t::nband, raster_loader_config::nband_count, rasterinfo_t::nband_count, raster_loader_config::nodataval, rasterinfo_t::nodataval, raster_loader_config::out_srid, raster_loader_config::outdb, raster_loader_config::pad_tile, PT_END, rtpixdump::rast, raster_loader_config::raster_column, raster_destroy(), rt_band_check_is_nodata(), rt_band_destroy(), rt_band_new_offline(), raster_loader_config::rt_file, raster_loader_config::rt_filename, rt_pixtype_size(), rt_raster_add_band(), rt_raster_from_gdal_dataset(), rt_raster_get_band(), rt_raster_get_num_bands(), rt_raster_new(), rt_raster_set_geotransform_matrix(), rt_raster_set_srid(), rt_raster_to_hexwkb(), rt_util_gdal_datatype_to_pixtype(), rtalloc(), rtdealloc_stringbuffer(), rterror(), rtinfo(), rtwarn(), raster_loader_config::schema, raster_loader_config::skip_nodataval_check, raster_loader_config::srid, rasterinfo_t::srid, SRID_UNKNOWN, rasterinfo_t::srs, raster_loader_config::table, raster_loader_config::tile_size, and rasterinfo_t::tile_size.

Referenced by process_rasters().

1563  {
1564  GDALDatasetH hdsSrc;
1565  GDALRasterBandH hbandSrc;
1566  int nband = 0;
1567  int i = 0;
1568  int ntiles[2] = {1, 1};
1569  int _tile_size[2] = {0, 0};
1570  int xtile = 0;
1571  int ytile = 0;
1572  double gt[6] = {0.};
1573  const char* pszProjectionRef = NULL;
1574  int tilesize = 0;
1575 
1576  rt_raster rast = NULL;
1577  int numbands = 0;
1578  rt_band band = NULL;
1579  char *hex;
1580  uint32_t hexlen = 0;
1581 
1582  info->srid = config->srid;
1583 
1584  hdsSrc = GDALOpenShared(config->rt_file[idx], GA_ReadOnly);
1585  if (hdsSrc == NULL) {
1586  rterror(_("convert_raster: Could not open raster: %s"), config->rt_file[idx]);
1587  return 0;
1588  }
1589 
1590  nband = GDALGetRasterCount(hdsSrc);
1591  if (!nband) {
1592  rterror(_("convert_raster: No bands found in raster: %s"), config->rt_file[idx]);
1593  GDALClose(hdsSrc);
1594  return 0;
1595  }
1596 
1597  /* check that bands specified are available */
1598  for (i = 0; i < config->nband_count; i++) {
1599  if (config->nband[i] > nband) {
1600  rterror(_("convert_raster: Band %d not found in raster: %s"), config->nband[i], config->rt_file[idx]);
1601  GDALClose(hdsSrc);
1602  return 0;
1603  }
1604  }
1605 
1606  /* record srs */
1607  pszProjectionRef = GDALGetProjectionRef(hdsSrc);
1608  if (pszProjectionRef != NULL && pszProjectionRef[0] != '\0') {
1609  info->srs = rtalloc(sizeof(char) * (strlen(pszProjectionRef) + 1));
1610  if (info->srs == NULL) {
1611  rterror(_("convert_raster: Could not allocate memory for storing SRS"));
1612  GDALClose(hdsSrc);
1613  return 0;
1614  }
1615  strcpy(info->srs, pszProjectionRef);
1616 
1617  if (info->srid == SRID_UNKNOWN) {
1618  OGRSpatialReferenceH hSRS = OSRNewSpatialReference(NULL);
1619  if (OSRSetFromUserInput(hSRS, pszProjectionRef) == OGRERR_NONE) {
1620  const char* pszAuthorityName = OSRGetAuthorityName(hSRS, NULL);
1621  const char* pszAuthorityCode = OSRGetAuthorityCode(hSRS, NULL);
1622  if (
1623  pszAuthorityName != NULL &&
1624  strcmp(pszAuthorityName, "EPSG") == 0 &&
1625  pszAuthorityCode != NULL
1626  ) {
1627  info->srid = atoi(pszAuthorityCode);
1628  }
1629  }
1630  OSRDestroySpatialReference(hSRS);
1631  }
1632  }
1633 
1634  if ( info->srid == SRID_UNKNOWN && config->out_srid != SRID_UNKNOWN ) {
1635  rterror(_("convert_raster: could not determine source srid, cannot transform to target srid %d"), config->out_srid);
1636  GDALClose(hdsSrc);
1637  return 0;
1638  }
1639 
1640  /* record geotransform matrix */
1641  if (GDALGetGeoTransform(hdsSrc, info->gt) != CE_None) {
1642  rtinfo(_("Using default geotransform matrix (0, 1, 0, 0, 0, -1) for raster: %s"), config->rt_file[idx]);
1643  info->gt[0] = 0;
1644  info->gt[1] = 1;
1645  info->gt[2] = 0;
1646  info->gt[3] = 0;
1647  info->gt[4] = 0;
1648  info->gt[5] = -1;
1649  }
1650  memcpy(gt, info->gt, sizeof(double) * 6);
1651 
1652  /* record # of bands */
1653  /* user-specified bands */
1654  if (config->nband_count > 0) {
1655  info->nband_count = config->nband_count;
1656  info->nband = rtalloc(sizeof(int) * info->nband_count);
1657  if (info->nband == NULL) {
1658  rterror(_("convert_raster: Could not allocate memory for storing band indices"));
1659  GDALClose(hdsSrc);
1660  return 0;
1661  }
1662  memcpy(info->nband, config->nband, sizeof(int) * info->nband_count);
1663  }
1664  /* all bands */
1665  else {
1666  info->nband_count = nband;
1667  info->nband = rtalloc(sizeof(int) * info->nband_count);
1668  if (info->nband == NULL) {
1669  rterror(_("convert_raster: Could not allocate memory for storing band indices"));
1670  GDALClose(hdsSrc);
1671  return 0;
1672  }
1673  for (i = 0; i < info->nband_count; i++)
1674  info->nband[i] = i + 1;
1675  }
1676 
1677  /* initialize parameters dependent on nband */
1678  info->gdalbandtype = rtalloc(sizeof(GDALDataType) * info->nband_count);
1679  if (info->gdalbandtype == NULL) {
1680  rterror(_("convert_raster: Could not allocate memory for storing GDAL data type"));
1681  GDALClose(hdsSrc);
1682  return 0;
1683  }
1684  info->bandtype = rtalloc(sizeof(rt_pixtype) * info->nband_count);
1685  if (info->bandtype == NULL) {
1686  rterror(_("convert_raster: Could not allocate memory for storing pixel type"));
1687  GDALClose(hdsSrc);
1688  return 0;
1689  }
1690  info->hasnodata = rtalloc(sizeof(int) * info->nband_count);
1691  if (info->hasnodata == NULL) {
1692  rterror(_("convert_raster: Could not allocate memory for storing hasnodata flag"));
1693  GDALClose(hdsSrc);
1694  return 0;
1695  }
1696  info->nodataval = rtalloc(sizeof(double) * info->nband_count);
1697  if (info->nodataval == NULL) {
1698  rterror(_("convert_raster: Could not allocate memory for storing nodata value"));
1699  GDALClose(hdsSrc);
1700  return 0;
1701  }
1702  memset(info->gdalbandtype, GDT_Unknown, sizeof(GDALDataType) * info->nband_count);
1703  memset(info->bandtype, PT_END, sizeof(rt_pixtype) * info->nband_count);
1704  memset(info->hasnodata, 0, sizeof(int) * info->nband_count);
1705  memset(info->nodataval, 0, sizeof(double) * info->nband_count);
1706 
1707  /* dimensions of raster */
1708  info->dim[0] = GDALGetRasterXSize(hdsSrc);
1709  info->dim[1] = GDALGetRasterYSize(hdsSrc);
1710 
1711  /* tile size is "auto" */
1712  if (
1713  config->tile_size[0] == -1 &&
1714  config->tile_size[1] == -1
1715  ) {
1717  info->dim[0], info->dim[1],
1718  &(config->tile_size[0]), &(config->tile_size[1])
1719  );
1720 
1721  rtinfo(_("Using computed tile size: %dx%d"), config->tile_size[0], config->tile_size[1]);
1722  }
1723 
1724  /* decide on tile size */
1725  if (!config->tile_size[0])
1726  info->tile_size[0] = info->dim[0];
1727  else
1728  info->tile_size[0] = config->tile_size[0];
1729  if (!config->tile_size[1])
1730  info->tile_size[1] = info->dim[1];
1731  else
1732  info->tile_size[1] = config->tile_size[1];
1733 
1734  /* number of tiles */
1735  if (info->tile_size[0] != info->dim[0])
1736  ntiles[0] = (info->dim[0] + info->tile_size[0] - 1) / info->tile_size[0];
1737  if (info->tile_size[1] != info->dim[1])
1738  ntiles[1] = (info->dim[1] + info->tile_size[1] - 1) / info->tile_size[1];
1739 
1740  /* estimate size of 1 tile */
1741  tilesize = info->tile_size[0] * info->tile_size[1];
1742 
1743  /* go through bands for attributes */
1744  for (i = 0; i < info->nband_count; i++) {
1745  hbandSrc = GDALGetRasterBand(hdsSrc, info->nband[i]);
1746 
1747  /* datatype */
1748  info->gdalbandtype[i] = GDALGetRasterDataType(hbandSrc);
1749 
1750  /* complex data type? */
1751  if (GDALDataTypeIsComplex(info->gdalbandtype[i])) {
1752  rterror(_("convert_raster: The pixel type of band %d is a complex data type. PostGIS raster does not support complex data types"), i + 1);
1753  GDALClose(hdsSrc);
1754  return 0;
1755  }
1756 
1757  /* convert data type to that of postgis raster */
1759 
1760  /* hasnodata and nodataval */
1761  info->nodataval[i] = GDALGetRasterNoDataValue(hbandSrc, &(info->hasnodata[i]));
1762  if (!info->hasnodata[i]) {
1763  /* does NOT have nodata value, but user-specified */
1764  if (config->hasnodata) {
1765  info->hasnodata[i] = 1;
1766  info->nodataval[i] = config->nodataval;
1767  }
1768  else
1769  info->nodataval[i] = 0;
1770  }
1771 
1772  /* update estimated size of 1 tile */
1773  tilesize *= rt_pixtype_size(info->bandtype[i]);
1774  }
1775 
1776  /* roughly estimate size of one tile and all bands */
1777  tilesize *= 1.1;
1778  if (tilesize > MAXTILESIZE)
1779  rtwarn(_("The size of each output tile may exceed 1 GB. Use -t to specify a reasonable tile size"));
1780 
1781  /* out-db raster */
1782  if (config->outdb) {
1783  GDALClose(hdsSrc);
1784 
1785  /* each tile is a raster */
1786  for (ytile = 0; ytile < ntiles[1]; ytile++) {
1787  /* edge y tile */
1788  if (!config->pad_tile && ntiles[1] > 1 && (ytile + 1) == ntiles[1])
1789  _tile_size[1] = info->dim[1] - (ytile * info->tile_size[1]);
1790  else
1791  _tile_size[1] = info->tile_size[1];
1792 
1793  for (xtile = 0; xtile < ntiles[0]; xtile++) {
1794 
1795  /* edge x tile */
1796  if (!config->pad_tile && ntiles[0] > 1 && (xtile + 1) == ntiles[0])
1797  _tile_size[0] = info->dim[0] - (xtile * info->tile_size[0]);
1798  else
1799  _tile_size[0] = info->tile_size[0];
1800 
1801  /* compute tile's upper-left corner */
1802  GDALApplyGeoTransform(
1803  info->gt,
1804  xtile * info->tile_size[0], ytile * info->tile_size[1],
1805  &(gt[0]), &(gt[3])
1806  );
1807 
1808  /* create raster object */
1809  rast = rt_raster_new(_tile_size[0], _tile_size[1]);
1810  if (rast == NULL) {
1811  rterror(_("convert_raster: Could not create raster"));
1812  return 0;
1813  }
1814 
1815  /* set raster attributes */
1816  rt_raster_set_srid(rast, info->srid);
1818 
1819  /* add bands */
1820  for (i = 0; i < info->nband_count; i++) {
1821  band = rt_band_new_offline(
1822  _tile_size[0], _tile_size[1],
1823  info->bandtype[i],
1824  info->hasnodata[i], info->nodataval[i],
1825  info->nband[i] - 1,
1826  config->rt_file[idx]
1827  );
1828  if (band == NULL) {
1829  rterror(_("convert_raster: Could not create offline band"));
1830  raster_destroy(rast);
1831  return 0;
1832  }
1833 
1834  /* add band to raster */
1835  if (rt_raster_add_band(rast, band, rt_raster_get_num_bands(rast)) == -1) {
1836  rterror(_("convert_raster: Could not add offlineband to raster"));
1837  rt_band_destroy(band);
1838  raster_destroy(rast);
1839  return 0;
1840  }
1841 
1842  /* inspect each band of raster where band is NODATA */
1843  if (!config->skip_nodataval_check)
1845  }
1846 
1847  /* convert rt_raster to hexwkb */
1848  hex = rt_raster_to_hexwkb(rast, FALSE, &hexlen);
1849  raster_destroy(rast);
1850 
1851  if (hex == NULL) {
1852  rterror(_("convert_raster: Could not convert PostGIS raster to hex WKB"));
1853  return 0;
1854  }
1855 
1856  /* add hexwkb to tileset */
1857  append_stringbuffer(tileset, hex);
1858 
1859  /* flush if tileset gets too big */
1860  if (tileset->length > 10) {
1861  if (!insert_records(
1862  config->schema, config->table, config->raster_column,
1863  (config->file_column ? config->rt_filename[idx] : NULL), config->file_column_name,
1864  config->copy_statements, config->out_srid,
1865  tileset, buffer
1866  )) {
1867  rterror(_("convert_raster: Could not convert raster tiles into INSERT or COPY statements"));
1868  return 0;
1869  }
1870 
1871  rtdealloc_stringbuffer(tileset, 0);
1872  }
1873  }
1874  }
1875  }
1876  /* in-db raster */
1877  else {
1878  VRTDatasetH hdsDst;
1879  VRTSourcedRasterBandH hbandDst;
1880 
1881  /* each tile is a VRT with constraints set for just the data required for the tile */
1882  for (ytile = 0; ytile < ntiles[1]; ytile++) {
1883 
1884  /* edge y tile */
1885  if (!config->pad_tile && ntiles[1] > 1 && (ytile + 1) == ntiles[1])
1886  _tile_size[1] = info->dim[1] - (ytile * info->tile_size[1]);
1887  else
1888  _tile_size[1] = info->tile_size[1];
1889 
1890  for (xtile = 0; xtile < ntiles[0]; xtile++) {
1891  /*
1892  char fn[100];
1893  sprintf(fn, "/tmp/tile%d.vrt", (ytile * ntiles[0]) + xtile);
1894  */
1895 
1896  /* edge x tile */
1897  if (!config->pad_tile && ntiles[0] > 1 && (xtile + 1) == ntiles[0])
1898  _tile_size[0] = info->dim[0] - (xtile * info->tile_size[0]);
1899  else
1900  _tile_size[0] = info->tile_size[0];
1901 
1902  /* compute tile's upper-left corner */
1903  GDALApplyGeoTransform(
1904  info->gt,
1905  xtile * info->tile_size[0], ytile * info->tile_size[1],
1906  &(gt[0]), &(gt[3])
1907  );
1908  /*
1909  rtinfo(_("tile (%d, %d) gt = (%f, %f, %f, %f, %f, %f)"),
1910  xtile, ytile,
1911  gt[0], gt[1], gt[2], gt[3], gt[4], gt[5]
1912  );
1913  */
1914 
1915  /* create VRT dataset */
1916  hdsDst = VRTCreate(_tile_size[0], _tile_size[1]);
1917  /*
1918  GDALSetDescription(hdsDst, fn);
1919  */
1920  GDALSetProjection(hdsDst, info->srs);
1921  GDALSetGeoTransform(hdsDst, gt);
1922 
1923  /* add bands as simple sources */
1924  for (i = 0; i < info->nband_count; i++) {
1925  GDALAddBand(hdsDst, info->gdalbandtype[i], NULL);
1926  hbandDst = (VRTSourcedRasterBandH) GDALGetRasterBand(hdsDst, i + 1);
1927 
1928  if (info->hasnodata[i])
1929  GDALSetRasterNoDataValue(hbandDst, info->nodataval[i]);
1930 
1931  VRTAddSimpleSource(
1932  hbandDst, GDALGetRasterBand(hdsSrc, info->nband[i]),
1933  xtile * info->tile_size[0], ytile * info->tile_size[1],
1934  _tile_size[0], _tile_size[1],
1935  0, 0,
1936  _tile_size[0], _tile_size[1],
1937  "near", VRT_NODATA_UNSET
1938  );
1939  }
1940 
1941  /* make sure VRT reflects all changes */
1942  VRTFlushCache(hdsDst);
1943 
1944  /* convert VRT dataset to rt_raster */
1945  rast = rt_raster_from_gdal_dataset(hdsDst);
1946  if (rast == NULL) {
1947  rterror(_("convert_raster: Could not convert VRT dataset to PostGIS raster"));
1948  GDALClose(hdsDst);
1949  return 0;
1950  }
1951 
1952  /* set srid if provided */
1953  rt_raster_set_srid(rast, info->srid);
1954 
1955  /* inspect each band of raster where band is NODATA */
1956  numbands = rt_raster_get_num_bands(rast);
1957  for (i = 0; i < numbands; i++) {
1958  band = rt_raster_get_band(rast, i);
1959  if (band != NULL && !config->skip_nodataval_check)
1961  }
1962 
1963  /* convert rt_raster to hexwkb */
1964  hex = rt_raster_to_hexwkb(rast, FALSE, &hexlen);
1965  raster_destroy(rast);
1966 
1967  if (hex == NULL) {
1968  rterror(_("convert_raster: Could not convert PostGIS raster to hex WKB"));
1969  GDALClose(hdsDst);
1970  return 0;
1971  }
1972 
1973  /* add hexwkb to tileset */
1974  append_stringbuffer(tileset, hex);
1975 
1976  GDALClose(hdsDst);
1977 
1978  /* flush if tileset gets too big */
1979  if (tileset->length > 10) {
1980  if (!insert_records(
1981  config->schema, config->table, config->raster_column,
1982  (config->file_column ? config->rt_filename[idx] : NULL), config->file_column_name,
1983  config->copy_statements, config->out_srid,
1984  tileset, buffer
1985  )) {
1986  rterror(_("convert_raster: Could not convert raster tiles into INSERT or COPY statements"));
1987  GDALClose(hdsSrc);
1988  return 0;
1989  }
1990 
1991  rtdealloc_stringbuffer(tileset, 0);
1992  }
1993  }
1994  }
1995 
1996  GDALClose(hdsSrc);
1997  }
1998 
1999  return 1;
2000 }
GDALDataType * gdalbandtype
Definition: raster2pgsql.h:179
int rt_raster_get_num_bands(rt_raster raster)
Definition: rt_raster.c:372
#define MAXTILESIZE
Definition: raster2pgsql.h:68
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
rt_pixtype * bandtype
Definition: raster2pgsql.h:180
static void rtdealloc_stringbuffer(STRINGBUFFER *buffer, int freebuffer)
Definition: raster2pgsql.c:774
#define _(String)
Definition: shpcommon.h:24
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
uint32_t length
Definition: raster2pgsql.h:196
static int append_stringbuffer(STRINGBUFFER *buffer, const char *str)
Definition: raster2pgsql.c:807
gt
Definition: window.py:77
int tile_size[2]
Definition: raster2pgsql.h:191
rt_pixtype
Definition: librtcore.h:185
void rt_band_destroy(rt_band band)
Destroy a raster band.
Definition: rt_band.c:242
unsigned int uint32_t
Definition: uthash.h:78
static int insert_records(const char *schema, const char *table, const char *column, const char *filename, const char *file_column_name, int copy_statements, int out_srid, STRINGBUFFER *tileset, STRINGBUFFER *buffer)
Definition: raster2pgsql.c:877
nband
Definition: pixval.py:52
void rtwarn(const char *fmt,...)
Definition: rt_context.c:224
char * rt_raster_to_hexwkb(rt_raster raster, int outasin, uint32_t *hexwkbsize)
Return this raster in HEXWKB form (null-terminated hex)
Definition: rt_wkb.c:669
#define SRID_UNKNOWN
Unknown SRID value.
Definition: liblwgeom.h:188
uint32_t dim[2]
Definition: raster2pgsql.h:172
rt_band rt_raster_get_band(rt_raster raster, int bandNum)
Return Nth band, or NULL if unavailable.
Definition: rt_raster.c:381
rt_pixtype rt_util_gdal_datatype_to_pixtype(GDALDataType gdt)
Convert GDALDataType to rt_pixtype.
Definition: rt_util.c:153
void rtinfo(const char *fmt,...)
Definition: rt_context.c:211
int rt_band_check_is_nodata(rt_band band)
Returns TRUE if the band is only nodata values.
Definition: rt_band.c:1619
void rt_raster_set_srid(rt_raster raster, int32_t srid)
Set raster&#39;s SRID.
Definition: rt_raster.c:363
rt_raster rt_raster_new(uint32_t width, uint32_t height)
Construct a raster with given dimensions.
Definition: rt_raster.c:48
static void calc_tile_size(int dimX, int dimY, int *tileX, int *tileY)
Definition: raster2pgsql.c:450
int rt_pixtype_size(rt_pixtype pixtype)
Return size in bytes of a value in the given pixtype.
Definition: rt_pixel.c:39
#define FALSE
Definition: dbfopen.c:168
static void raster_destroy(rt_raster raster)
Definition: raster2pgsql.c:76
int rt_raster_add_band(rt_raster raster, rt_band band, int index)
Add band data to a raster.
Definition: rt_raster.c:405
rt_raster rt_raster_from_gdal_dataset(GDALDatasetH ds)
Return a raster from a GDAL dataset.
Definition: rt_raster.c:2165
double * nodataval
Definition: raster2pgsql.h:185
double gt[6]
Definition: raster2pgsql.h:188
if(!(yy_init))
rt_band rt_band_new_offline(uint16_t width, uint16_t height, rt_pixtype pixtype, uint32_t hasnodata, double nodataval, uint8_t bandNum, const char *path)
Create an out-db rt_band.
Definition: rt_band.c:119
Here is the call graph for this function:
Here is the caller graph for this function: