diff -rupN postgis-old/raster/rt_core/rt_api.c postgis-new/raster/rt_core/rt_api.c --- postgis-old/raster/rt_core/rt_api.c 2011-04-29 15:54:37.000000000 -0700 +++ postgis-new/raster/rt_core/rt_api.c 2011-04-29 15:54:50.000000000 -0700 @@ -4296,3 +4296,397 @@ rt_band_reclass(rt_context ctx, rt_band return band; } + +/** + * Return formatted GDAL raster from raster + * + * @param ctx : context, for thread safety + * @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_context ctx, 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 != ctx); + 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(ctx, raster, srs, &src_drv); + if (NULL == src_ds) { + ctx->err("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) { + ctx->err("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) { + ctx->err("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) { + ctx->err("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 ctx : context, for thread safety + * @param drv_count : number of GDAL drivers available + * @return set of "gdaldriver" values of available GDAL drivers + * + */ +rt_gdaldriver +rt_raster_gdal_drivers(rt_context ctx, 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; + + assert(NULL != ctx); + + GDALAllRegister(); + count = GDALGetDriverCount(); + rtn = (rt_gdaldriver) ctx->alloc(count * sizeof(struct rt_gdaldriver_t)); + if (NULL == rtn) { + ctx->err("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 *) ctx->alloc(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 *) ctx->alloc(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 *) ctx->alloc(txt_len); + memcpy(rtn[j].create_options, txt, txt_len); + + j++; + } + + /* free unused memory */ + ctx->realloc(rtn, j * sizeof(struct rt_gdaldriver_t)); + *drv_count = j; + + return rtn; +} + + +/** + * Return GDAL datasource using GDAL MEM driver from raster + * @param ctx : context, for thread safety + * @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_context ctx, 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 != ctx); + 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) { + ctx->err("rt_raster_to_gdal_mem: Unable to load the MEM GDAL driver\n"); + return 0; + } + } + + ds = GDALCreate(drv, "", rt_raster_get_width(ctx, raster), + rt_raster_get_height(ctx, raster), 0, GDT_Byte, NULL); + if (NULL == ds) { + ctx->err("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(ctx, raster); + gt[1] = rt_raster_get_x_scale(ctx, raster); + gt[2] = rt_raster_get_x_skew(ctx, raster); + gt[3] = rt_raster_get_y_offset(ctx, raster); + gt[4] = rt_raster_get_y_skew(ctx, raster); + gt[5] = rt_raster_get_y_scale(ctx, raster); + cplerr = GDALSetGeoTransform(ds, gt); + if (cplerr != CE_None) { + ctx->err("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) { + ctx->err("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(ctx, raster); + for (i = 0; i < numBands; i++) { + rtband = rt_raster_get_band(ctx, raster, i); + if (0 == rtband) { + ctx->err("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(ctx, 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: + ctx->warn("rt_raster_to_gdal_mem: Unknown pixel type for band\n"); + gdal_pt = GDT_Unknown; + break; + } + + pVoid = rt_band_get_data(ctx, rtband); + RASTER_DEBUGF(4, "Band data is at pos %p", pVoid); + + pszDataPointer = (char *) ctx->alloc(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 */ + ctx->dealloc(pszDataPointer); + + if (GDALAddBand(ds, gdal_pt, apszOptions) == CE_Failure) { + ctx->err("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) { + ctx->err("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) { + ctx->err("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(ctx, rtband) != FALSE) { + nodata = rt_band_get_nodata(ctx, rtband); + if (GDALSetRasterNoDataValue(band, nodata) != CE_None) + ctx->warn("rt_raster_to_gdal_mem: Couldn't set nodata value for band\n"); + } + } + + return ds; +} diff -rupN postgis-old/raster/rt_core/rt_api.h postgis-new/raster/rt_core/rt_api.h --- postgis-old/raster/rt_core/rt_api.h 2011-04-29 15:54:37.000000000 -0700 +++ postgis-new/raster/rt_core/rt_api.h 2011-04-29 15:54:50.000000000 -0700 @@ -93,6 +93,8 @@ #include "gdal_frmts.h" #include "gdal.h" #include "ogr_api.h" +#include "ogr_srs_api.h" +#include "cpl_vsi.h" #include "../../postgis_config.h" @@ -798,6 +800,44 @@ rt_band rt_band_reclass(rt_context ctx, rt_pixtype pixtype, uint32_t hasnodata, double nodataval, rt_reclassexpr *exprset, int exprcount); +/** + * Return formatted GDAL raster from raster + * + * @param ctx : context, for thread safety + * @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_context ctx, rt_raster raster, char *srs, + char *format, char **options, uint64_t *gdalsize); + +/** + * Returns a set of available GDAL drivers + * + * @param ctx : context, for thread safety + * @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(rt_context ctx, uint32_t *drv_count); + +/** + * Return GDAL datasource using GDAL MEM driver from raster + * @param ctx : context, for thread safety + * @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_context ctx, rt_raster raster, + char *srs, GDALDriverH *rtn_drv); + /*- utilities -------------------------------------------------------*/ /* Set of functions to clamp double to int of different size diff -rupN postgis-old/raster/rt_pg/rt_pg.c postgis-new/raster/rt_pg/rt_pg.c --- postgis-old/raster/rt_pg/rt_pg.c 2011-04-29 15:54:37.000000000 -0700 +++ postgis-new/raster/rt_pg/rt_pg.c 2011-04-29 15:54:50.000000000 -0700 @@ -76,6 +76,7 @@ static char *strtoupper(char *str); 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... @@ -177,6 +178,10 @@ Datum RASTER_summaryStats(PG_FUNCTION_AR /* 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 */ /* --------------------------------------------------------------------------- Name : replace - Search & replace a substring by another one. @@ -370,6 +375,38 @@ removespaces(char *str) { 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; +} + static void * rt_pgalloc(size_t size) { @@ -3598,6 +3635,347 @@ Datum RASTER_reclass(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; + rt_context ctx = get_rt_context(fcinfo); + + 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(ctx, pgraster); + if (!raster) { + elog(ERROR, "RASTER_asGDALRaster: Could not deserialize raster"); + rt_context_destroy(ctx); + PG_RETURN_NULL(); + } + + /* format is required */ + if (PG_ARGISNULL(1)) { + elog(ERROR, "RASTER_asGDALRaster: Format must be provided"); + rt_raster_destroy(ctx, raster); + rt_context_destroy(ctx); + 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(ctx, 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(ctx, raster); + rt_context_destroy(ctx); + + 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; + + rt_context ctx; + 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); + + /** + * Create a new context. We don't call get_rt_context, because this + * function changes the memory context to the one pointed by + * fcinfo->flinfo->fn_mcxt, and we want to keep the current one + */ + ctx = rt_context_new(rt_pgalloc, rt_pgrealloc, rt_pgfree); + rt_context_set_message_handlers(ctx, rt_pgerr, rt_pgwarn, rt_pginfo); + + drv_set = rt_raster_gdal_drivers(ctx, &drv_count); + if (NULL == drv_set || !drv_count) { + elog(ERROR, "RASTER_gdal_drivers: No GDAL drivers found"); + PG_RETURN_NULL(); + } + + 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 *)); + + POSTGIS_RT_DEBUGF(3, "idx length %d", snprintf(NULL, 0, "%d", drv_set2[call_cntr].idx)); + + 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(4, "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 */ /* ---------------------------------------------------------------- */ diff -rupN postgis-old/raster/rt_pg/rtpostgis.sql.in.c postgis-new/raster/rt_pg/rtpostgis.sql.in.c --- postgis-old/raster/rt_pg/rtpostgis.sql.in.c 2011-04-29 15:54:37.000000000 -0700 +++ postgis-new/raster/rt_pg/rtpostgis.sql.in.c 2011-04-29 15:54:50.000000000 -0700 @@ -437,6 +437,65 @@ CREATE OR REPLACE FUNCTION st_reclass(ra 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 ----------------------------------------------------------------------- -- This function can not be STRICT, because nodatavalueexpr can be NULL (could be just '' though) diff -rupN postgis-old/raster/scripts/python/rtgdalraster.py postgis-new/raster/scripts/python/rtgdalraster.py --- postgis-old/raster/scripts/python/rtgdalraster.py 1969-12-31 16:00:00.000000000 -0800 +++ postgis-new/raster/scripts/python/rtgdalraster.py 2011-04-29 15:54:50.000000000 -0700 @@ -0,0 +1,98 @@ +#! /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 +# +# 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 -r -o [-f ] [-c ]", + 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 diff -rupN postgis-old/raster/test/core/testapi.c postgis-new/raster/test/core/testapi.c --- postgis-old/raster/test/core/testapi.c 2011-04-29 15:54:37.000000000 -0700 +++ postgis-new/raster/test/core/testapi.c 2011-04-29 15:54:50.000000000 -0700 @@ -1211,6 +1211,51 @@ static void testBandReclass(rt_context c rt_raster_destroy(ctx, raster); } +struct rt_gdaldriver_t { + int idx; + char *short_name; + char *long_name; + char *create_options; +}; +static void testGDALDrivers(rt_context ctx) { + int i; + uint32_t size; + rt_gdaldriver drv; + + drv = (rt_gdaldriver) rt_raster_gdal_drivers(ctx, &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_context ctx, 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(ctx, raster, bandNums, lenBandNums); + assert(rast); + + rtn = rt_raster_to_gdal(ctx, 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(ctx, rast); +} + int main() { @@ -1537,6 +1582,14 @@ main() testBandReclass(ctx); printf("Successfully tested rt_band_reclass\n"); + printf("Testing rt_raster_to_gdal\n"); + testRasterToGDAL(ctx, raster); + printf("Successfully tested rt_raster_to_gdal\n"); + + printf("Testing rt_raster_gdal_drivers\n"); + testGDALDrivers(ctx); + printf("Successfully tested rt_raster_gdal_drivers\n"); + deepRelease(ctx, raster); rt_context_destroy(ctx); diff -rupN postgis-old/raster/test/regress/Makefile.in postgis-new/raster/test/regress/Makefile.in --- postgis-old/raster/test/regress/Makefile.in 2011-04-29 15:54:37.000000000 -0700 +++ postgis-new/raster/test/regress/Makefile.in 2011-04-29 15:54:50.000000000 -0700 @@ -43,6 +43,7 @@ TEST_FUNC = \ rt_addband.sql \ rt_band.sql \ rt_reclass.sql \ + rt_asgdalraster.sql \ $(NULL) TEST_PROPS = \ diff -rupN postgis-old/raster/test/regress/rt_asgdalraster.sql postgis-new/raster/test/regress/rt_asgdalraster.sql --- postgis-old/raster/test/regress/rt_asgdalraster.sql 1969-12-31 16:00:00.000000000 -0800 +++ postgis-new/raster/test/regress/rt_asgdalraster.sql 2011-04-29 15:54:50.000000000 -0700 @@ -0,0 +1,147 @@ +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'] + ) +); diff -rupN postgis-old/raster/test/regress/rt_asgdalraster_expected postgis-new/raster/test/regress/rt_asgdalraster_expected --- postgis-old/raster/test/regress/rt_asgdalraster_expected 1969-12-31 16:00:00.000000000 -0800 +++ postgis-new/raster/test/regress/rt_asgdalraster_expected 2011-04-29 15:54:50.000000000 -0700 @@ -0,0 +1,15 @@ +2fd5de5f94416f2999fc24fbc3974a4b +4af1ca04adf045690fb60764cc93de75 +81de10108dd8eede75ef47c8f6f09bc5 +55279950e29968bcf36b2c11ce8bf88b +55279950e29968bcf36b2c11ce8bf88b +4ae84e77caa1a35aa4da233e83feb746 +24188762b5745acda4aa7b92776e7280 +01714fc4d54e679372bd4616585851a1 +5faf1a59559d5d2950629fb2bf276911 +f0de16a21de438249ec0bbd28eafe63c +f0de16a21de438249ec0bbd28eafe63c +ed4e4d60d2ccc82c4597e80c1f97272f +83b6012757444ff7786b5e8473e6cea3 +b699be9f7045656edbb14ecec1c75b5a +96a96f9a6d31544dfc5ba5cd25f40934