PostGIS  2.1.10dev-r@@SVN_REVISION@@
static int convert_raster ( int  idx,
RTLOADERCFG config,
RASTERINFO info,
STRINGBUFFER tileset,
STRINGBUFFER buffer 
)
static

Definition at line 1546 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::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(), 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().

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

Here is the call graph for this function:

Here is the caller graph for this function: