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

◆ rt_raster_to_gdal_mem()

GDALDatasetH rt_raster_to_gdal_mem ( rt_raster  raster,
const char *  srs,
uint32_t *  bandNums,
int *  excludeNodataValues,
int  count,
GDALDriverH *  rtn_drv,
int *  destroy_rtn_drv 
)

Return GDAL dataset using GDAL MEM driver from raster.

Parameters
raster: raster to convert to GDAL MEM
srs: the raster's coordinate system in OGC WKT
bandNums: array of band numbers to extract from raster and include in the GDAL dataset (0 based)
excludeNodataValues: array of zero, nonzero where if non-zero, ignore nodata values for the band to check for pixels with value
count: number of elements in bandNums and exclude_nodata_values
rtn_drv: is set to the GDAL driver object
destroy_rtn_drv: if non-zero, caller must destroy the MEM driver
Returns
GDAL dataset using GDAL MEM driver
Parameters
raster: raster to convert to GDAL MEM
srs: the raster's coordinate system in OGC WKT
bandNums: array of band numbers to extract from raster and include in the GDAL dataset (0 based)
excludeNodataValues: array of zero, nonzero where if non-zero, ignore nodata values for the band
count: number of elements in bandNums
rtn_drv: is set to the GDAL driver object
destroy_rtn_drv: if non-zero, caller must destroy the MEM driver
Returns
GDAL dataset using GDAL MEM driver

Definition at line 1813 of file rt_raster.c.

1820 {
1821 GDALDriverH drv = NULL;
1822 GDALDatasetH ds = NULL;
1823 double gt[6] = {0.0};
1824 CPLErr cplerr;
1825 GDALDataType gdal_pt = GDT_Unknown;
1826 GDALRasterBandH band;
1827 void *pVoid;
1828 char *pszDataPointer;
1829 char szGDALOption[50];
1830 char *apszOptions[4];
1831 double nodata = 0.0;
1832 int allocBandNums = 0;
1833 int allocNodataValues = 0;
1834
1835 int i;
1836 uint32_t numBands;
1837 uint32_t width = 0;
1838 uint32_t height = 0;
1839 rt_band rtband = NULL;
1840 rt_pixtype pt = PT_END;
1841
1842 assert(NULL != raster);
1843 assert(NULL != rtn_drv);
1844 assert(NULL != destroy_rtn_drv);
1845
1846 *destroy_rtn_drv = 0;
1847
1848 /* store raster in GDAL MEM raster */
1849 if (!rt_util_gdal_driver_registered("MEM")) {
1850 RASTER_DEBUG(4, "Registering MEM driver");
1851 GDALRegister_MEM();
1852 *destroy_rtn_drv = 1;
1853 }
1854 drv = GDALGetDriverByName("MEM");
1855 if (NULL == drv) {
1856 rterror("rt_raster_to_gdal_mem: Could not load the MEM GDAL driver");
1857 return 0;
1858 }
1859 *rtn_drv = drv;
1860
1861 /* unload driver from GDAL driver manager */
1862 if (*destroy_rtn_drv) {
1863 RASTER_DEBUG(4, "Deregistering MEM driver");
1864 GDALDeregisterDriver(drv);
1865 }
1866
1867 width = rt_raster_get_width(raster);
1868 height = rt_raster_get_height(raster);
1869 ds = GDALCreate(
1870 drv, "",
1871 width, height,
1872 0, GDT_Byte, NULL
1873 );
1874 if (NULL == ds) {
1875 rterror("rt_raster_to_gdal_mem: Could not create a GDALDataset to convert into");
1876 return 0;
1877 }
1878
1879 /* add geotransform */
1880 rt_raster_get_geotransform_matrix(raster, gt);
1881 cplerr = GDALSetGeoTransform(ds, gt);
1882 if (cplerr != CE_None) {
1883 rterror("rt_raster_to_gdal_mem: Could not set geotransformation");
1884 GDALClose(ds);
1885 return 0;
1886 }
1887
1888 /* set spatial reference */
1889 if (NULL != srs && strlen(srs)) {
1890 char *_srs = rt_util_gdal_convert_sr(srs, 0);
1891 if (_srs == NULL) {
1892 rterror("rt_raster_to_gdal_mem: Could not convert srs to GDAL accepted format");
1893 GDALClose(ds);
1894 return 0;
1895 }
1896
1897 cplerr = GDALSetProjection(ds, _srs);
1898 CPLFree(_srs);
1899 if (cplerr != CE_None) {
1900 rterror("rt_raster_to_gdal_mem: Could not set projection");
1901 GDALClose(ds);
1902 return 0;
1903 }
1904 RASTER_DEBUGF(3, "Projection set to: %s", GDALGetProjectionRef(ds));
1905 }
1906
1907 /* process bandNums */
1908 numBands = rt_raster_get_num_bands(raster);
1909 if (NULL != bandNums && count > 0) {
1910 for (i = 0; i < count; i++) {
1911 if (bandNums[i] >= numBands) {
1912 rterror("rt_raster_to_gdal_mem: The band index %d is invalid", bandNums[i]);
1913 GDALClose(ds);
1914 return 0;
1915 }
1916 }
1917 }
1918 else {
1919 count = numBands;
1920 bandNums = (uint32_t *) rtalloc(sizeof(uint32_t) * count);
1921 if (NULL == bandNums) {
1922 rterror("rt_raster_to_gdal_mem: Could not allocate memory for band indices");
1923 GDALClose(ds);
1924 return 0;
1925 }
1926 allocBandNums = 1;
1927 for (i = 0; i < count; i++) bandNums[i] = i;
1928 }
1929
1930 /* process exclude_nodata_values */
1931 if (NULL == excludeNodataValues) {
1932 excludeNodataValues = (int *) rtalloc(sizeof(int) * count);
1933 if (NULL == excludeNodataValues) {
1934 rterror("rt_raster_to_gdal_mem: Could not allocate memory for NODATA flags");
1935 GDALClose(ds);
1936 return 0;
1937 }
1938 allocNodataValues = 1;
1939 for (i = 0; i < count; i++) excludeNodataValues[i] = 1;
1940 }
1941
1942 /* add band(s) */
1943 for (i = 0; i < count; i++) {
1944 rtband = rt_raster_get_band(raster, bandNums[i]);
1945 if (NULL == rtband) {
1946 rterror("rt_raster_to_gdal_mem: Could not get requested band index %d", bandNums[i]);
1947 if (allocBandNums) rtdealloc(bandNums);
1948 if (allocNodataValues) rtdealloc(excludeNodataValues);
1949 GDALClose(ds);
1950 return 0;
1951 }
1952
1953 pt = rt_band_get_pixtype(rtband);
1954 gdal_pt = rt_util_pixtype_to_gdal_datatype(pt);
1955 if (gdal_pt == GDT_Unknown)
1956 rtwarn("rt_raster_to_gdal_mem: Unknown pixel type for band");
1957
1958 /*
1959 For all pixel types other than PT_8BSI, set pointer to start of data
1960 */
1961#if POSTGIS_GDAL_VERSION < 30700
1962 if (pt != PT_8BSI) {
1963#endif
1964 pVoid = rt_band_get_data(rtband);
1965 RASTER_DEBUGF(4, "Band data is at pos %p", pVoid);
1966
1967 pszDataPointer = (char *) rtalloc(20 * sizeof (char));
1968 sprintf(pszDataPointer, "%p", pVoid);
1969 RASTER_DEBUGF(4, "rt_raster_to_gdal_mem: szDatapointer is %p",
1970 pszDataPointer);
1971
1972 if (strncasecmp(pszDataPointer, "0x", 2) == 0)
1973 sprintf(szGDALOption, "DATAPOINTER=%s", pszDataPointer);
1974 else
1975 sprintf(szGDALOption, "DATAPOINTER=0x%s", pszDataPointer);
1976
1977 RASTER_DEBUG(3, "Storing info for GDAL MEM raster band");
1978
1979 apszOptions[0] = szGDALOption;
1980 apszOptions[1] = NULL; /* pixel offset, not needed */
1981 apszOptions[2] = NULL; /* line offset, not needed */
1982 apszOptions[3] = NULL;
1983
1984 /* free */
1985 rtdealloc(pszDataPointer);
1986
1987 /* add band */
1988 if (GDALAddBand(ds, gdal_pt, apszOptions) == CE_Failure) {
1989 rterror("rt_raster_to_gdal_mem: Could not add GDAL raster band");
1990 if (allocBandNums) rtdealloc(bandNums);
1991 GDALClose(ds);
1992 return 0;
1993 }
1994#if POSTGIS_GDAL_VERSION < 30700
1995 /*
1996 PT_8BSI is special as GDAL (prior to 3.7) has no equivalent pixel type.
1997 Must convert 8BSI to 16BSI so create basic band
1998 */
1999 } else {
2000 /* add band */
2001 if (GDALAddBand(ds, gdal_pt, NULL) == CE_Failure) {
2002 rterror("rt_raster_to_gdal_mem: Could not add GDAL raster band");
2003 if (allocBandNums) rtdealloc(bandNums);
2004 if (allocNodataValues) rtdealloc(excludeNodataValues);
2005 GDALClose(ds);
2006 return 0;
2007 }
2008 }
2009#endif
2010
2011 /* check band count */
2012 if (GDALGetRasterCount(ds) != i + 1) {
2013 rterror("rt_raster_to_gdal_mem: Error creating GDAL MEM raster band");
2014 if (allocBandNums) rtdealloc(bandNums);
2015 if (allocNodataValues) rtdealloc(excludeNodataValues);
2016 GDALClose(ds);
2017 return 0;
2018 }
2019
2020 /* get new band */
2021 band = NULL;
2022 band = GDALGetRasterBand(ds, i + 1);
2023 if (NULL == band) {
2024 rterror("rt_raster_to_gdal_mem: Could not get GDAL band for additional processing");
2025 if (allocBandNums) rtdealloc(bandNums);
2026 if (allocNodataValues) rtdealloc(excludeNodataValues);
2027 GDALClose(ds);
2028 return 0;
2029 }
2030
2031#if POSTGIS_GDAL_VERSION < 30700
2032
2033 /* PT_8BSI requires manual setting of pixels */
2034 if (pt == PT_8BSI) {
2035 uint32_t nXBlocks, nYBlocks;
2036 int nXBlockSize, nYBlockSize;
2037 uint32_t iXBlock, iYBlock;
2038 int nXValid, nYValid;
2039 int iX, iY;
2040 int iXMax, iYMax;
2041
2042 int x, y, z;
2043 uint32_t valueslen = 0;
2044 int16_t *values = NULL;
2045 double value = 0.;
2046
2047 /* this makes use of GDAL's "natural" blocks */
2048 GDALGetBlockSize(band, &nXBlockSize, &nYBlockSize);
2049 nXBlocks = (width + nXBlockSize - 1) / nXBlockSize;
2050 nYBlocks = (height + nYBlockSize - 1) / nYBlockSize;
2051 RASTER_DEBUGF(4, "(nXBlockSize, nYBlockSize) = (%d, %d)", nXBlockSize, nYBlockSize);
2052 RASTER_DEBUGF(4, "(nXBlocks, nYBlocks) = (%d, %d)", nXBlocks, nYBlocks);
2053
2054 /* length is for the desired pixel type */
2055 valueslen = rt_pixtype_size(PT_16BSI) * nXBlockSize * nYBlockSize;
2056 values = rtalloc(valueslen);
2057 if (NULL == values) {
2058 rterror("rt_raster_to_gdal_mem: Could not allocate memory for GDAL band pixel values");
2059 if (allocBandNums) rtdealloc(bandNums);
2060 if (allocNodataValues) rtdealloc(excludeNodataValues);
2061 GDALClose(ds);
2062 return 0;
2063 }
2064
2065 for (iYBlock = 0; iYBlock < nYBlocks; iYBlock++) {
2066 for (iXBlock = 0; iXBlock < nXBlocks; iXBlock++) {
2067 memset(values, 0, valueslen);
2068
2069 x = iXBlock * nXBlockSize;
2070 y = iYBlock * nYBlockSize;
2071 RASTER_DEBUGF(4, "(iXBlock, iYBlock) = (%d, %d)", iXBlock, iYBlock);
2072 RASTER_DEBUGF(4, "(x, y) = (%d, %d)", x, y);
2073
2074 /* valid block width */
2075 if ((iXBlock + 1) * nXBlockSize > width)
2076 nXValid = width - (iXBlock * nXBlockSize);
2077 else
2078 nXValid = nXBlockSize;
2079
2080 /* valid block height */
2081 if ((iYBlock + 1) * nYBlockSize > height)
2082 nYValid = height - (iYBlock * nYBlockSize);
2083 else
2084 nYValid = nYBlockSize;
2085
2086 RASTER_DEBUGF(4, "(nXValid, nYValid) = (%d, %d)", nXValid, nYValid);
2087
2088 /* convert 8BSI values to 16BSI */
2089 z = 0;
2090 iYMax = y + nYValid;
2091 iXMax = x + nXValid;
2092 for (iY = y; iY < iYMax; iY++) {
2093 for (iX = x; iX < iXMax; iX++) {
2094 if (rt_band_get_pixel(rtband, iX, iY, &value, NULL) != ES_NONE) {
2095 rterror("rt_raster_to_gdal_mem: Could not get pixel value to convert from 8BSI to 16BSI");
2096 rtdealloc(values);
2097 if (allocBandNums) rtdealloc(bandNums);
2098 if (allocNodataValues) rtdealloc(excludeNodataValues);
2099 GDALClose(ds);
2100 return 0;
2101 }
2102
2103 values[z++] = rt_util_clamp_to_16BSI(value);
2104 }
2105 }
2106
2107 /* burn values */
2108 if (GDALRasterIO(
2109 band, GF_Write,
2110 x, y,
2111 nXValid, nYValid,
2112 values, nXValid, nYValid,
2113 gdal_pt,
2114 0, 0
2115 ) != CE_None) {
2116 rterror("rt_raster_to_gdal_mem: Could not write converted 8BSI to 16BSI values to GDAL band");
2117 rtdealloc(values);
2118 if (allocBandNums) rtdealloc(bandNums);
2119 if (allocNodataValues) rtdealloc(excludeNodataValues);
2120 GDALClose(ds);
2121 return 0;
2122 }
2123 }
2124 }
2125
2126 rtdealloc(values);
2127 }
2128#endif
2129
2130 /* Add nodata value for band */
2131 if (rt_band_get_hasnodata_flag(rtband) != FALSE && excludeNodataValues[i]) {
2132 rt_band_get_nodata(rtband, &nodata);
2133 if (GDALSetRasterNoDataValue(band, nodata) != CE_None)
2134 rtwarn("rt_raster_to_gdal_mem: Could not set nodata value for band");
2135 RASTER_DEBUGF(3, "nodata value set to %f", GDALGetRasterNoDataValue(band, NULL));
2136 }
2137
2138#if POSTGIS_DEBUG_LEVEL > 3
2139 {
2140 GDALRasterBandH _grb = NULL;
2141 double _min;
2142 double _max;
2143 double _mean;
2144 double _stddev;
2145
2146 _grb = GDALGetRasterBand(ds, i + 1);
2147 GDALComputeRasterStatistics(_grb, FALSE, &_min, &_max, &_mean, &_stddev, NULL, NULL);
2148 RASTER_DEBUGF(4, "GDAL Band %d stats: %f, %f, %f, %f", i + 1, _min, _max, _mean, _stddev);
2149 }
2150#endif
2151
2152 }
2153
2154 /* necessary??? */
2155 GDALFlushCache(ds);
2156
2157 if (allocBandNums) rtdealloc(bandNums);
2158 if (allocNodataValues) rtdealloc(excludeNodataValues);
2159
2160 return ds;
2161}
FALSE
#define FALSE
Definition dbfopen.c:72
rterror
void rterror(const char *fmt,...) __attribute__((format(printf
Wrappers used for reporting errors and info.
rtalloc
void * rtalloc(size_t size)
Wrappers used for managing memory.
Definition rt_context.c:191
RASTER_DEBUG
#define RASTER_DEBUG(level, msg)
Definition librtcore.h:304
RASTER_DEBUGF
#define RASTER_DEBUGF(level, msg,...)
Definition librtcore.h:308
rtwarn
void void void rtwarn(const char *fmt,...) __attribute__((format(printf
rt_band_get_hasnodata_flag
int rt_band_get_hasnodata_flag(rt_band band)
Get hasnodata flag value.
Definition rt_band.c:833
rt_band_get_data
void * rt_band_get_data(rt_band band)
Get pointer to raster band data.
Definition rt_band.c:559
rt_band_get_pixel
rt_errorstate rt_band_get_pixel(rt_band band, int x, int y, double *value, int *nodata)
Get pixel value.
Definition rt_band.c:1551
rt_pixtype
rt_pixtype
Definition librtcore.h:188
PT_END
@ PT_END
Definition librtcore.h:201
PT_8BSI
@ PT_8BSI
Definition librtcore.h:192
PT_16BSI
@ PT_16BSI
Definition librtcore.h:194
rt_util_gdal_convert_sr
char * rt_util_gdal_convert_sr(const char *srs, int proj4)
Definition rt_util.c:322
rt_util_pixtype_to_gdal_datatype
GDALDataType rt_util_pixtype_to_gdal_datatype(rt_pixtype pt)
Convert rt_pixtype to GDALDataType.
Definition rt_util.c:213
ES_NONE
@ ES_NONE
Definition librtcore.h:182
rt_util_clamp_to_16BSI
int16_t rt_util_clamp_to_16BSI(double value)
Definition rt_util.c:61
rt_util_gdal_driver_registered
int rt_util_gdal_driver_registered(const char *drv)
Definition rt_util.c:463
rt_band_get_nodata
rt_errorstate rt_band_get_nodata(rt_band band, double *nodata)
Get NODATA value.
Definition rt_band.c:2067
rt_band_get_pixtype
rt_pixtype rt_band_get_pixtype(rt_band band)
Return pixeltype of this band.
Definition rt_band.c:790
rtdealloc
void rtdealloc(void *mem)
Definition rt_context.c:206
rt_pixtype_size
int rt_pixtype_size(rt_pixtype pixtype)
Return size in bytes of a value in the given pixtype.
Definition rt_pixel.c:40
genraster.value
int value
Definition genraster.py:62
genraster.count
int count
Definition genraster.py:57
rt_raster_get_num_bands
uint16_t rt_raster_get_num_bands(rt_raster raster)
Definition rt_raster.c:376
rt_raster_get_height
uint16_t rt_raster_get_height(rt_raster raster)
Definition rt_raster.c:133
rt_raster_get_band
rt_band rt_raster_get_band(rt_raster raster, int n)
Return Nth band, or NULL if unavailable.
Definition rt_raster.c:385
rt_raster_get_width
uint16_t rt_raster_get_width(rt_raster raster)
Definition rt_raster.c:125
rt_raster_get_geotransform_matrix
void rt_raster_get_geotransform_matrix(rt_raster raster, double *gt)
Get 6-element array of raster geotransform matrix.
Definition rt_raster.c:588
rt_band_t
Definition librtcore.h:2515

References ES_NONE, FALSE, PT_16BSI, PT_8BSI, PT_END, RASTER_DEBUG, RASTER_DEBUGF, rt_band_get_data(), rt_band_get_hasnodata_flag(), rt_band_get_nodata(), rt_band_get_pixel(), rt_band_get_pixtype(), rt_pixtype_size(), rt_raster_get_band(), rt_raster_get_geotransform_matrix(), rt_raster_get_height(), rt_raster_get_num_bands(), rt_raster_get_width(), rt_util_clamp_to_16BSI(), rt_util_gdal_convert_sr(), rt_util_gdal_driver_registered(), rt_util_pixtype_to_gdal_datatype(), rtalloc(), rtdealloc(), rterror(), and rtwarn().

Referenced by rt_raster_gdal_contour(), rt_raster_gdal_polygonize(), rt_raster_gdal_warp(), rt_raster_to_gdal(), and test_gdal_to_raster().

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