Commit a10cbcfa authored by Bborie Park's avatar Bborie Park

Add ST_AsGDALRaster function and helper functions ST_GDALDrivers and ST_srtext

	- added rt_raster_to_gdal, rt_raster_gdal_drivers and rt_raster_to_gdal_mem functions to rt_core/rt_api.c and rt_api.h
	- added test cases to test/core/testapi.c
	- added RASTER_asGDALRaster and RASTER_getGDALDrivers to rt_pg/rt_pg.c
	- added SQL functions
	- added regression tests

Associated ticket is #901


git-svn-id: http://svn.osgeo.org/postgis/trunk@7155 b70326c6-7e19-0410-871a-916f4a2858ee
parent 6eb30cd2
This diff is collapsed.
......@@ -91,6 +91,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"
/**
......@@ -845,6 +847,41 @@ rt_raster rt_raster_from_band(rt_raster raster, uint32_t *bandNums,
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 -------------------------------------------------------*/
......
......@@ -72,6 +72,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...
......@@ -205,6 +206,9 @@ Datum RASTER_quantile(PG_FUNCTION_ARGS);
/* 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
......@@ -387,6 +391,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;
}
PG_FUNCTION_INFO_V1(RASTER_lib_version);
Datum RASTER_lib_version(PG_FUNCTION_ARGS)
{
......@@ -4117,6 +4153,329 @@ 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;
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 */
/* ---------------------------------------------------------------- */
......
......@@ -1022,6 +1022,65 @@ CREATE OR REPLACE FUNCTION st_reclass(rast raster, reclassexpr text, pixeltype t
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
-----------------------------------------------------------------------
......
#! /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
......@@ -1261,6 +1261,51 @@ static void testBandReclass() {
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()
{
......@@ -1583,6 +1628,13 @@ 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);
......
......@@ -42,6 +42,7 @@ TEST_FUNC = \
rt_box2d.sql \
rt_addband.sql \
rt_band.sql \
rt_asgdalraster.sql \
$(NULL)
TEST_PROPS = \
......
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),