PostGIS 3.7.0dev-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 1553 of file raster2pgsql.c.

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