#include <math.h>
#include <stdio.h> /* for printf (default message handler) */
#include <stdarg.h> /* for va_list, va_start etc */
-#include <string.h> /* for memcpy */
+#include <string.h> /* for memcpy and strlen */
#include <assert.h>
#include <float.h> /* for FLT_EPSILON and float type limits */
#include <limits.h> /* for integer type limits */
return 1;
}
+/**
+ * Return formatted GDAL raster from raster
+ *
+ * @param raster : the raster to convert
+ * @param srs : the raster's coordinate system in OGC WKT or PROJ.4
+ * @param format : format to convert to. GDAL driver short name
+ * @param options : list of format creation options. array of strings
+ * @param gdalsize : will be set to the size of returned bytea
+ *
+ * @return formatted GDAL raster
+ */
+uint8_t*
+rt_raster_to_gdal(rt_raster raster, char *srs,
+ char *format, char **options, uint64_t *gdalsize) {
+ GDALDriverH src_drv = NULL;
+ GDALDatasetH src_ds = NULL;
+
+ vsi_l_offset rtn_lenvsi;
+ uint64_t rtn_len = 0;
+
+ GDALDriverH rtn_drv = NULL;
+ GDALDatasetH rtn_ds = NULL;
+ uint8_t *rtn = NULL;
+
+ assert(NULL != raster);
+
+ /* any supported format is possible */
+ GDALAllRegister();
+ RASTER_DEBUG(3, "rt_raster_to_gdal: loaded all supported GDAL formats");
+
+ /* output format not specified */
+ if (format == NULL || !strlen(format))
+ format = "GTiff";
+ RASTER_DEBUGF(3, "rt_raster_to_gdal: output format is %s", format);
+
+ /* load raster into a GDAL MEM raster */
+ src_ds = rt_raster_to_gdal_mem(raster, srs, &src_drv);
+ if (NULL == src_ds) {
+ rterror("rt_raster_to_gdal: Unable to convert raster to GDAL MEM format\n");
+ if (NULL != src_drv) {
+ GDALDeregisterDriver(src_drv);
+ GDALDestroyDriver(src_drv);
+ }
+ return 0;
+ }
+
+ /* load driver */
+ rtn_drv = GDALGetDriverByName(format);
+ if (NULL == rtn_drv) {
+ rterror("rt_raster_to_gdal: Unable to load the output GDAL driver\n");
+ GDALClose(src_ds);
+ GDALDeregisterDriver(src_drv);
+ GDALDestroyDriver(src_drv);
+ return 0;
+ }
+ RASTER_DEBUG(3, "rt_raster_to_gdal: Output driver loaded");
+
+ /* convert GDAL MEM raster to output format */
+ RASTER_DEBUG(3, "rt_raster_to_gdal: Copying GDAL MEM raster to memory file in output format");
+ rtn_ds = GDALCreateCopy(
+ rtn_drv,
+ "/vsimem/out.dat", /* should be fine assuming this is in a process */
+ src_ds,
+ FALSE, /* should copy be strictly equivelent? */
+ options, /* format options */
+ NULL, /* progress function */
+ NULL /* progress data */
+ );
+ if (NULL == rtn_ds) {
+ rterror("rt_raster_to_gdal: Unable to create the output GDAL dataset\n");
+ GDALClose(src_ds);
+ GDALDeregisterDriver(rtn_drv);
+ GDALDestroyDriver(rtn_drv);
+ GDALDeregisterDriver(src_drv);
+ GDALDestroyDriver(src_drv);
+ return 0;
+ }
+
+ /* close source dataset */
+ GDALClose(src_ds);
+ GDALDeregisterDriver(src_drv);
+ GDALDestroyDriver(src_drv);
+ RASTER_DEBUG(3, "rt_raster_to_gdal: Closed GDAL MEM raster");
+
+ /* close dataset, this also flushes any pending writes */
+ GDALClose(rtn_ds);
+ GDALDeregisterDriver(rtn_drv);
+ GDALDestroyDriver(rtn_drv);
+ RASTER_DEBUG(3, "rt_raster_to_gdal: Closed GDAL output raster");
+
+ RASTER_DEBUG(3, "rt_raster_to_gdal: Done copying GDAL MEM raster to memory file in output format");
+
+ /* from memory file to buffer */
+ RASTER_DEBUG(3, "rt_raster_to_gdal: Copying GDAL memory file to buffer");
+ rtn = VSIGetMemFileBuffer("/vsimem/out.dat", &rtn_lenvsi, TRUE);
+ RASTER_DEBUG(3, "rt_raster_to_gdal: Done copying GDAL memory file to buffer");
+ if (NULL == rtn) {
+ rterror("rt_raster_to_gdal: Unable to create the output GDAL raster\n");
+ return 0;
+ }
+
+ rtn_len = (uint64_t) rtn_lenvsi;
+ *gdalsize = rtn_len;
+
+ return rtn;
+}
+
+struct rt_gdaldriver_t {
+ int idx;
+ char *short_name;
+ char *long_name;
+ char *create_options;
+};
+
+/**
+ * Returns a set of available GDAL drivers
+ *
+ * @param drv_count : number of GDAL drivers available
+ *
+ * @return set of "gdaldriver" values of available GDAL drivers
+ */
+rt_gdaldriver
+rt_raster_gdal_drivers(uint32_t *drv_count) {
+ const char *state;
+ const char *txt;
+ int txt_len;
+ GDALDriverH *drv = NULL;
+ rt_gdaldriver rtn = NULL;
+ int count;
+ int i;
+ uint32_t j;
+
+ GDALAllRegister();
+ count = GDALGetDriverCount();
+ rtn = (rt_gdaldriver) rtalloc(count * sizeof(struct rt_gdaldriver_t));
+ if (NULL == rtn) {
+ rterror("rt_raster_gdal_drivers: Unable to allocate memory for gdaldriver structure");
+ return 0;
+ }
+
+ for (i = 0, j = 0; i < count; i++) {
+ drv = GDALGetDriver(i);
+
+ /* CreateCopy support */
+ state = GDALGetMetadataItem(drv, GDAL_DCAP_CREATECOPY, NULL);
+ if (state == NULL) continue;
+
+ /* VirtualIO support */
+ state = GDALGetMetadataItem(drv, GDAL_DCAP_VIRTUALIO, NULL);
+ if (state == NULL) continue;
+
+ /* index of driver */
+ rtn[j].idx = i;
+
+ /* short name */
+ txt = GDALGetDriverShortName(drv);
+ txt_len = strlen(txt);
+
+ RASTER_DEBUGF(3, "rt_raster_gdal_driver: driver %s (%d) supports CreateCopy() and VirtualIO()", txt, i);
+
+ txt_len = (txt_len + 1) * sizeof(char);
+ rtn[j].short_name = (char *) rtalloc(txt_len);
+ memcpy(rtn[j].short_name, txt, txt_len);
+
+ /* long name */
+ txt = GDALGetDriverLongName(drv);
+ txt_len = strlen(txt);
+
+ txt_len = (txt_len + 1) * sizeof(char);
+ rtn[j].long_name = (char *) rtalloc(txt_len);
+ memcpy(rtn[j].long_name, txt, txt_len);
+
+ /* creation options */
+ txt = GDALGetDriverCreationOptionList(drv);
+ txt_len = strlen(txt);
+
+ txt_len = (txt_len + 1) * sizeof(char);
+ rtn[j].create_options = (char *) rtalloc(txt_len);
+ memcpy(rtn[j].create_options, txt, txt_len);
+
+ j++;
+ }
+
+ /* free unused memory */
+ rtrealloc(rtn, j * sizeof(struct rt_gdaldriver_t));
+ *drv_count = j;
+
+ return rtn;
+}
+
+/**
+ * Return GDAL datasource using GDAL MEM driver from raster
+ *
+ * @param raster : raster to convert to GDAL MEM
+ * @param srs : the raster's coordinate system in OGC WKT or PROJ.4
+ * @param rtn_drv : is set to the GDAL driver object
+ *
+ * @return GDAL datasource using GDAL MEM driver
+ */
+GDALDatasetH
+rt_raster_to_gdal_mem(rt_raster raster, char *srs,
+ GDALDriverH *rtn_drv) {
+ GDALDriverH drv = NULL;
+ int drv_gen = 0;
+ GDALDatasetH ds = NULL;
+ double gt[] = {0.0, 1.0, 0.0, 0.0, 0.0, 1.0};
+ CPLErr cplerr;
+ GDALDataType gdal_pt = GDT_Unknown;
+ GDALRasterBandH band;
+ void *pVoid;
+ char *pszDataPointer;
+ char szGDALOption[50];
+ char *apszOptions[4];
+ double nodata = 0.0;
+
+ int i;
+ int numBands;
+ rt_band rtband = NULL;
+ rt_pixtype pt = PT_END;
+
+ assert(NULL != raster);
+
+ /* store raster in GDAL MEM raster */
+ GDALRegister_MEM();
+ if (NULL == *rtn_drv) {
+ drv = GDALGetDriverByName("MEM");
+ *rtn_drv = drv;
+ drv_gen = 1;
+ if (NULL == drv) {
+ rterror("rt_raster_to_gdal_mem: Unable to load the MEM GDAL driver\n");
+ return 0;
+ }
+ }
+
+ ds = GDALCreate(drv, "", rt_raster_get_width(raster),
+ rt_raster_get_height(raster), 0, GDT_Byte, NULL);
+ if (NULL == ds) {
+ rterror("rt_raster_to_gdal_mem: Couldn't create a GDALDataset to convert it\n");
+ if (drv_gen) {
+ GDALDeregisterDriver(drv);
+ GDALDestroyDriver(drv);
+ }
+ return 0;
+ }
+
+ /* add geotransform */
+ gt[0] = rt_raster_get_x_offset(raster);
+ gt[1] = rt_raster_get_x_scale(raster);
+ gt[2] = rt_raster_get_x_skew(raster);
+ gt[3] = rt_raster_get_y_offset(raster);
+ gt[4] = rt_raster_get_y_skew(raster);
+ gt[5] = rt_raster_get_y_scale(raster);
+ cplerr = GDALSetGeoTransform(ds, gt);
+ if (cplerr != CE_None) {
+ rterror("rt_raster_to_gdal_mem: Unable to set geotransformation\n");
+ GDALClose(ds);
+ if (drv_gen) {
+ GDALDeregisterDriver(drv);
+ GDALDestroyDriver(drv);
+ }
+ return 0;
+ }
+
+ /* set spatial reference */
+ if (NULL != srs && strlen(srs)) {
+ cplerr = GDALSetProjection(ds, srs);
+ if (cplerr != CE_None) {
+ rterror("rt_raster_to_gdal_mem: Unable to set projection\n");
+ GDALClose(ds);
+ if (drv_gen) {
+ GDALDeregisterDriver(drv);
+ GDALDestroyDriver(drv);
+ }
+ return 0;
+ }
+ }
+
+ /* add band(s) */
+ numBands = rt_raster_get_num_bands(raster);
+ for (i = 0; i < numBands; i++) {
+ rtband = rt_raster_get_band(raster, i);
+ if (0 == rtband) {
+ rterror("rt_raster_to_gdal_mem: Unable to get requested band\n");
+ GDALClose(ds);
+ if (drv_gen) {
+ GDALDeregisterDriver(drv);
+ GDALDestroyDriver(drv);
+ }
+ return 0;
+ }
+
+ pt = rt_band_get_pixtype(rtband);
+ switch (pt) {
+ case PT_1BB: case PT_2BUI: case PT_4BUI: case PT_8BSI: case PT_8BUI:
+ gdal_pt = GDT_Byte;
+ break;
+ case PT_16BSI: case PT_16BUI:
+ if (pt == PT_16BSI)
+ gdal_pt = GDT_Int16;
+ else
+ gdal_pt = GDT_UInt16;
+ break;
+ case PT_32BSI: case PT_32BUI: case PT_32BF:
+ if (pt == PT_32BSI)
+ gdal_pt = GDT_Int32;
+ else if (pt == PT_32BUI)
+ gdal_pt = GDT_UInt32;
+ else
+ gdal_pt = GDT_Float32;
+ break;
+ case PT_64BF:
+ gdal_pt = GDT_Float64;
+ break;
+ default:
+ rtwarn("rt_raster_to_gdal_mem: Unknown pixel type for band\n");
+ gdal_pt = GDT_Unknown;
+ break;
+ }
+
+ pVoid = rt_band_get_data(rtband);
+ RASTER_DEBUGF(4, "Band data is at pos %p", pVoid);
+
+ pszDataPointer = (char *) rtalloc(20 * sizeof (char));
+ sprintf(pszDataPointer, "%p", pVoid);
+ RASTER_DEBUGF(4, "rt_raster_to_gdal_mem: szDatapointer is %p",
+ pszDataPointer);
+
+ if (strnicmp(pszDataPointer, "0x", 2) == 0)
+ sprintf(szGDALOption, "DATAPOINTER=%s", pszDataPointer);
+ else
+ sprintf(szGDALOption, "DATAPOINTER=0x%s", pszDataPointer);
+
+ RASTER_DEBUG(3, "Storing info for GDAL MEM raster band");
+
+ apszOptions[0] = szGDALOption;
+ apszOptions[1] = NULL; /* pixel offset, not needed */
+ apszOptions[2] = NULL; /* line offset, not needed */
+ apszOptions[3] = NULL;
+
+ /* free */
+ rtdealloc(pszDataPointer);
+
+ if (GDALAddBand(ds, gdal_pt, apszOptions) == CE_Failure) {
+ rterror("rt_raster_to_gdal_mem: Couldn't add GDALRasterBand\n");
+ GDALClose(ds);
+ if (drv_gen) {
+ GDALDeregisterDriver(drv);
+ GDALDestroyDriver(drv);
+ }
+ return 0;
+ }
+
+ /* check band count */
+ if (GDALGetRasterCount(ds) != i + 1) {
+ rterror("rt_raster_to_gdal_mem: Error creating GDAL MEM raster band\n");
+ GDALClose(ds);
+ if (drv_gen) {
+ GDALDeregisterDriver(drv);
+ GDALDestroyDriver(drv);
+ }
+ return 0;
+ }
+
+ /* get band to write data to */
+ band = NULL;
+ band = GDALGetRasterBand(ds, i + 1);
+ if (NULL == band) {
+ rterror("rt_raster_to_gdal_mem: Couldn't get GDAL band for additional processing\n");
+ GDALClose(ds);
+ if (drv_gen) {
+ GDALDeregisterDriver(drv);
+ GDALDestroyDriver(drv);
+ }
+ return 0;
+ }
+
+ /* Add nodata value for band */
+ if (rt_band_get_hasnodata_flag(rtband) != FALSE) {
+ nodata = rt_band_get_nodata(rtband);
+ if (GDALSetRasterNoDataValue(band, nodata) != CE_None)
+ rtwarn("rt_raster_to_gdal_mem: Couldn't set nodata value for band\n");
+ }
+ }
+
+ return ds;
+}
#include "gdal_frmts.h"
#include "gdal.h"
#include "ogr_api.h"
+#include "ogr_srs_api.h"
+#include "cpl_vsi.h"
#include "../../postgis_config.h"
/**
int rt_raster_replace_band(rt_raster raster, rt_band band,
int index);
+/**
+ * Return formatted GDAL raster from raster
+ *
+ * @param raster : the raster to convert
+ * @param srs : the raster's coordinate system in OGC WKT or PROJ.4
+ * @param format : format to convert to. GDAL driver short name
+ * @param options : list of format creation options. array of strings
+ * @param gdalsize : will be set to the size of returned bytea
+ *
+ * @return formatted GDAL raster
+ */
+uint8_t *rt_raster_to_gdal(rt_raster raster, char *srs,
+ char *format, char **options, uint64_t *gdalsize);
+
+/**
+ * Returns a set of available GDAL drivers
+ *
+ * @param drv_count : number of GDAL drivers available
+ *
+ * @return set of "gdaldriver" values of available GDAL drivers
+ */
+typedef struct rt_gdaldriver_t* rt_gdaldriver;
+rt_gdaldriver rt_raster_gdal_drivers(uint32_t *drv_count);
+
+/**
+ * Return GDAL datasource using GDAL MEM driver from raster
+ *
+ * @param raster : raster to convert to GDAL MEM
+ * @param srs : the raster's coordinate system in OGC WKT or PROJ.4
+ * @param rtn_drv : is set to the GDAL driver object
+ *
+ * @return GDAL datasource using GDAL MEM driver
+ */
+GDALDatasetH rt_raster_to_gdal_mem(rt_raster raster, char *srs,
+ GDALDriverH *rtn_drv);
/*- utilities -------------------------------------------------------*/
static char *chartrim(char* input, char *remove); /* for RASTER_reclass */
static char **strsplit(const char *str, const char *delimiter, int *n); /* for RASTER_reclass */
static char *removespaces(char *str); /* for RASTER_reclass */
+static char *trim(char* input); /* for RASTER_asGDALRaster */
/***************************************************************
* Some rules for returning NOTICE or ERROR...
/* reclassify specified bands of a raster */
Datum RASTER_reclass(PG_FUNCTION_ARGS);
+/* convert raster to GDAL raster */
+Datum RASTER_asGDALRaster(PG_FUNCTION_ARGS);
+Datum RASTER_getGDALDrivers(PG_FUNCTION_ARGS);
/* Replace function taken from
* http://ubuntuforums.org/showthread.php?s=aa6f015109fd7e4c7e30d2fd8b717497&t=141670&page=3
return rtn;
}
+static char*
+trim(char *input) {
+ char *start;
+ char *ptr;
+
+ if (!input) {
+ return NULL;
+ }
+ else if (!*input) {
+ return input;
+ }
+
+ start = (char *) palloc(sizeof(char) * strlen(input) + 1);
+ if (NULL == start) {
+ fprintf(stderr, "Not enough memory\n");
+ return NULL;
+ }
+ strcpy(start, input);
+
+ /* trim left */
+ while (isspace(*start)) {
+ start++;
+ }
+
+ /* trim right */
+ ptr = start + strlen(start);
+ while (isspace(*--ptr));
+ *(++ptr) = '\0';
+
+ return start;
+}
+
PG_FUNCTION_INFO_V1(RASTER_lib_version);
Datum RASTER_lib_version(PG_FUNCTION_ARGS)
{
PG_RETURN_POINTER(pgrast);
}
+/**
+ * Returns formatted GDAL raster as bytea object of raster
+ */
+PG_FUNCTION_INFO_V1(RASTER_asGDALRaster);
+Datum RASTER_asGDALRaster(PG_FUNCTION_ARGS)
+{
+ rt_pgraster *pgraster = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
+ rt_raster raster;
+
+ text *formattext = NULL;
+ char *format = NULL;
+ char **options = NULL;
+ text *optiontext = NULL;
+ char *option = NULL;
+ text *srstext = NULL;
+ char *srs = NULL;
+
+ ArrayType *array;
+ Oid etype;
+ Datum *e;
+ bool *nulls;
+ int16 typlen;
+ bool typbyval;
+ char typalign;
+ int ndims = 1;
+ int *dims;
+ int *lbs;
+ int n = 0;
+ int i = 0;
+ int j = 0;
+
+ uint8_t *gdal = NULL;
+ uint64_t gdal_size = 0;
+ bytea *result = NULL;
+ uint64_t result_size = 0;
+
+ POSTGIS_RT_DEBUG(3, "RASTER_asGDALRaster: Starting");
+
+ raster = rt_raster_deserialize(pgraster);
+ if (!raster) {
+ elog(ERROR, "RASTER_asGDALRaster: Could not deserialize raster");
+ PG_RETURN_NULL();
+ }
+
+ /* format is required */
+ if (PG_ARGISNULL(1)) {
+ elog(NOTICE, "Format must be provided");
+ rt_raster_destroy(raster);
+ PG_RETURN_NULL();
+ }
+ else {
+ formattext = PG_GETARG_TEXT_P(1);
+ format = text_to_cstring(formattext);
+ }
+
+ POSTGIS_RT_DEBUGF(3, "RASTER_asGDALRaster: Arg 1 (format) is %s", format);
+
+ /* process options */
+ if (!PG_ARGISNULL(2)) {
+ POSTGIS_RT_DEBUG(3, "RASTER_asGDALRaster: Processing Arg 2 (options)");
+ array = PG_GETARG_ARRAYTYPE_P(2);
+ etype = ARR_ELEMTYPE(array);
+ get_typlenbyvalalign(etype, &typlen, &typbyval, &typalign);
+
+ switch (etype) {
+ case TEXTOID:
+ break;
+ default:
+ elog(ERROR, "RASTER_asGDALRaster: Invalid data type for options");
+ PG_RETURN_NULL();
+ break;
+ }
+
+ ndims = ARR_NDIM(array);
+ dims = ARR_DIMS(array);
+ lbs = ARR_LBOUND(array);
+
+ deconstruct_array(array, etype, typlen, typbyval, typalign, &e,
+ &nulls, &n);
+
+ if (n) {
+ options = (char **) palloc(sizeof(char *) * (n + 1));
+ /* clean each option */
+ for (i = 0, j = 0; i < n; i++) {
+ if (nulls[i]) continue;
+
+ option = NULL;
+ switch (etype) {
+ case TEXTOID:
+ optiontext = (text *) DatumGetPointer(e[i]);
+ if (NULL == optiontext) break;
+ option = text_to_cstring(optiontext);
+
+ /* trim string */
+ option = trim(option);
+ POSTGIS_RT_DEBUGF(3, "RASTER_asGDALRaster: option is '%s'", option);
+ break;
+ }
+
+ if (strlen(option)) {
+ options[j] = (char *) palloc(sizeof(char) * (strlen(option) + 1));
+ options[j] = option;
+ j++;
+ }
+ }
+
+ if (j > 0) {
+ /* add NULL to end */
+ options[j] = (char *) palloc(sizeof(char));
+ options[j] = '\0';
+
+ /* trim allocation */
+ options = repalloc(options, j * sizeof(char *));
+ }
+ else {
+ pfree(options);
+ options = NULL;
+ }
+ }
+ }
+
+ /* process srs */
+ if (!PG_ARGISNULL(3)) {
+ srstext = PG_GETARG_TEXT_P(3);
+ srs = text_to_cstring(srstext);
+ POSTGIS_RT_DEBUGF(3, "RASTER_asGDALRaster: Arg 3 (srs) is %s", srs);
+ }
+
+ POSTGIS_RT_DEBUG(3, "RASTER_asGDALRaster: Generating GDAL raster");
+ gdal = rt_raster_to_gdal(raster, srs, format, options, &gdal_size);
+
+ /* free memory */
+ if (NULL != options) {
+ for (i = j - 1; i >= 0; i--)
+ pfree(options[i]);
+ pfree(options);
+ }
+ rt_raster_destroy(raster);
+
+ if (!gdal) {
+ elog(ERROR, "RASTER_asGDALRaster: Could not allocate and generate GDAL raster");
+ PG_RETURN_NULL();
+ }
+ POSTGIS_RT_DEBUGF(3, "RASTER_asGDALRaster: GDAL raster generated with %d bytes", (int) gdal_size);
+
+ result_size = gdal_size + VARHDRSZ;
+ result = (bytea *) palloc(result_size);
+ if (NULL == result) {
+ elog(ERROR, "RASTER_asGDALRaster: Insufficient virtual memory for GDAL raster");
+ PG_RETURN_NULL();
+ }
+ SET_VARSIZE(result, result_size);
+ memcpy(VARDATA(result), gdal, VARSIZE(result) - VARHDRSZ);
+
+ /* for test output
+ FILE *fh = NULL;
+ fh = fopen("/tmp/out.dat", "w");
+ fwrite(gdal, sizeof(uint8_t), gdal_size, fh);
+ fclose(fh);
+ */
+
+ POSTGIS_RT_DEBUG(3, "RASTER_asGDALRaster: Returning pointer to GDAL raster");
+ PG_RETURN_POINTER(result);
+}
+
+/**
+ * Needed for sizeof
+ */
+struct rt_gdaldriver_t {
+ int idx;
+ char *short_name;
+ char *long_name;
+ char *create_options;
+};
+
+/**
+ * Returns available GDAL drivers
+ */
+PG_FUNCTION_INFO_V1(RASTER_getGDALDrivers);
+Datum RASTER_getGDALDrivers(PG_FUNCTION_ARGS)
+{
+ FuncCallContext *funcctx;
+ TupleDesc tupdesc;
+ AttInMetadata *attinmeta;
+
+ uint32_t drv_count;
+ rt_gdaldriver drv_set;
+ rt_gdaldriver drv_set2;
+ int call_cntr;
+ int max_calls;
+
+ /* first call of function */
+ if (SRF_IS_FIRSTCALL()) {
+ MemoryContext oldcontext;
+
+ /* create a function context for cross-call persistence */
+ funcctx = SRF_FIRSTCALL_INIT();
+
+ /* switch to memory context appropriate for multiple function calls */
+ oldcontext = MemoryContextSwitchTo(funcctx->multi_call_memory_ctx);
+
+ drv_set = rt_raster_gdal_drivers(&drv_count);
+ if (NULL == drv_set || !drv_count) {
+ elog(NOTICE, "No GDAL drivers found");
+ SRF_RETURN_DONE(funcctx);
+ }
+
+ POSTGIS_RT_DEBUGF(3, "%d drivers returned", (int) drv_count);
+
+ /* Store needed information */
+ funcctx->user_fctx = drv_set;
+
+ /* total number of tuples to be returned */
+ funcctx->max_calls = drv_count;
+
+ /* Build a tuple descriptor for our result type */
+ if (get_call_result_type(fcinfo, NULL, &tupdesc) != TYPEFUNC_COMPOSITE) {
+ ereport(ERROR, (
+ errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
+ errmsg(
+ "function returning record called in context "
+ "that cannot accept type record"
+ )
+ ));
+ }
+
+ /*
+ * generate attribute metadata needed later to produce tuples from raw
+ * C strings
+ */
+ attinmeta = TupleDescGetAttInMetadata(tupdesc);
+ funcctx->attinmeta = attinmeta;
+ MemoryContextSwitchTo(oldcontext);
+ }
+
+ /* stuff done on every call of the function */
+ funcctx = SRF_PERCALL_SETUP();
+
+ call_cntr = funcctx->call_cntr;
+ max_calls = funcctx->max_calls;
+ attinmeta = funcctx->attinmeta;
+ drv_set2 = funcctx->user_fctx;
+
+ /* do when there is more left to send */
+ if (call_cntr < max_calls) {
+ char **values;
+ HeapTuple tuple;
+ Datum result;
+
+ POSTGIS_RT_DEBUGF(3, "Result %d", call_cntr);
+
+ /*
+ * Prepare a values array for building the returned tuple.
+ * This should be an array of C strings which will
+ * be processed later by the type input functions.
+ */
+ values = (char **) palloc(4 * sizeof(char *));
+
+ values[0] = (char *) palloc(
+ (snprintf(NULL, 0, "%d", drv_set2[call_cntr].idx) + 1) * sizeof(char)
+ );
+ values[1] = (char *) palloc(
+ (strlen(drv_set2[call_cntr].short_name) + 1) * sizeof(char)
+ );
+ values[2] = (char *) palloc(
+ (strlen(drv_set2[call_cntr].long_name) + 1) * sizeof(char)
+ );
+ values[3] = (char *) palloc(
+ (strlen(drv_set2[call_cntr].create_options) + 1) * sizeof(char)
+ );
+
+ snprintf(
+ values[0],
+ (snprintf(NULL, 0, "%d", drv_set2[call_cntr].idx) + 1) * sizeof(char),
+ "%d",
+ drv_set2[call_cntr].idx
+ );
+ snprintf(
+ values[1],
+ (strlen(drv_set2[call_cntr].short_name) + 1) * sizeof(char),
+ "%s",
+ drv_set2[call_cntr].short_name
+ );
+ snprintf(
+ values[2],
+ (strlen(drv_set2[call_cntr].long_name) + 1) * sizeof(char),
+ "%s",
+ drv_set2[call_cntr].long_name
+ );
+ snprintf(
+ values[3],
+ (strlen(drv_set2[call_cntr].create_options) + 1) * sizeof(char),
+ "%s",
+ drv_set2[call_cntr].create_options
+ );
+
+ POSTGIS_RT_DEBUGF(4, "Result %d, Index %s", call_cntr, values[0]);
+ POSTGIS_RT_DEBUGF(4, "Result %d, Short Name %s", call_cntr, values[1]);
+ POSTGIS_RT_DEBUGF(4, "Result %d, Full Name %s", call_cntr, values[2]);
+ POSTGIS_RT_DEBUGF(5, "Result %d, Create Options %s", call_cntr, values[3]);
+
+ /* build a tuple */
+ tuple = BuildTupleFromCStrings(attinmeta, values);
+
+ /* make the tuple into a datum */
+ result = HeapTupleGetDatum(tuple);
+
+ /* clean up (this is not really necessary) */
+ pfree(values[3]);
+ pfree(values[2]);
+ pfree(values[1]);
+ pfree(values[0]);
+ pfree(values);
+
+ SRF_RETURN_NEXT(funcctx, result);
+ }
+ /* do when there is no more left */
+ else {
+ pfree(drv_set2);
+ SRF_RETURN_DONE(funcctx);
+ }
+}
+
/* ---------------------------------------------------------------- */
/* Memory allocation / error reporting hooks */
/* ---------------------------------------------------------------- */
AS $$ SELECT st_reclass($1, ROW(1, $2, $3, NULL)) $$
LANGUAGE 'SQL' IMMUTABLE STRICT;
+-----------------------------------------------------------------------
+-- ST_AsGDALRaster and supporting functions
+-----------------------------------------------------------------------
+-- returns set of available and usable GDAL drivers
+CREATE OR REPLACE FUNCTION st_gdaldrivers(OUT idx int, OUT short_name text, OUT long_name text, OUT create_options text)
+ RETURNS SETOF record
+ AS 'MODULE_PATHNAME', 'RASTER_getGDALDrivers'
+ LANGUAGE 'C' IMMUTABLE STRICT;
+
+-- return the srtext or proj4text of an srid
+CREATE OR REPLACE FUNCTION st_srtext(q raster)
+ RETURNS text
+ AS $$
+ DECLARE
+ q_srid int;
+ q_srs text;
+ BEGIN
+ q_srid := st_srid($1);
+
+ IF q_srid != -1 THEN
+ SELECT
+ CASE
+ WHEN length(srtext) > 0
+ THEN srtext
+ WHEN length(proj4text) > 0
+ THEN proj4text
+ ELSE NULL
+ END
+ INTO q_srs
+ FROM spatial_ref_sys
+ WHERE srid = q_srid;
+
+ IF NOT FOUND THEN
+ q_srs := NULL;
+ END IF;
+ ELSE
+ q_srs := NULL;
+ END IF;
+
+ RETURN q_srs;
+ END;
+ $$ LANGUAGE 'plpgsql' IMMUTABLE STRICT;
+
+-- Cannot be strict as "options" and "srs" can be NULL
+CREATE OR REPLACE FUNCTION st_asgdalraster(rast raster, format text, options text[], srs text)
+ RETURNS bytea
+ AS 'MODULE_PATHNAME', 'RASTER_asGDALRaster'
+ LANGUAGE 'C' IMMUTABLE;
+
+CREATE OR REPLACE FUNCTION st_asgdalraster(rast raster, format text, options text[])
+ RETURNS bytea
+ AS $$ SELECT st_asgdalraster($1, $2, $3, st_srtext($1)) $$
+ LANGUAGE 'SQL' IMMUTABLE STRICT;
+
+CREATE OR REPLACE FUNCTION st_asgdalraster(rast raster, format text)
+ RETURNS bytea
+ AS $$ SELECT st_asgdalraster($1, $2, NULL, st_srtext($1)) $$
+ LANGUAGE 'SQL' IMMUTABLE STRICT;
+
-----------------------------------------------------------------------
-- MapAlgebra
-----------------------------------------------------------------------
--- /dev/null
+#! /usr/bin/env python
+#
+# $Id$
+#
+# A simple utility for passing parameters to ST_AsGDALRaster
+# for a GDAL raster output
+# This utility is handy for debugging purposes.
+#
+# Example:
+# rtgdalraster.py -d "dbname=postgres" -r "(SELECT rast FROM pele LIMIT 1)" -f JPEG -o pele.jpg
+#
+#
+###############################################################################
+# Copyright (C) 2009 Mateusz Loskot <mateusz@loskot.net>
+#
+# This program is free software; you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation; either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program; if not, write to the Free Software
+# Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+#
+###############################################################################
+from optparse import OptionParser
+import sys
+import os
+import psycopg2
+
+def logit(msg):
+ if VERBOSE is True:
+ sys.stderr.write("LOG - " + str(msg) + "\n")
+
+###############################################################################
+try:
+ prs = OptionParser(version="%prog",
+ usage="%prog -d <DB> -r <RASTER> -o <FILE> [-f <FORMAT>] [-c <OPTIONS>]",
+ description="Output raster in a GDAL raster format")
+ prs.add_option("-d", "--db", dest="db", action="store", default=None,
+ help="PostgreSQL database connection string, required")
+ prs.add_option("-r", "--raster", dest="raster", action="store", default=None,
+ help="sql expression to create or access a existing raster, required")
+ prs.add_option("-o", "--output", dest="output", action="store", default=None,
+ help="file to put the output of ST_AsGDALRaster if successful, required")
+ prs.add_option("-f", "--format", dest="format", action="store", default="GTiff",
+ help="format of GDAL raster output, optional, default=GTiff")
+ prs.add_option("-c", "--config", dest="config", action="store", default=None,
+ help="comma separated list of GDAL raster creation options, optional")
+ prs.add_option("-v", "--verbose", dest="verbose", action="store_true", default=False,
+ help="be excessively verbose and useful for debugging")
+
+ (opts, args) = prs.parse_args()
+
+ if opts.db is None:
+ prs.error("use -d option to specify database connection string")
+ if opts.raster is None:
+ prs.error("use -r option to specify a sql expression for a raster")
+ if opts.output is None:
+ prs.error("use -o option to specify raster output file")
+
+ if opts.config is not None:
+ opts.cfg = opts.config.split(',')
+ else:
+ opts.cfg = None
+
+ global VERBOSE
+ VERBOSE = opts.verbose
+
+ conn = psycopg2.connect(opts.db)
+ logit("Connected to %s" % opts.db)
+ logit("Raster expression is %s" % opts.raster)
+
+ cur = conn.cursor()
+ sql = "SELECT ST_AsGDALRaster(%s, %%s, %%s::text[], NULL::text)" % (opts.raster)
+ cur.execute(sql, (opts.format, opts.cfg))
+ logit("Number of rows: %i" % cur.rowcount)
+ rec = cur.fetchone()[0];
+ logit("size of raster (bytes): %i" % len(rec))
+
+ fh = open(opts.output, 'wb', -1)
+ fh.write(rec);
+ fh.flush();
+ fh.close();
+ logit("size of %s (bytes): %i" % (opts.output, os.stat(opts.output)[6]));
+
+ cur.close();
+ conn.close();
+
+ print "raster outputted to %s" % opts.output;
+
+except Exception, e:
+ print "ERROR: ", e
rt_raster_destroy(raster);
}
+struct rt_gdaldriver_t {
+ int idx;
+ char *short_name;
+ char *long_name;
+ char *create_options;
+};
+static void testGDALDrivers() {
+ int i;
+ uint32_t size;
+ rt_gdaldriver drv;
+
+ drv = (rt_gdaldriver) rt_raster_gdal_drivers(&size);
+ /*printf("size: %d\n", size);*/
+ CHECK(size);
+
+ for (i = 0; i < size; i++) {
+ /*printf("gdal_driver: %s\n", drv[i].short_name);*/
+ CHECK(drv[i].short_name);
+ }
+}
+
+static void testRasterToGDAL(rt_raster raster) {
+ uint32_t bandNums[] = {1};
+ int lenBandNums = 1;
+ uint64_t gdalSize;
+ rt_raster rast;
+ uint8_t *rtn = NULL;
+
+ rast = rt_raster_from_band(raster, bandNums, lenBandNums);
+ assert(rast);
+
+ rtn = rt_raster_to_gdal(rast, NULL, "GTiff", NULL, &gdalSize);
+ /*printf("gdalSize: %d\n", (int) gdalSize);*/
+ CHECK(gdalSize);
+
+ /*
+ FILE *fh = NULL;
+ fh = fopen("/tmp/out.tif", "w");
+ fwrite(rtn, sizeof(uint8_t), gdalSize, fh);
+ fclose(fh);
+ */
+
+ rt_raster_destroy(rast);
+}
+
int
main()
{
testBandReclass();
printf("Successfully tested rt_band_reclass\n");
+ printf("Testing rt_raster_to_gdal\n");
+ testRasterToGDAL(raster);
+ printf("Successfully tested rt_raster_to_gdal\n");
+
+ printf("Testing rt_raster_gdal_drivers\n");
+ testGDALDrivers();
+ printf("Successfully tested rt_raster_gdal_drivers\n");
deepRelease(raster);
rt_box2d.sql \
rt_addband.sql \
rt_band.sql \
+ rt_asgdalraster.sql \
$(NULL)
TEST_PROPS = \
--- /dev/null
+SELECT md5(
+ ST_AsGDALRaster(
+ ST_AddBand(ST_MakeEmptyRaster(200, 200, 10, 10, 2, 2, 0, 0, -1), 1, '64BF', 123.4567, NULL),
+ 'GTiff'
+ )
+);
+SELECT md5(
+ ST_AsGDALRaster(
+ ST_AddBand(
+ ST_AddBand(
+ ST_AddBand(
+ ST_MakeEmptyRaster(200, 200, 10, 10, 2, 2, 0, 0,-1)
+ , 1, '64BF', 1234.5678, NULL
+ )
+ , '64BF', 987.654321, NULL
+ )
+ , '64BF', 9876.54321, NULL
+ ),
+ 'GTiff'
+ )
+);
+SELECT md5(
+ ST_AsGDALRaster(
+ ST_AddBand(
+ ST_AddBand(
+ ST_AddBand(
+ ST_MakeEmptyRaster(200, 200, 10, 10, 2, 2, 0, 0, -1)
+ , 1, '64BF', 1234.5678, -9999
+ )
+ , '64BF', 987.654321, -9999
+ )
+ , '64BF', 9876.54321, -9999
+ ),
+ 'GTiff'
+ )
+);
+SELECT md5(
+ ST_AsGDALRaster(
+ ST_AddBand(ST_MakeEmptyRaster(200, 200, 10, 10, 2, 2, 0, 0, -1), 1, '8BSI', 123, NULL),
+ 'PNG'
+ )
+);
+SELECT md5(
+ ST_AsGDALRaster(
+ ST_AddBand(ST_MakeEmptyRaster(200, 200, 10, 10, 2, 2, 0, 0, -1), 1, '8BUI', 123, NULL),
+ 'PNG'
+ )
+);
+SELECT md5(
+ ST_AsGDALRaster(
+ ST_AddBand(ST_MakeEmptyRaster(200, 200, 10, 10, 2, 2, 0, 0, -1), 1, '8BSI', -123, NULL),
+ 'PNG'
+ )
+);
+SELECT md5(
+ ST_AsGDALRaster(
+ ST_AddBand(ST_MakeEmptyRaster(200, 200, 10, 10, 2, 2, 0, 0, -1), 1, '8BUI', 254, NULL),
+ 'PNG'
+ )
+);
+SELECT md5(
+ ST_AsGDALRaster(
+ ST_AddBand(
+ ST_AddBand(
+ ST_AddBand(
+ ST_MakeEmptyRaster(200, 200, 10, 10, 2, 2, 0, 0, -1)
+ , 1, '8BSI', 1, -1
+ )
+ , 2, '8BSI', 11, -1
+ )
+ , 3, '8BSI', 111, -1
+ ),
+ 'PNG',
+ ARRAY['ZLEVEL=1']
+ )
+);
+SELECT md5(
+ ST_AsGDALRaster(
+ ST_AddBand(
+ ST_AddBand(
+ ST_AddBand(
+ ST_MakeEmptyRaster(200, 200, 10, 10, 2, 2, 0, 0, -1)
+ , 1, '8BSI', 1, -1
+ )
+ , 2, '8BSI', 11, -1
+ )
+ , 3, '8BSI', 111, -1
+ ),
+ 'PNG',
+ ARRAY['ZLEVEL=9']
+ )
+);
+SELECT md5(
+ ST_AsGDALRaster(
+ ST_AddBand(ST_MakeEmptyRaster(200, 200, 10, 10, 2, 2, 0, 0, -1), 1, '8BSI', 123, NULL),
+ 'JPEG'
+ )
+);
+SELECT md5(
+ ST_AsGDALRaster(
+ ST_AddBand(ST_MakeEmptyRaster(200, 200, 10, 10, 2, 2, 0, 0, -1), 1, '8BUI', 123, NULL),
+ 'JPEG'
+ )
+);
+SELECT md5(
+ ST_AsGDALRaster(
+ ST_AddBand(ST_MakeEmptyRaster(200, 200, 10, 10, 2, 2, 0, 0, -1), 1, '8BSI', -123, NULL),
+ 'JPEG'
+ )
+);
+SELECT md5(
+ ST_AsGDALRaster(
+ ST_AddBand(ST_MakeEmptyRaster(200, 200, 10, 10, 2, 2, 0, 0, -1), 1, '8BUI', 254, NULL),
+ 'JPEG'
+ )
+);
+SELECT md5(
+ ST_AsGDALRaster(
+ ST_AddBand(
+ ST_AddBand(
+ ST_AddBand(
+ ST_MakeEmptyRaster(200, 200, 10, 10, 2, 2, 0, 0, -1)
+ , 1, '8BSI', 1, -1
+ )
+ , 2, '8BSI', 11, -1
+ )
+ , 3, '8BSI', 111, -1
+ ),
+ 'JPEG'
+ )
+);
+SELECT md5(
+ ST_AsGDALRaster(
+ ST_AddBand(
+ ST_AddBand(
+ ST_AddBand(
+ ST_MakeEmptyRaster(200, 200, 10, 10, 2, 2, 0, 0, -1)
+ , 1, '8BSI', 1, -1
+ )
+ , 2, '8BSI', 11, -1
+ )
+ , 3, '8BSI', 111, -1
+ ),
+ 'JPEG',
+ ARRAY['QUALITY=90','PROGRESSIVE=ON']
+ )
+);
--- /dev/null
+2fd5de5f94416f2999fc24fbc3974a4b
+4af1ca04adf045690fb60764cc93de75
+81de10108dd8eede75ef47c8f6f09bc5
+55279950e29968bcf36b2c11ce8bf88b
+55279950e29968bcf36b2c11ce8bf88b
+4ae84e77caa1a35aa4da233e83feb746
+24188762b5745acda4aa7b92776e7280
+da3c4e827b617f9d7cd83d5759f1d781
+0d6d88030359474c2e86280a65e5e90a
+f0de16a21de438249ec0bbd28eafe63c
+f0de16a21de438249ec0bbd28eafe63c
+ed4e4d60d2ccc82c4597e80c1f97272f
+83b6012757444ff7786b5e8473e6cea3
+b699be9f7045656edbb14ecec1c75b5a
+96a96f9a6d31544dfc5ba5cd25f40934
projects / postgis / commitdiff