PostGIS 3.0.6dev-r@@SVN_REVISION@@
Loading...
Searching...
No Matches

◆ convert_raster()

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

Definition at line 1541 of file raster2pgsql.c.

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

References _, append_stringbuffer(), rasterinfo_t::bandtype, buffer(), 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, rasterinfo_t::gt, raster_loader_config::hasnodata, rasterinfo_t::hasnodata, insert_records(), stringbuffer_t::length, MAXTILESIZE, 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, 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().

Here is the call graph for this function:
Here is the caller graph for this function: