PostGIS  2.5.0beta2dev-r@@SVN_REVISION@@

◆ convert_raster()

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

Definition at line 1566 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().

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