36 #include <executor/spi.h>
37 #include <utils/lsyscache.h>
38 #include <utils/array.h>
39 #include <utils/builtins.h>
40 #include <catalog/pg_type.h>
41 #include <executor/executor.h>
43 #include "../../postgis_config.h"
44 #include "lwgeom_pg.h"
80 #if defined(__clang__)
81 # pragma clang diagnostic push
82 # pragma clang diagnostic ignored "-Wgnu-variable-sized-type-not-at-end"
89 #if POSTGIS_PGSQL_VERSION < 120
94 FunctionCallInfoBaseData fcinfo;
95 char fcinfo_data[SizeForFunctionCallInfo(FUNC_MAX_ARGS)];
97 FunctionCallInfo ufc_info;
101 #if defined(__clang__)
102 # pragma clang diagnostic pop
134 elog(ERROR,
"rtpg_nmapalgebra_arg_init: Could not allocate memory for arguments");
158 #if POSTGIS_PGSQL_VERSION >= 120
170 if (arg->
raster != NULL) {
187 if( arg->
mask != NULL )
210 if (arg == NULL || array == NULL) {
211 elog(ERROR,
"rtpg_nmapalgebra_rastbandarg_process: NULL values not permitted for parameters");
215 etype = ARR_ELEMTYPE(array);
216 get_typlenbyvalalign(etype, &typlen, &typbyval, &typalign);
221 typlen, typbyval, typalign,
226 elog(ERROR,
"rtpg_nmapalgebra_rastbandarg_process: Invalid argument for rastbandarg");
244 arg->
nband == NULL ||
248 elog(ERROR,
"rtpg_nmapalgebra_rastbandarg_process: Could not allocate memory for processing rastbandarg");
257 for (i = 0; i < n; i++) {
272 tup = (HeapTupleHeader) DatumGetPointer(e[i]);
274 elog(ERROR,
"rtpg_nmapalgebra_rastbandarg_process: Invalid argument for rastbandarg at index %d", i);
280 tupv = GetAttributeByName(tup,
"rast", &isnull);
282 elog(NOTICE,
"First argument (nband) of rastbandarg at index %d is NULL. Assuming NULL raster", i);
296 for (j = 0; j < i; j++) {
308 if (arg->
raster[i] == NULL) {
309 elog(ERROR,
"rtpg_nmapalgebra_rastbandarg_process: Could not deserialize raster at index %d", i);
325 tupv = GetAttributeByName(tup,
"nband", &isnull);
328 elog(NOTICE,
"First argument (nband) of rastbandarg at index %d is NULL. Assuming nband = %d", i,
nband);
331 nband = DatumGetInt32(tupv);
334 elog(ERROR,
"rtpg_nmapalgebra_rastbandarg_process: Band number provided for rastbandarg at index %d must be greater than zero (1-based)", i);
358 arg->
nband == NULL ||
361 elog(ERROR,
"rtpg_nmapalgebra_rastbandarg_process: Could not reallocate memory for processed rastbandarg");
376 double *
value,
int *nodata
384 ArrayType *mdValues = NULL;
385 Datum *_values = NULL;
386 bool *_nodata = NULL;
388 ArrayType *mdPos = NULL;
397 int lbound[3] = {1, 1, 1};
398 Datum datum = (Datum) NULL;
412 if (_values == NULL || _nodata == NULL) {
413 elog(ERROR,
"rtpg_nmapalgebra_callback: Could not allocate memory for values array");
420 for (z = 0; z < arg->
rasters; z++) {
422 for (
y = 0;
y < arg->
rows;
y++) {
428 _nodata[i] = (bool) arg->
nodata[z][
y][
x];
430 _values[i] = Float8GetDatum(arg->
values[z][
y][
x]);
432 _values[i] = (Datum) NULL;
440 get_typlenbyvalalign(FLOAT8OID, &typlen, &typbyval, &typalign);
443 mdValues = construct_md_array(
447 typlen, typbyval, typalign
452 _pos = palloc(
sizeof(Datum) * (arg->
rasters + 1) * 2);
453 _null = palloc(
sizeof(
bool) * (arg->
rasters + 1) * 2);
454 if (_pos == NULL || _null == NULL) {
456 elog(ERROR,
"rtpg_nmapalgebra_callback: Could not allocate memory for position array");
459 memset(_null, 0,
sizeof(
bool) * (arg->
rasters + 1) * 2);
468 for (z = 0; z < arg->
rasters; z++) {
477 get_typlenbyvalalign(INT4OID, &typlen, &typbyval, &typalign);
485 mdPos = construct_md_array(
489 typlen, typbyval, typalign
494 #if POSTGIS_PGSQL_VERSION < 120
495 callback->
ufc_info.arg[0] = PointerGetDatum(mdValues);
496 callback->
ufc_info.arg[1] = PointerGetDatum(mdPos);
498 callback->
ufc_info->args[0].value = PointerGetDatum(mdValues);
499 callback->
ufc_info->args[1].value = PointerGetDatum(mdPos);
503 #if POSTGIS_PGSQL_VERSION < 120
504 datum = FunctionCallInvoke(&(callback->
ufc_info));
506 datum = FunctionCallInvoke(callback->
ufc_info);
512 #if POSTGIS_PGSQL_VERSION < 120
520 *
value = DatumGetFloat8(datum);
523 *
value = (double) DatumGetFloat4(datum);
526 *
value = (double) DatumGetInt32(datum);
529 *
value = (double) DatumGetInt16(datum);
547 ArrayType *maskArray;
578 elog(ERROR,
"RASTER_nMapAlgebra: Could not initialize argument structure");
585 elog(ERROR,
"RASTER_nMapAlgebra: Could not process rastbandarg");
589 POSTGIS_RT_DEBUGF(4,
"allnull, allempty, noband = %d, %d, %d", allnull, allempty, noband);
593 elog(NOTICE,
"All input rasters are NULL. Returning NULL");
599 if (!PG_ARGISNULL(2)) {
606 elog(ERROR,
"RASTER_nMapAlgebra: Invalid pixel type: %s", pixtypename);
612 if (!PG_ARGISNULL(3)){
613 arg->
distance[0] = PG_GETARG_INT32(3);
618 if (!PG_ARGISNULL(4)){
619 arg->
distance[1] = PG_GETARG_INT32(4);
625 elog(ERROR,
"RASTER_nMapAlgebra: Distance for X and Y axis must be greater than or equal to zero");
630 if (!PG_ARGISNULL(5)) {
638 if (PG_ARGISNULL(6)) {
639 elog(NOTICE,
"Custom extent is NULL. Returning NULL");
650 elog(ERROR,
"RASTER_nMapAlgebra: Could not deserialize custom extent");
654 elog(NOTICE,
"Custom extent is an empty raster. Returning empty raster");
659 elog(ERROR,
"RASTER_nMapAlgebra: Could not create empty raster");
665 if (!pgraster) PG_RETURN_NULL();
667 SET_VARSIZE(pgraster, pgraster->
size);
668 PG_RETURN_POINTER(pgraster);
674 if( PG_ARGISNULL(7) ){
679 maskArray = PG_GETARG_ARRAYTYPE_P(7);
680 etype = ARR_ELEMTYPE(maskArray);
681 get_typlenbyvalalign(etype,&typlen,&typbyval,&typalign);
689 elog(ERROR,
"RASTER_nMapAlgebra: Mask data type must be FLOAT8 or FLOAT4");
693 ndims = ARR_NDIM(maskArray);
696 elog(ERROR,
"RASTER_nMapAlgebra: Mask Must be a 2D array");
701 maskDims = ARR_DIMS(maskArray);
703 if (maskDims[0] % 2 == 0 || maskDims[1] % 2 == 0) {
704 elog(ERROR,
"RASTER_nMapAlgebra: Mask dimensions must be odd");
712 typlen, typbyval,typalign,
713 &maskElements,&maskNulls,&num
716 if (num < 1 || num != (maskDims[0] * maskDims[1])) {
721 elog(ERROR,
"RASTER_nMapAlgebra: Could not deconstruct new values array");
727 arg->
mask->
values = palloc(
sizeof(
double*)* maskDims[0]);
728 arg->
mask->
nodata = palloc(
sizeof(
int*)*maskDims[0]);
729 for (i = 0; i < maskDims[0]; i++) {
730 arg->
mask->
values[i] = (
double*) palloc(
sizeof(
double) * maskDims[1]);
731 arg->
mask->
nodata[i] = (
int*) palloc(
sizeof(
int) * maskDims[1]);
736 for (
y = 0;
y < maskDims[0];
y++) {
737 for (
x = 0;
x < maskDims[1];
x++) {
745 arg->
mask->
values[
y][
x] = (double) DatumGetFloat4(maskElements[i]);
749 arg->
mask->
values[
y][
x] = (double) DatumGetFloat8(maskElements[i]);
760 if (maskDims[0] == 1 && maskDims[1] == 1) {
771 if (PG_ARGISNULL(8) || !PG_GETARG_BOOL(8)) {
772 if (arg->
mask != NULL)
783 elog(NOTICE,
"All input rasters are empty. Returning empty raster");
788 elog(NOTICE,
"All input rasters do not have bands at indicated indexes. Returning empty raster");
796 elog(ERROR,
"RASTER_nMapAlgebra: Could not create empty raster");
802 if (!pgraster) PG_RETURN_NULL();
804 SET_VARSIZE(pgraster, pgraster->
size);
805 PG_RETURN_POINTER(pgraster);
809 if (!PG_ARGISNULL(1) || get_fn_expr_argtype(fcinfo->flinfo, 1) == REGPROCEDUREOID) {
828 get_func_result_type(
855 elog(ERROR,
"RASTER_nMapAlgebra: Function provided must return a double precision, float, int or smallint");
858 elog(ERROR,
"RASTER_nMapAlgebra: Function provided must return scalar (double precision, float, int, smallint)");
861 elog(ERROR,
"RASTER_nMapAlgebra: Function provided must have three input parameters");
864 elog(ERROR,
"RASTER_nMapAlgebra: Function provided must return double precision, not resultset");
871 elog(NOTICE,
"Function provided is VOLATILE. Unless required and for best performance, function should be IMMUTABLE or STABLE");
874 #if POSTGIS_PGSQL_VERSION < 120
892 if (!PG_ARGISNULL(9))
893 #if POSTGIS_PGSQL_VERSION < 120
902 #if POSTGIS_PGSQL_VERSION < 120
904 construct_empty_array(TEXTOID)
908 arg->
callback.
ufc_info->args[2].value = PointerGetDatum(construct_empty_array(TEXTOID));
913 #if POSTGIS_PGSQL_VERSION < 120
925 elog(ERROR,
"RASTER_nMapAlgebra: callbackfunc must be provided");
970 if (itrset == NULL) {
972 elog(ERROR,
"RASTER_nMapAlgebra: Could not allocate memory for iterator arguments");
1001 elog(ERROR,
"RASTER_nMapAlgebra: Could not run raster iterator function");
1015 SET_VARSIZE(pgraster, pgraster->
size);
1016 PG_RETURN_POINTER(pgraster);
1060 elog(ERROR,
"rtpg_nmapalgebraexpr_arg_init: Could not allocate memory for arguments");
1066 elog(ERROR,
"rtpg_nmapalgebraexpr_arg_init: Could not allocate memory for arg->bandarg");
1079 elog(ERROR,
"rtpg_nmapalgebraexpr_arg_init: Could not allocate memory for spi_argpos");
1110 double *
value,
int *nodata
1113 SPIPlanPtr plan = NULL;
1133 else if (arg->
nodata[0][0][0] && !arg->
nodata[1][0][0]) {
1143 else if (!arg->
nodata[0][0][0] && arg->
nodata[1][0][0]) {
1170 if (arg->
nodata[0][0][0]) {
1206 SPITupleTable *tuptable = NULL;
1209 bool isnull =
FALSE;
1214 memset(values, (Datum) NULL,
sizeof(Datum) * callback->
kw.
count);
1215 memset(nulls,
FALSE,
sizeof(
char) * callback->
kw.
count);
1220 for (i = 0; i < callback->
kw.
count; i++) {
1222 if (idx < 1)
continue;
1228 values[idx] = Int32GetDatum(arg->
src_pixel[0][0] + 1);
1232 values[idx] = Int32GetDatum(arg->
src_pixel[0][1] + 1);
1238 if (!arg->
nodata[0][0][0])
1239 values[idx] = Float8GetDatum(arg->
values[0][0][0]);
1246 values[idx] = Int32GetDatum(arg->
src_pixel[0][0] + 1);
1250 values[idx] = Int32GetDatum(arg->
src_pixel[0][1] + 1);
1256 if (!arg->
nodata[0][0][0])
1257 values[idx] = Float8GetDatum(arg->
values[0][0][0]);
1264 values[idx] = Int32GetDatum(arg->
src_pixel[1][0] + 1);
1268 values[idx] = Int32GetDatum(arg->
src_pixel[1][1] + 1);
1274 if (!arg->
nodata[1][0][0])
1275 values[idx] = Float8GetDatum(arg->
values[1][0][0]);
1285 err = SPI_execute_plan(plan, values, nulls,
TRUE, 1);
1286 if (err != SPI_OK_SELECT || SPI_tuptable == NULL || SPI_processed != 1) {
1287 elog(ERROR,
"rtpg_nmapalgebraexpr_callback: Unexpected error when running prepared statement %d",
id);
1292 tupdesc = SPI_tuptable->tupdesc;
1293 tuptable = SPI_tuptable;
1294 tuple = tuptable->vals[0];
1296 datum = SPI_getbinval(tuple, tupdesc, 1, &isnull);
1297 if (SPI_result == SPI_ERROR_NOATTRIBUTE) {
1298 if (SPI_tuptable) SPI_freetuptable(tuptable);
1299 elog(ERROR,
"rtpg_nmapalgebraexpr_callback: Could not get result of prepared statement %d",
id);
1304 *
value = DatumGetFloat8(datum);
1324 if (SPI_tuptable) SPI_freetuptable(tuptable);
1334 MemoryContext mainMemCtx = CurrentMemoryContext;
1337 uint16_t exprpos[3] = {1, 4, 5};
1351 SPITupleTable *tuptable = NULL;
1354 bool isnull =
FALSE;
1360 const int argkwcount = 12;
1376 if (PG_ARGISNULL(0))
1382 elog(ERROR,
"RASTER_nMapAlgebraExpr: Could not initialize argument structure");
1389 elog(ERROR,
"RASTER_nMapAlgebra: Could not process rastbandarg");
1393 POSTGIS_RT_DEBUGF(4,
"allnull, allempty, noband = %d, %d, %d", allnull, allempty, noband);
1397 elog(NOTICE,
"All input rasters are NULL. Returning NULL");
1409 if (!PG_ARGISNULL(2)) {
1416 elog(ERROR,
"RASTER_nMapAlgebraExpr: Invalid pixel type: %s", pixtypename);
1423 if (!PG_ARGISNULL(3)) {
1429 if (numraster < 2) {
1430 elog(NOTICE,
"CUSTOM extent type not supported. Defaulting to FIRST");
1434 elog(NOTICE,
"CUSTOM extent type not supported. Defaulting to INTERSECTION");
1438 else if (numraster < 2)
1444 if (!PG_ARGISNULL(6)) {
1452 elog(NOTICE,
"All input rasters are empty. Returning empty raster");
1457 elog(NOTICE,
"All input rasters do not have bands at indicated indexes. Returning empty raster");
1465 elog(ERROR,
"RASTER_nMapAlgebraExpr: Could not create empty raster");
1471 if (!pgraster) PG_RETURN_NULL();
1473 SET_VARSIZE(pgraster, pgraster->
size);
1474 PG_RETURN_POINTER(pgraster);
1478 if (SPI_connect() != SPI_OK_CONNECT) {
1480 elog(ERROR,
"RASTER_nMapAlgebraExpr: Could not connect to the SPI manager");
1496 char place[12] =
"$1";
1498 if (PG_ARGISNULL(exprpos[i]))
1504 for (j = 0, k = 1; j < argkwcount; j++) {
1516 sprintf(place,
"$%d", k);
1522 len = strlen(
"SELECT (") + strlen(expr) + strlen(
")::double precision");
1523 sql = (
char *) palloc(len + 1);
1527 elog(ERROR,
"RASTER_nMapAlgebraExpr: Could not allocate memory for expression parameter %d", exprpos[i]);
1531 memcpy(
sql,
"SELECT (", strlen(
"SELECT ("));
1532 memcpy(
sql + strlen(
"SELECT ("), expr, strlen(expr));
1533 memcpy(
sql + strlen(
"SELECT (") + strlen(expr),
")::double precision", strlen(
")::double precision"));
1542 if (argtype == NULL) {
1546 elog(ERROR,
"RASTER_nMapAlgebraExpr: Could not allocate memory for prepared plan argtypes of expression parameter %d", exprpos[i]);
1551 for (j = 0, k = 0; j < argkwcount; j++) {
1556 (strstr(argkw[j],
"[rast.x]") != NULL) ||
1557 (strstr(argkw[j],
"[rast.y]") != NULL) ||
1558 (strstr(argkw[j],
"[rast1.x]") != NULL) ||
1559 (strstr(argkw[j],
"[rast1.y]") != NULL) ||
1560 (strstr(argkw[j],
"[rast2.x]") != NULL) ||
1561 (strstr(argkw[j],
"[rast2.y]") != NULL)
1563 argtype[k] = INT4OID;
1566 argtype[k] = FLOAT8OID;
1578 elog(ERROR,
"RASTER_nMapAlgebraExpr: Could not create prepared plan of expression parameter %d", exprpos[i]);
1584 POSTGIS_RT_DEBUGF(3,
"expression parameter %d has no args, simply executing", exprpos[i]);
1585 err = SPI_execute(
sql,
TRUE, 0);
1588 if (err != SPI_OK_SELECT || SPI_tuptable == NULL || SPI_processed != 1) {
1591 elog(ERROR,
"RASTER_nMapAlgebraExpr: Could not evaluate expression parameter %d", exprpos[i]);
1596 tupdesc = SPI_tuptable->tupdesc;
1597 tuptable = SPI_tuptable;
1598 tuple = tuptable->vals[0];
1600 datum = SPI_getbinval(tuple, tupdesc, 1, &isnull);
1601 if (SPI_result == SPI_ERROR_NOATTRIBUTE) {
1602 if (SPI_tuptable) SPI_freetuptable(tuptable);
1605 elog(ERROR,
"RASTER_nMapAlgebraExpr: Could not get result of expression parameter %d", exprpos[i]);
1614 if (SPI_tuptable) SPI_freetuptable(tuptable);
1634 for (i = 0; i < numraster; i++) {
1658 if (itrset == NULL) {
1661 elog(ERROR,
"RASTER_nMapAlgebra: Could not allocate memory for iterator arguments");
1666 for (i = 0; i < numraster; i++) {
1690 elog(ERROR,
"RASTER_nMapAlgebraExpr: Could not run raster iterator function");
1693 else if (
raster == NULL) {
1699 MemoryContextSwitchTo(mainMemCtx);
1710 SET_VARSIZE(pgraster, pgraster->
size);
1711 PG_RETURN_POINTER(pgraster);
1731 assert(cutype && strlen(cutype) > 0);
1733 if (strcmp(cutype,
"LAST") == 0)
1735 else if (strcmp(cutype,
"FIRST") == 0)
1737 else if (strcmp(cutype,
"MIN") == 0)
1739 else if (strcmp(cutype,
"MAX") == 0)
1741 else if (strcmp(cutype,
"COUNT") == 0)
1743 else if (strcmp(cutype,
"SUM") == 0)
1745 else if (strcmp(cutype,
"MEAN") == 0)
1747 else if (strcmp(cutype,
"RANGE") == 0)
1774 for (i = 0; i < arg->
numband; i++) {
1798 double *
value,
int *nodata
1810 elog(ERROR,
"rtpg_union_callback: Invalid arguments passed to callback");
1826 else if (!arg->
nodata[0][0][0] && arg->
nodata[1][0][0]) {
1832 else if (arg->
nodata[0][0][0] && !arg->
nodata[1][0][0]) {
1860 else if (!arg->
nodata[0][0][0] && arg->
nodata[1][0][0])
1863 else if (arg->
nodata[0][0][0] && !arg->
nodata[1][0][0])
1889 double *
value,
int *nodata
1899 elog(ERROR,
"rtpg_union_mean_callback: Invalid arguments passed to callback");
1925 double *
value,
int *nodata
1935 elog(ERROR,
"rtpg_union_range_callback: Invalid arguments passed to callback");
1968 HeapTupleHeader tup;
1974 char *utypename = NULL;
1977 etype = ARR_ELEMTYPE(array);
1978 get_typlenbyvalalign(etype, &typlen, &typbyval, &typalign);
1983 typlen, typbyval, typalign,
1988 elog(ERROR,
"rtpg_union_unionarg_process: Invalid argument for unionarg");
1996 elog(ERROR,
"rtpg_union_unionarg_process: Could not allocate memory for band information");
2001 for (i = 0; i < n; i++) {
2010 tup = (HeapTupleHeader) DatumGetPointer(e[i]);
2012 elog(ERROR,
"rtpg_union_unionarg_process: Invalid argument for unionarg");
2017 tupv = GetAttributeByName(tup,
"nband", &isnull);
2020 elog(NOTICE,
"First argument (nband) of unionarg is NULL. Assuming nband = %d",
nband);
2023 nband = DatumGetInt32(tupv);
2026 elog(ERROR,
"rtpg_union_unionarg_process: Band number must be greater than zero (1-based)");
2031 tupv = GetAttributeByName(tup,
"uniontype", &isnull);
2033 elog(NOTICE,
"Second argument (uniontype) of unionarg is NULL. Assuming uniontype = LAST");
2058 elog(ERROR,
"rtpg_union_unionarg_process: Could not reallocate memory for band information");
2075 if (numbands <= arg->numband)
2085 elog(ERROR,
"rtpg_union_noarg: Could not reallocate memory for band information");
2091 for (; i < arg->
numband; i++) {
2099 elog(ERROR,
"rtpg_union_noarg: Could not allocate memory for working rasters");
2108 elog(ERROR,
"rtpg_union_noarg: Could not create working raster");
2121 MemoryContext aggcontext;
2122 MemoryContext oldcontext;
2132 int isempty[2] = {0};
2133 int hasband[2] = {0};
2135 double _offset[4] = {0.};
2143 char *utypename = NULL;
2147 double nodataval = 0;
2153 uint16_t _dim[2] = {0};
2160 if (!AggCheckCallContext(fcinfo, &aggcontext)) {
2161 elog(ERROR,
"RASTER_union_transfn: Cannot be called in a non-aggregate context");
2166 oldcontext = MemoryContextSwitchTo(aggcontext);
2168 if (PG_ARGISNULL(0)) {
2173 MemoryContextSwitchTo(oldcontext);
2174 elog(ERROR,
"RASTER_union_transfn: Could not allocate memory for state variable");
2190 if (!PG_ARGISNULL(1)) {
2192 pgraster = (
rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
2199 PG_FREE_IF_COPY(pgraster, 1);
2201 MemoryContextSwitchTo(oldcontext);
2202 elog(ERROR,
"RASTER_union_transfn: Could not deserialize raster");
2215 if (!PG_ARGISNULL(2)) {
2216 Oid calltype = get_fn_expr_argtype(fcinfo->flinfo, 2);
2251 if (numband > idx) {
2269 PG_FREE_IF_COPY(pgraster, 1);
2272 MemoryContextSwitchTo(oldcontext);
2273 elog(ERROR,
"RASTER_union_transfn: Could not allocate memory for band information");
2278 for (i = idx; i < iwr->
numband; i++) {
2302 nband = PG_GETARG_INT32(2);
2308 PG_FREE_IF_COPY(pgraster, 1);
2311 MemoryContextSwitchTo(oldcontext);
2312 elog(ERROR,
"RASTER_union_transfn: Band number must be greater than zero (1-based)");
2323 PG_FREE_IF_COPY(pgraster, 1);
2326 MemoryContextSwitchTo(oldcontext);
2327 elog(ERROR,
"RASTER_union_transfn: Could not allocate memory for band information");
2348 PG_FREE_IF_COPY(pgraster, 1);
2351 MemoryContextSwitchTo(oldcontext);
2352 elog(ERROR,
"RASTER_union_transfn: Could not process unionarg");
2361 if (nargs > 3 && !PG_ARGISNULL(3)) {
2376 for (i = 0; i < iwr->
numband; i++) {
2391 PG_FREE_IF_COPY(pgraster, 1);
2394 MemoryContextSwitchTo(oldcontext);
2395 elog(ERROR,
"RASTER_union_transfn: Could not allocate memory for working raster(s)");
2410 PG_FREE_IF_COPY(pgraster, 1);
2413 MemoryContextSwitchTo(oldcontext);
2414 elog(ERROR,
"RASTER_union_transfn: Could not create working raster");
2431 PG_FREE_IF_COPY(pgraster, 1);
2434 MemoryContextSwitchTo(oldcontext);
2435 elog(ERROR,
"RASTER_union_transfn: Could not check and balance number of bands");
2442 if (itrset == NULL) {
2447 PG_FREE_IF_COPY(pgraster, 1);
2450 MemoryContextSwitchTo(oldcontext);
2451 elog(ERROR,
"RASTER_union_transfn: Could not allocate memory for iterator arguments");
2456 for (i = 0; i < iwr->
numband; i++) {
2476 if (!isempty[0] && hasband[0])
2478 else if (!isempty[1] && hasband[1])
2485 if (_band != NULL) {
2523 itrset[0].
nband = 0;
2539 if (!isempty[0] && !isempty[1] && hasband[0] && hasband[1]) {
2553 PG_FREE_IF_COPY(pgraster, 1);
2556 MemoryContextSwitchTo(oldcontext);
2557 elog(ERROR,
"RASTER_union_transfn: Could not create internal raster");
2561 _offset[0], _offset[1], _offset[2], _offset[3]);
2568 double igt[6] = {0};
2584 hasnodata, nodataval,
2593 PG_FREE_IF_COPY(pgraster, 1);
2596 MemoryContextSwitchTo(oldcontext);
2597 elog(ERROR,
"RASTER_union_transfn: Could not add new band to internal raster");
2605 for (
y = 0;
y < _dim[1];
y++) {
2606 POSTGIS_RT_DEBUGF(4,
"Getting pixel line of working raster at (x, y, length) = (0, %d, %d)",
y, _dim[0]);
2620 PG_FREE_IF_COPY(pgraster, 1);
2623 MemoryContextSwitchTo(oldcontext);
2624 elog(ERROR,
"RASTER_union_transfn: Could not get pixel line from band of working raster");
2628 POSTGIS_RT_DEBUGF(4,
"Setting pixel line at (x, y, length) = (%d, %d, %d)", (
int) _offset[0], (
int) _offset[1] +
y, nvals);
2631 (
int) _offset[0], (
int) _offset[1] +
y,
2641 PG_FREE_IF_COPY(pgraster, 1);
2644 MemoryContextSwitchTo(oldcontext);
2645 elog(ERROR,
"RASTER_union_transfn: Could not set pixel line to band of internal raster");
2661 hasnodata, nodataval,
2678 PG_FREE_IF_COPY(pgraster, 1);
2681 MemoryContextSwitchTo(oldcontext);
2682 elog(ERROR,
"RASTER_union_transfn: Could not run raster iterator function");
2690 for (
y = 0;
y < _dim[1];
y++) {
2691 POSTGIS_RT_DEBUGF(4,
"Getting pixel line of iterator raster at (x, y, length) = (0, %d, %d)",
y, _dim[0]);
2707 PG_FREE_IF_COPY(pgraster, 1);
2710 MemoryContextSwitchTo(oldcontext);
2711 elog(ERROR,
"RASTER_union_transfn: Could not get pixel line from band of working raster");
2715 POSTGIS_RT_DEBUGF(4,
"Setting pixel line at (x, y, length) = (%d, %d, %d)", (
int) _offset[2], (
int) _offset[3] +
y, nvals);
2718 (
int) _offset[2], (
int) _offset[3] +
y,
2730 PG_FREE_IF_COPY(pgraster, 1);
2733 MemoryContextSwitchTo(oldcontext);
2734 elog(ERROR,
"RASTER_union_transfn: Could not set pixel line to band of internal raster");
2754 hasnodata, nodataval,
2768 PG_FREE_IF_COPY(pgraster, 1);
2771 MemoryContextSwitchTo(oldcontext);
2772 elog(ERROR,
"RASTER_union_transfn: Could not run raster iterator function");
2791 PG_FREE_IF_COPY(pgraster, 1);
2795 MemoryContextSwitchTo(oldcontext);
2799 PG_RETURN_POINTER(iwr);
2819 double nodataval = 0;
2824 if (!AggCheckCallContext(fcinfo, NULL)) {
2825 elog(ERROR,
"RASTER_union_finalfn: Cannot be called in a non-aggregate context");
2830 if (PG_ARGISNULL(0))
2837 if (itrset == NULL) {
2839 elog(ERROR,
"RASTER_union_finalfn: Could not allocate memory for iterator arguments");
2843 for (i = 0; i < iwr->
numband; i++) {
2858 itrset[0].
nband = 0;
2860 itrset[1].
nband = 0;
2868 hasnodata, nodataval,
2881 hasnodata, nodataval,
2895 elog(ERROR,
"RASTER_union_finalfn: Could not run raster iterator function");
2901 if (_raster == NULL)
2909 status = (_rtn == NULL) ? -1 : 0;
2934 elog(ERROR,
"RASTER_union_finalfn: Could not add band to final raster");
2943 if (!_rtn) PG_RETURN_NULL();
2953 SET_VARSIZE(pgraster, pgraster->
size);
2954 PG_RETURN_POINTER(pgraster);
2983 elog(ERROR,
"rtpg_clip_arg_init: Could not allocate memory for function arguments");
2997 if (arg->
band != NULL)
3002 if (arg->
mask != NULL)
3010 double *
value,
int *nodata
3038 unsigned char *wkb = NULL;
3067 if (PG_ARGISNULL(0) || PG_ARGISNULL(2))
3073 elog(ERROR,
"RASTER_clip: Could not initialize argument structure");
3078 pgraster = (
rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
3082 if (arg->
raster == NULL) {
3084 PG_FREE_IF_COPY(pgraster, 0);
3085 elog(ERROR,
"RASTER_clip: Could not deserialize raster");
3091 elog(NOTICE,
"Input raster is empty or has no bands. Returning empty raster");
3094 PG_FREE_IF_COPY(pgraster, 0);
3098 elog(ERROR,
"RASTER_clip: Could not create empty raster");
3107 SET_VARSIZE(pgrtn, pgrtn->
size);
3108 PG_RETURN_POINTER(pgrtn);
3116 gser = PG_GETARG_GSERIALIZED_P(2);
3128 elog(NOTICE,
"Geometry provided does not have the same SRID as the raster. Returning NULL");
3131 PG_FREE_IF_COPY(pgraster, 0);
3133 PG_FREE_IF_COPY(gser, 2);
3139 if (!PG_ARGISNULL(4) && !PG_GETARG_BOOL(4))
3146 PG_FREE_IF_COPY(pgraster, 0);
3148 PG_FREE_IF_COPY(gser, 2);
3150 elog(ERROR,
"RASTER_clip: Could not get convex hull of raster");
3157 PG_FREE_IF_COPY(gser, 2);
3162 elog(NOTICE,
"The input raster and input geometry do not intersect. Returning empty raster");
3165 PG_FREE_IF_COPY(pgraster, 0);
3170 elog(ERROR,
"RASTER_clip: Could not create empty raster");
3179 SET_VARSIZE(pgrtn, pgrtn->
size);
3180 PG_RETURN_POINTER(pgrtn);
3184 if (!PG_ARGISNULL(1)) {
3185 array = PG_GETARG_ARRAYTYPE_P(1);
3186 etype = ARR_ELEMTYPE(array);
3187 get_typlenbyvalalign(etype, &typlen, &typbyval, &typalign);
3195 PG_FREE_IF_COPY(pgraster, 0);
3197 elog(ERROR,
"RASTER_clip: Invalid data type for band indexes");
3204 typlen, typbyval, typalign,
3209 if (arg->
band == NULL) {
3211 PG_FREE_IF_COPY(pgraster, 0);
3213 elog(ERROR,
"RASTER_clip: Could not allocate memory for band arguments");
3217 for (i = 0, j = 0; i < arg->
numbands; i++) {
3218 if (nulls[i])
continue;
3222 arg->
band[j].
nband = DatumGetInt16(e[i]) - 1;
3225 arg->
band[j].
nband = DatumGetInt32(e[i]) - 1;
3232 if (j < arg->numbands) {
3234 if (arg->
band == NULL) {
3236 PG_FREE_IF_COPY(pgraster, 0);
3238 elog(ERROR,
"RASTER_clip: Could not reallocate memory for band arguments");
3246 for (i = 0; i < arg->
numbands; i++) {
3248 elog(NOTICE,
"Band at index %d not found in raster", arg->
band[i].
nband + 1);
3250 PG_FREE_IF_COPY(pgraster, 0);
3265 if (arg->
band == NULL) {
3268 PG_FREE_IF_COPY(pgraster, 0);
3271 elog(ERROR,
"RASTER_clip: Could not allocate memory for band arguments");
3275 for (i = 0; i < arg->
numbands; i++) {
3284 if (!PG_ARGISNULL(3)) {
3285 array = PG_GETARG_ARRAYTYPE_P(3);
3286 etype = ARR_ELEMTYPE(array);
3287 get_typlenbyvalalign(etype, &typlen, &typbyval, &typalign);
3295 PG_FREE_IF_COPY(pgraster, 0);
3297 elog(ERROR,
"RASTER_clip: Invalid data type for NODATA values");
3304 typlen, typbyval, typalign,
3309 for (i = 0, j = 0; i < arg->
numbands; i++, j++) {
3350 if (arg->
mask == NULL) {
3352 PG_FREE_IF_COPY(pgraster, 0);
3353 elog(ERROR,
"RASTER_clip: Could not rasterize intersection geometry");
3364 if (itrset == NULL) {
3366 PG_FREE_IF_COPY(pgraster, 0);
3367 elog(ERROR,
"RASTER_clip: Could not allocate memory for iterator arguments");
3372 for (i = 0; i < arg->
numbands; i++) {
3373 POSTGIS_RT_DEBUGF(4,
"band arg %d (nband, hasnodata, nodataval) = (%d, %d, %f)",
3401 PG_FREE_IF_COPY(pgraster, 0);
3402 elog(ERROR,
"RASTER_clip: Could not create output raster");
3411 PG_FREE_IF_COPY(pgraster, 0);
3412 elog(ERROR,
"RASTER_clip: Could not add NODATA band to output raster");
3426 itrset[1].
nband = 0;
3434 hasnodata, nodataval,
3446 PG_FREE_IF_COPY(pgraster, 0);
3447 elog(ERROR,
"RASTER_clip: Could not run raster iterator function");
3462 PG_FREE_IF_COPY(pgraster, 0);
3463 elog(NOTICE,
"RASTER_clip: Could not get band from working raster");
3472 PG_FREE_IF_COPY(pgraster, 0);
3473 elog(ERROR,
"RASTER_clip: Could not add new band to output raster");
3483 PG_FREE_IF_COPY(pgraster, 0);
3493 SET_VARSIZE(pgrtn, pgrtn->
size);
3494 PG_RETURN_POINTER(pgrtn);
3527 HeapTupleHeader tup;
3532 text *exprtext = NULL;
3537 char *pixeltype = NULL;
3538 text *pixeltypetext = NULL;
3540 double nodataval = 0;
3541 bool hasnodata =
FALSE;
3543 char **comma_set = NULL;
3545 char **colon_set = NULL;
3547 char **dash_set = NULL;
3553 if (PG_ARGISNULL(0))
3555 pgraster = (
rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
3560 PG_FREE_IF_COPY(pgraster, 0);
3561 elog(ERROR,
"RASTER_reclass: Could not deserialize raster");
3565 POSTGIS_RT_DEBUGF(3,
"RASTER_reclass: %d possible bands to be reclassified", numBands);
3569 array = PG_GETARG_ARRAYTYPE_P(1);
3570 etype = ARR_ELEMTYPE(array);
3571 get_typlenbyvalalign(etype, &typlen, &typbyval, &typalign);
3573 deconstruct_array(array, etype, typlen, typbyval, typalign, &e,
3577 elog(NOTICE,
"Invalid argument for reclassargset. Returning original raster");
3581 PG_FREE_IF_COPY(pgraster, 0);
3585 SET_VARSIZE(pgrtn, pgrtn->
size);
3586 PG_RETURN_POINTER(pgrtn);
3594 for (i = 0; i < n; i++) {
3595 if (nulls[i])
continue;
3598 tup = (HeapTupleHeader) DatumGetPointer(e[i]);
3600 elog(NOTICE,
"Invalid argument for reclassargset. Returning original raster");
3604 PG_FREE_IF_COPY(pgraster, 0);
3608 SET_VARSIZE(pgrtn, pgrtn->
size);
3609 PG_RETURN_POINTER(pgrtn);
3613 tupv = GetAttributeByName(tup,
"nband", &isnull);
3615 elog(NOTICE,
"Invalid argument for reclassargset. Missing value of nband for reclassarg of index %d . Returning original raster", i);
3619 PG_FREE_IF_COPY(pgraster, 0);
3623 SET_VARSIZE(pgrtn, pgrtn->
size);
3624 PG_RETURN_POINTER(pgrtn);
3626 nband = DatumGetInt32(tupv);
3630 if (nband < 1 || nband > numBands) {
3631 elog(NOTICE,
"Invalid argument for reclassargset. Invalid band index (must use 1-based) for reclassarg of index %d . Returning original raster", i);
3635 PG_FREE_IF_COPY(pgraster, 0);
3639 SET_VARSIZE(pgrtn, pgrtn->
size);
3640 PG_RETURN_POINTER(pgrtn);
3644 tupv = GetAttributeByName(tup,
"reclassexpr", &isnull);
3646 elog(NOTICE,
"Invalid argument for reclassargset. Missing value of reclassexpr for reclassarg of index %d . Returning original raster", i);
3650 PG_FREE_IF_COPY(pgraster, 0);
3654 SET_VARSIZE(pgrtn, pgrtn->
size);
3655 PG_RETURN_POINTER(pgrtn);
3657 exprtext = (text *) DatumGetPointer(tupv);
3658 if (NULL == exprtext) {
3659 elog(NOTICE,
"Invalid argument for reclassargset. Missing value of reclassexpr for reclassarg of index %d . Returning original raster", i);
3663 PG_FREE_IF_COPY(pgraster, 0);
3667 SET_VARSIZE(pgrtn, pgrtn->
size);
3668 PG_RETURN_POINTER(pgrtn);
3679 elog(NOTICE,
"Invalid argument for reclassargset. Invalid expression of reclassexpr for reclassarg of index %d . Returning original raster", i);
3683 PG_FREE_IF_COPY(pgraster, 0);
3687 SET_VARSIZE(pgrtn, pgrtn->
size);
3688 PG_RETURN_POINTER(pgrtn);
3695 for (a = 0, j = 0; a < comma_n; a++) {
3701 elog(NOTICE,
"Invalid argument for reclassargset. Invalid expression of reclassexpr for reclassarg of index %d . Returning original raster", i);
3702 for (k = 0; k < j; k++) pfree(exprset[k]);
3707 PG_FREE_IF_COPY(pgraster, 0);
3711 SET_VARSIZE(pgrtn, pgrtn->
size);
3712 PG_RETURN_POINTER(pgrtn);
3718 for (b = 0; b < colon_n; b++) {
3723 if (dash_n < 1 || dash_n > 3) {
3724 elog(NOTICE,
"Invalid argument for reclassargset. Invalid expression of reclassexpr for reclassarg of index %d . Returning original raster", i);
3725 for (k = 0; k < j; k++) pfree(exprset[k]);
3730 PG_FREE_IF_COPY(pgraster, 0);
3734 SET_VARSIZE(pgrtn, pgrtn->
size);
3735 PG_RETURN_POINTER(pgrtn);
3738 for (c = 0; c < dash_n; c++) {
3742 strlen(dash_set[c]) == 1 && (
3743 strchr(dash_set[c],
'(') != NULL ||
3744 strchr(dash_set[c],
'[') != NULL ||
3745 strchr(dash_set[c],
')') != NULL ||
3746 strchr(dash_set[c],
']') != NULL
3750 junk = palloc(
sizeof(
char) * (strlen(dash_set[c + 1]) + 2));
3752 for (k = 0; k <= j; k++) pfree(exprset[k]);
3755 PG_FREE_IF_COPY(pgraster, 0);
3757 elog(ERROR,
"RASTER_reclass: Could not allocate memory");
3761 sprintf(junk,
"%s%s", dash_set[c], dash_set[c + 1]);
3763 dash_set[c] = repalloc(dash_set[c],
sizeof(
char) * (strlen(junk) + 1));
3764 strcpy(dash_set[c], junk);
3768 for (dash_it = 1; dash_it < dash_n; dash_it++) {
3769 dash_set[dash_it - 1] = repalloc(dash_set[dash_it - 1], (strlen(dash_set[dash_it]) + 1) *
sizeof(
char));
3770 strcpy(dash_set[dash_it - 1], dash_set[dash_it]);
3774 pfree(dash_set[dash_n]);
3775 dash_set = repalloc(dash_set,
sizeof(
char *) * dash_n);
3779 if (c < 1 && dash_n > 2) {
3780 elog(NOTICE,
"Invalid argument for reclassargset. Invalid expression of reclassexpr for reclassarg of index %d . Returning original raster", i);
3781 for (k = 0; k < j; k++) pfree(exprset[k]);
3786 PG_FREE_IF_COPY(pgraster, 0);
3790 SET_VARSIZE(pgrtn, pgrtn->
size);
3791 PG_RETURN_POINTER(pgrtn);
3802 strchr(dash_set[c],
')') != NULL ||
3803 strchr(dash_set[c],
']') != NULL
3808 else if (strchr(dash_set[c],
'(') != NULL){
3818 strrchr(dash_set[c],
'(') != NULL ||
3819 strrchr(dash_set[c],
'[') != NULL
3824 else if (strrchr(dash_set[c],
']') != NULL) {
3832 POSTGIS_RT_DEBUGF(4,
"RASTER_reclass: exc_val %d inc_val %d", exc_val, inc_val);
3840 val = strtod(dash_set[c], &junk);
3841 if (errno != 0 || dash_set[c] == junk) {
3842 elog(NOTICE,
"Invalid argument for reclassargset. Invalid expression of reclassexpr for reclassarg of index %d . Returning original raster", i);
3843 for (k = 0; k < j; k++) pfree(exprset[k]);
3848 PG_FREE_IF_COPY(pgraster, 0);
3852 SET_VARSIZE(pgrtn, pgrtn->
size);
3853 PG_RETURN_POINTER(pgrtn);
3859 junk = strstr(colon_set[b], dash_set[c]);
3864 "(colon_set[%d], dash_set[%d], junk) = (%s, %s, %s)",
3865 b, c, colon_set[b], dash_set[c], junk
3868 if (junk != colon_set[b]) {
3870 if (*(junk - 1) ==
'-') {
3873 ((junk - 1) == colon_set[b]) ||
3874 (*(junk - 2) ==
'-') ||
3875 (*(junk - 2) ==
'[') ||
3876 (*(junk - 2) ==
'(')
3896 exprset[j]->
src.
min = val;
3902 exprset[j]->
src.
max = val;
3912 exprset[j]->
dst.
min = val;
3915 exprset[j]->
dst.
max = val;
3923 , exprset[j]->src.min
3933 tupv = GetAttributeByName(tup,
"pixeltype", &isnull);
3935 elog(NOTICE,
"Invalid argument for reclassargset. Missing value of pixeltype for reclassarg of index %d . Returning original raster", i);
3939 PG_FREE_IF_COPY(pgraster, 0);
3943 SET_VARSIZE(pgrtn, pgrtn->
size);
3944 PG_RETURN_POINTER(pgrtn);
3946 pixeltypetext = (text *) DatumGetPointer(tupv);
3947 if (NULL == pixeltypetext) {
3948 elog(NOTICE,
"Invalid argument for reclassargset. Missing value of pixeltype for reclassarg of index %d . Returning original raster", i);
3952 PG_FREE_IF_COPY(pgraster, 0);
3956 SET_VARSIZE(pgrtn, pgrtn->
size);
3957 PG_RETURN_POINTER(pgrtn);
3964 tupv = GetAttributeByName(tup,
"nodataval", &isnull);
3970 nodataval = DatumGetFloat8(tupv);
3979 elog(NOTICE,
"Could not find raster band of index %d. Returning original raster",
nband);
3980 for (k = 0; k < j; k++) pfree(exprset[k]);
3985 PG_FREE_IF_COPY(pgraster, 0);
3989 SET_VARSIZE(pgrtn, pgrtn->
size);
3990 PG_RETURN_POINTER(pgrtn);
3994 for (k = 0; k < j; k++) pfree(exprset[k]);
3998 PG_FREE_IF_COPY(pgraster, 0);
4000 elog(ERROR,
"RASTER_reclass: Could not reclassify raster band of index %d",
nband);
4006 for (k = 0; k < j; k++) pfree(exprset[k]);
4011 PG_FREE_IF_COPY(pgraster, 0);
4013 elog(ERROR,
"RASTER_reclass: Could not replace raster band of index %d with reclassified band",
nband);
4021 for (k = 0; k < j; k++) pfree(exprset[k]);
4027 PG_FREE_IF_COPY(pgraster, 0);
4033 SET_VARSIZE(pgrtn, pgrtn->
size);
4034 PG_RETURN_POINTER(pgrtn);
4063 elog(ERROR,
"rtpg_colormap_arg: Could not allocate memory for function arguments");
4074 elog(ERROR,
"rtpg_colormap_arg: Could not allocate memory for function arguments");
4107 for (i = 0; i < arg->
nentry; i++) {
4108 if (arg->
entry[i] != NULL)
4109 pfree(arg->
entry[i]);
4115 for (i = 0; i < arg->
nelement; i++)
4135 if (PG_ARGISNULL(0))
4141 elog(ERROR,
"RASTER_colorMap: Could not initialize argument structure");
4146 pgraster = (
rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
4152 PG_FREE_IF_COPY(pgraster, 0);
4153 elog(ERROR,
"RASTER_colorMap: Could not deserialize raster");
4158 if (!PG_ARGISNULL(1))
4159 arg->
nband = PG_GETARG_INT32(1);
4164 elog(NOTICE,
"Raster does not have band at index %d. Returning empty raster", arg->
nband);
4169 PG_FREE_IF_COPY(pgraster, 0);
4170 elog(ERROR,
"RASTER_colorMap: Could not create empty raster");
4175 PG_FREE_IF_COPY(pgraster, 0);
4179 if (pgraster == NULL)
4182 SET_VARSIZE(pgraster, ((
rt_pgraster*) pgraster)->size);
4183 PG_RETURN_POINTER(pgraster);
4188 if (arg->
band == NULL) {
4191 PG_FREE_IF_COPY(pgraster, 0);
4192 elog(ERROR,
"RASTER_colorMap: Could not get band at index %d",
nband);
4197 if (!PG_ARGISNULL(3)) {
4198 char *method = NULL;
4206 if (strcmp(method,
"INTERPOLATE") == 0)
4208 else if (strcmp(method,
"EXACT") == 0)
4210 else if (strcmp(method,
"NEAREST") == 0)
4213 elog(NOTICE,
"Unknown value provided for method. Defaulting to INTERPOLATE");
4223 if (PG_ARGISNULL(2)) {
4225 PG_FREE_IF_COPY(pgraster, 0);
4226 elog(ERROR,
"RASTER_colorMap: Value must be provided for colormap");
4240 if (!strlen(colormap)) {
4242 PG_FREE_IF_COPY(pgraster, 0);
4243 elog(ERROR,
"RASTER_colorMap: Value must be provided for colormap");
4251 PG_FREE_IF_COPY(pgraster, 0);
4252 elog(ERROR,
"RASTER_colorMap: Could not process the value provided for colormap");
4260 PG_FREE_IF_COPY(pgraster, 0);
4261 elog(ERROR,
"RASTER_colorMap: Could not allocate memory for colormap entries");
4267 for (i = 0; i < arg->
nentry; i++) {
4281 if (!strlen(_entry)) {
4291 PG_FREE_IF_COPY(pgraster, 0);
4292 elog(ERROR,
"RASTER_colorMap: Could not process colormap entry %d", i + 1);
4296 elog(NOTICE,
"More than five elements in colormap entry %d. Using at most five elements", i + 1);
4305 for (j = 0; j < arg->
nelement; j++) {
4314 char *percent = NULL;
4318 strcmp(_element,
"NV") == 0 ||
4319 strcmp(_element,
"NULL") == 0 ||
4320 strcmp(_element,
"NODATA") == 0
4325 elog(NOTICE,
"More than one NODATA entry found. Using only the first one");
4334 else if ((percent = strchr(_element,
'%')) != NULL) {
4346 PG_FREE_IF_COPY(pgraster, 0);
4347 elog(ERROR,
"RASTER_colorMap: Could not get band's summary stats to process percentages");
4353 tmp = palloc(
sizeof(
char) * (percent - _element + 1));
4357 PG_FREE_IF_COPY(pgraster, 0);
4358 elog(ERROR,
"RASTER_colorMap: Could not allocate memory for value of percentage");
4362 memcpy(tmp, _element, percent - _element);
4363 tmp[percent - _element] =
'\0';
4368 value = strtod(tmp, NULL);
4370 if (errno != 0 || _element == junk) {
4373 PG_FREE_IF_COPY(pgraster, 0);
4374 elog(ERROR,
"RASTER_colorMap: Could not process percent string to value");
4380 elog(NOTICE,
"Percentage values cannot be less than zero. Defaulting to zero");
4383 else if (
value > 100.) {
4384 elog(NOTICE,
"Percentage values cannot be greater than 100. Defaulting to 100");
4396 if (errno != 0 || _element == junk) {
4399 PG_FREE_IF_COPY(pgraster, 0);
4400 elog(ERROR,
"RASTER_colorMap: Could not process string to value");
4411 value = (int) strtod(_element, &junk);
4412 if (errno != 0 || _element == junk) {
4415 PG_FREE_IF_COPY(pgraster, 0);
4416 elog(ERROR,
"RASTER_colorMap: Could not process string to value");
4421 elog(NOTICE,
"RGBA value cannot be greater than 255. Defaulting to 255");
4424 else if (
value < 0) {
4425 elog(NOTICE,
"RGBA value cannot be less than zero. Defaulting to zero");
4434 POSTGIS_RT_DEBUGF(4,
"colormap->entry[%d] (isnodata, value, R, G, B, A) = (%d, %f, %d, %d, %d, %d)",
4455 PG_FREE_IF_COPY(pgraster, 0);
4456 elog(ERROR,
"RASTER_colorMap: Could not create new raster with applied colormap");
4461 PG_FREE_IF_COPY(pgraster, 0);
4467 if (pgraster == NULL)
4470 SET_VARSIZE(pgraster, ((
rt_pgraster*) pgraster)->size);
4471 PG_RETURN_POINTER(pgraster);
4483 int x,
y,
nband, width, height;
4485 double newnodatavalue = 0.0;
4486 double newinitialvalue = 0.0;
4487 double newval = 0.0;
4488 char *newexpr = NULL;
4489 char *initexpr = NULL;
4490 char *expression = NULL;
4491 int hasnodataval = 0;
4492 double nodataval = 0.;
4494 int skipcomputation = 0;
4496 const int argkwcount = 3;
4497 enum KEYWORDS { kVAL=0, kX=1, kY=2 };
4498 char *argkw[] = {
"[rast]",
"[rast.x]",
"[rast.y]"};
4499 Oid argkwtypes[] = { FLOAT8OID, INT4OID, INT4OID };
4501 Oid argtype[] = { FLOAT8OID, INT4OID, INT4OID };
4507 SPIPlanPtr spi_plan = NULL;
4508 SPITupleTable * tuptable = NULL;
4510 char * strFromText = NULL;
4511 Datum *values = NULL;
4512 Datum datum = (Datum)NULL;
4514 bool isnull =
FALSE;
4521 if (PG_ARGISNULL(0)) {
4522 elog(NOTICE,
"Raster is NULL. Returning NULL");
4528 pgraster = (
rt_pgraster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
4531 PG_FREE_IF_COPY(pgraster, 0);
4532 elog(ERROR,
"RASTER_mapAlgebraExpr: Could not deserialize raster");
4538 if (PG_ARGISNULL(1))
4541 nband = PG_GETARG_INT32(1);
4547 POSTGIS_RT_DEBUG(3,
"RASTER_mapAlgebraExpr: Creating new empty raster...");
4558 if ( NULL == newrast ) {
4559 PG_FREE_IF_COPY(pgraster, 0);
4560 elog(ERROR,
"RASTER_mapAlgebraExpr: Could not create a new raster");
4585 elog(NOTICE,
"Raster is empty. Returning an empty raster");
4587 PG_FREE_IF_COPY(pgraster, 0);
4591 if (NULL == pgrtn) {
4593 elog(ERROR,
"RASTER_mapAlgebraExpr: Could not serialize raster");
4597 SET_VARSIZE(pgrtn, pgrtn->
size);
4598 PG_RETURN_POINTER(pgrtn);
4609 elog(NOTICE,
"Raster does not have the required band. Returning a raster "
4612 PG_FREE_IF_COPY(pgraster, 0);
4616 if (NULL == pgrtn) {
4617 elog(ERROR,
"RASTER_mapAlgebraExpr: Could not serialize raster");
4621 SET_VARSIZE(pgrtn, pgrtn->
size);
4622 PG_RETURN_POINTER(pgrtn);
4627 if ( NULL ==
band ) {
4628 elog(NOTICE,
"Could not get the required band. Returning a raster "
4631 PG_FREE_IF_COPY(pgraster, 0);
4635 if (NULL == pgrtn) {
4636 elog(ERROR,
"RASTER_mapAlgebraExpr: Could not serialize raster");
4640 SET_VARSIZE(pgrtn, pgrtn->
size);
4641 PG_RETURN_POINTER(pgrtn);
4647 POSTGIS_RT_DEBUG(3,
"RASTER_mapAlgebraExpr: Getting NODATA value for band...");
4665 newinitialvalue = newnodatavalue;
4672 if (PG_ARGISNULL(2)) {
4680 if (newpixeltype ==
PT_END)
4684 if (newpixeltype ==
PT_END) {
4685 PG_FREE_IF_COPY(pgraster, 0);
4686 elog(ERROR,
"RASTER_mapAlgebraExpr: Invalid pixeltype");
4695 if (!PG_ARGISNULL(3)) {
4697 len = strlen(
"SELECT (") + strlen(expression) + strlen(
")::double precision");
4698 initexpr = (
char *)palloc(len + 1);
4700 strncpy(initexpr,
"SELECT (", strlen(
"SELECT ("));
4701 strncpy(initexpr + strlen(
"SELECT ("), expression, strlen(expression));
4702 strncpy(initexpr + strlen(
"SELECT (") + strlen(expression),
")::double precision", strlen(
")::double precision"));
4703 initexpr[len] =
'\0';
4721 if (!PG_ARGISNULL(4)) {
4723 nodataval = PG_GETARG_FLOAT8(4);
4724 newinitialvalue = nodataval;
4741 POSTGIS_RT_DEBUG(3,
"RASTER_mapAlgebraExpr: Band is a nodata band, returning "
4742 "a raster filled with nodata");
4745 newinitialvalue,
TRUE, newnodatavalue, 0);
4751 PG_FREE_IF_COPY(pgraster, 0);
4756 if (NULL == pgrtn) {
4757 elog(ERROR,
"RASTER_mapAlgebraExpr: Could not serialize raster");
4761 SET_VARSIZE(pgrtn, pgrtn->
size);
4762 PG_RETURN_POINTER(pgrtn);
4770 if (initexpr != NULL && ( !strcmp(initexpr,
"SELECT [rast]") || !strcmp(initexpr,
"SELECT [rast.val]") ) && !hasnodataval) {
4773 "Returning raster with band %d from original raster",
nband);
4785 PG_FREE_IF_COPY(pgraster, 0);
4790 if (NULL == pgrtn) {
4791 elog(ERROR,
"RASTER_mapAlgebraExpr: Could not serialize raster");
4795 SET_VARSIZE(pgrtn, pgrtn->
size);
4796 PG_RETURN_POINTER(pgrtn);
4803 if (initexpr != NULL && strstr(initexpr,
"[rast") == NULL) {
4804 ret = SPI_connect();
4805 if (ret != SPI_OK_CONNECT) {
4806 PG_FREE_IF_COPY(pgraster, 0);
4807 elog(ERROR,
"RASTER_mapAlgebraExpr: Could not connect to the SPI manager");
4812 ret = SPI_execute(initexpr,
FALSE, 0);
4814 if (ret != SPI_OK_SELECT || SPI_tuptable == NULL || SPI_processed != 1) {
4818 SPI_freetuptable(tuptable);
4819 PG_FREE_IF_COPY(pgraster, 0);
4822 elog(ERROR,
"RASTER_mapAlgebraExpr: Invalid construction for expression");
4826 tupdesc = SPI_tuptable->tupdesc;
4827 tuptable = SPI_tuptable;
4829 tuple = tuptable->vals[0];
4830 newexpr = SPI_getvalue(tuple, tupdesc, 1);
4832 POSTGIS_RT_DEBUG(3,
"Constant expression evaluated to NULL, keeping initvalue");
4833 newval = newinitialvalue;
4835 newval = atof(newexpr);
4838 SPI_freetuptable(tuptable);
4845 skipcomputation = 1;
4851 if (!hasnodataval) {
4852 newinitialvalue = newval;
4853 skipcomputation = 2;
4857 else if (
FLT_NEQ(newval, newinitialvalue)) {
4858 skipcomputation = 2;
4867 newinitialvalue,
TRUE, newnodatavalue, 0);
4874 if (expression == NULL || skipcomputation == 2) {
4880 PG_FREE_IF_COPY(pgraster, 0);
4885 if (NULL == pgrtn) {
4886 elog(ERROR,
"RASTER_mapAlgebraExpr: Could not serialize raster");
4890 SET_VARSIZE(pgrtn, pgrtn->
size);
4891 PG_RETURN_POINTER(pgrtn);
4894 RASTER_DEBUG(3,
"RASTER_mapAlgebraExpr: Creating new raster band...");
4898 if ( NULL == newband ) {
4899 elog(NOTICE,
"Could not modify band for new raster. Returning new "
4900 "raster with the original band");
4905 PG_FREE_IF_COPY(pgraster, 0);
4910 if (NULL == pgrtn) {
4911 elog(ERROR,
"RASTER_mapAlgebraExpr: Could not serialize raster");
4915 SET_VARSIZE(pgrtn, pgrtn->
size);
4916 PG_RETURN_POINTER(pgrtn);
4922 if (initexpr != NULL) {
4925 pfree(initexpr); initexpr=newexpr;
4927 sprintf(place,
"$1");
4928 for (i = 0, j = 1; i < argkwcount; i++) {
4931 pfree(initexpr); initexpr=newexpr;
4933 argtype[argcount] = argkwtypes[i];
4937 sprintf(place,
"$%d", j);
4947 values = (Datum *) palloc(
sizeof(Datum) * argcount);
4948 if (values == NULL) {
4953 PG_FREE_IF_COPY(pgraster, 0);
4956 elog(ERROR,
"RASTER_mapAlgebraExpr: Could not allocate memory for value parameters of prepared statement");
4961 nulls = (
char *)palloc(argcount);
4962 if (nulls == NULL) {
4967 PG_FREE_IF_COPY(pgraster, 0);
4970 elog(ERROR,
"RASTER_mapAlgebraExpr: Could not allocate memory for null parameters of prepared statement");
4975 ret = SPI_connect();
4976 if (ret != SPI_OK_CONNECT) {
4981 PG_FREE_IF_COPY(pgraster, 0);
4984 elog(ERROR,
"RASTER_mapAlgebraExpr: Could not connect to the SPI manager");
4989 spi_plan = SPI_prepare(initexpr, argcount, argtype);
4991 if (spi_plan == NULL) {
4994 PG_FREE_IF_COPY(pgraster, 0);
5001 elog(ERROR,
"RASTER_mapAlgebraExpr: Could not prepare expression");
5006 for (
x = 0;
x < width;
x++) {
5007 for(
y = 0;
y < height;
y++) {
5015 if (skipcomputation == 0) {
5016 if (initexpr != NULL) {
5018 memset(nulls,
'n', argcount);
5020 for (i = 0; i < argkwcount; i++) {
5022 if (idx < 1)
continue;
5027 values[idx] = Int32GetDatum(
x+1);
5032 values[idx] = Int32GetDatum(
y+1);
5035 else if (i == kVAL ) {
5036 values[idx] = Float8GetDatum(
r);
5042 ret = SPI_execute_plan(spi_plan, values, nulls,
FALSE, 0);
5043 if (ret != SPI_OK_SELECT || SPI_tuptable == NULL ||
5044 SPI_processed != 1) {
5046 SPI_freetuptable(tuptable);
5048 SPI_freeplan(spi_plan);
5056 PG_FREE_IF_COPY(pgraster, 0);
5059 elog(ERROR,
"RASTER_mapAlgebraExpr: Error executing prepared plan");
5064 tupdesc = SPI_tuptable->tupdesc;
5065 tuptable = SPI_tuptable;
5067 tuple = tuptable->vals[0];
5068 datum = SPI_getbinval(tuple, tupdesc, 1, &isnull);
5069 if ( SPI_result == SPI_ERROR_NOATTRIBUTE ) {
5070 POSTGIS_RT_DEBUGF(3,
"Expression for pixel %d,%d (value %g) errored, skip setting",
x+1,
y+1,
r);
5071 newval = newinitialvalue;
5073 else if ( isnull ) {
5074 POSTGIS_RT_DEBUGF(3,
"Expression for pixel %d,%d (value %g) evaluated to NULL, skip setting",
x+1,
y+1,
r);
5075 newval = newinitialvalue;
5077 newval = DatumGetFloat8(datum);
5080 SPI_freetuptable(tuptable);
5084 newval = newinitialvalue;
5097 if (initexpr != NULL) {
5098 SPI_freeplan(spi_plan);
5112 POSTGIS_RT_DEBUG(3,
"RASTER_mapAlgebraExpr: raster modified, serializing it.");
5116 PG_FREE_IF_COPY(pgraster, 0);
5123 SET_VARSIZE(pgrtn, pgrtn->
size);
5131 PG_RETURN_POINTER(pgrtn);
5146 int x,
y,
nband, width, height;
5148 double newnodatavalue = 0.0;
5149 double newinitialvalue = 0.0;
5150 double newval = 0.0;
5155 #if POSTGIS_PGSQL_VERSION < 120
5156 FunctionCallInfoData cbdata;
5158 LOCAL_FCINFO(cbdata, FUNC_MAX_ARGS);
5161 char * strFromText = NULL;
5167 if (PG_ARGISNULL(0)) {
5168 elog(WARNING,
"Raster is NULL. Returning NULL");
5174 pgraster = (
rt_pgraster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
5177 PG_FREE_IF_COPY(pgraster, 0);
5178 elog(ERROR,
"RASTER_mapAlgebraFct: Could not deserialize raster");
5186 if (PG_ARGISNULL(1))
5189 nband = PG_GETARG_INT32(1);
5205 if ( NULL == newrast ) {
5208 PG_FREE_IF_COPY(pgraster, 0);
5210 elog(ERROR,
"RASTER_mapAlgebraFct: Could not create a new raster");
5235 elog(NOTICE,
"Raster is empty. Returning an empty raster");
5237 PG_FREE_IF_COPY(pgraster, 0);
5241 if (NULL == pgrtn) {
5242 elog(ERROR,
"RASTER_mapAlgebraFct: Could not serialize raster");
5246 SET_VARSIZE(pgrtn, pgrtn->
size);
5247 PG_RETURN_POINTER(pgrtn);
5257 elog(NOTICE,
"Raster does not have the required band. Returning a raster "
5260 PG_FREE_IF_COPY(pgraster, 0);
5264 if (NULL == pgrtn) {
5265 elog(ERROR,
"RASTER_mapAlgebraFct: Could not serialize raster");
5269 SET_VARSIZE(pgrtn, pgrtn->
size);
5270 PG_RETURN_POINTER(pgrtn);
5275 if ( NULL ==
band ) {
5276 elog(NOTICE,
"Could not get the required band. Returning a raster "
5279 PG_FREE_IF_COPY(pgraster, 0);
5283 if (NULL == pgrtn) {
5284 elog(ERROR,
"RASTER_mapAlgebraFct: Could not serialize raster");
5288 SET_VARSIZE(pgrtn, pgrtn->
size);
5289 PG_RETURN_POINTER(pgrtn);
5295 POSTGIS_RT_DEBUG(3,
"RASTER_mapAlgebraFct: Getting NODATA value for band...");
5312 newinitialvalue = newnodatavalue;
5319 if (PG_ARGISNULL(2)) {
5327 if (newpixeltype ==
PT_END)
5331 if (newpixeltype ==
PT_END) {
5334 PG_FREE_IF_COPY(pgraster, 0);
5337 elog(ERROR,
"RASTER_mapAlgebraFct: Invalid pixeltype");
5345 if (PG_ARGISNULL(3)) {
5348 PG_FREE_IF_COPY(pgraster, 0);
5351 elog(ERROR,
"RASTER_mapAlgebraFct: Required function is missing. Returning NULL");
5355 oid = PG_GETARG_OID(3);
5356 if (oid == InvalidOid) {
5359 PG_FREE_IF_COPY(pgraster, 0);
5362 elog(ERROR,
"RASTER_mapAlgebraFct: Got invalid function object id. Returning NULL");
5366 fmgr_info(oid, &cbinfo);
5369 if (cbinfo.fn_retset) {
5372 PG_FREE_IF_COPY(pgraster, 0);
5375 elog(ERROR,
"RASTER_mapAlgebraFct: Function provided must return double precision not resultset");
5379 else if (cbinfo.fn_nargs < 2 || cbinfo.fn_nargs > 3) {
5382 PG_FREE_IF_COPY(pgraster, 0);
5385 elog(ERROR,
"RASTER_mapAlgebraFct: Function does not have two or three input parameters");
5389 if (cbinfo.fn_nargs == 2)
5394 if (func_volatile(oid) ==
'v') {
5395 elog(NOTICE,
"Function provided is VOLATILE. Unless required and for best performance, function should be IMMUTABLE or STABLE");
5399 #if POSTGIS_PGSQL_VERSION < 120
5400 InitFunctionCallInfoData(cbdata, &cbinfo, 2, InvalidOid, NULL, NULL);
5402 cbdata.argnull[0] =
FALSE;
5403 cbdata.argnull[1] =
FALSE;
5404 cbdata.argnull[2] =
FALSE;
5406 InitFunctionCallInfoData(*cbdata, &cbinfo, 2, InvalidOid, NULL, NULL);
5408 cbdata->args[0].isnull =
FALSE;
5409 cbdata->args[1].isnull =
FALSE;
5410 cbdata->args[2].isnull =
FALSE;
5414 if (PG_ARGISNULL(4)) {
5415 if (cbinfo.fn_strict) {
5418 PG_FREE_IF_COPY(pgraster, 0);
5421 elog(ERROR,
"RASTER_mapAlgebraFct: Strict callback functions cannot have null parameters");
5425 #if POSTGIS_PGSQL_VERSION < 120
5426 cbdata.arg[k] = (Datum)NULL;
5427 cbdata.argnull[k] =
TRUE;
5429 cbdata->args[k].value = (Datum)NULL;
5430 cbdata->args[k].isnull =
TRUE;
5434 #if POSTGIS_PGSQL_VERSION < 120
5435 cbdata.arg[k] = PG_GETARG_DATUM(4);
5437 cbdata->args[k].value = PG_GETARG_DATUM(4);
5448 POSTGIS_RT_DEBUG(3,
"RASTER_mapAlgebraFct: Band is a nodata band, returning "
5449 "a raster filled with nodata");
5452 newinitialvalue,
TRUE, newnodatavalue, 0);
5455 PG_FREE_IF_COPY(pgraster, 0);
5460 if (NULL == pgrtn) {
5461 elog(ERROR,
"RASTER_mapAlgebraFct: Could not serialize raster");
5465 SET_VARSIZE(pgrtn, pgrtn->
size);
5466 PG_RETURN_POINTER(pgrtn);
5475 newinitialvalue,
TRUE, newnodatavalue, 0);
5479 if ( NULL == newband ) {
5480 elog(NOTICE,
"Could not modify band for new raster. Returning new "
5481 "raster with the original band");
5484 PG_FREE_IF_COPY(pgraster, 0);
5489 if (NULL == pgrtn) {
5490 elog(ERROR,
"RASTER_mapAlgebraFct: Could not serialize raster");
5494 SET_VARSIZE(pgrtn, pgrtn->
size);
5495 PG_RETURN_POINTER(pgrtn);
5501 for (
x = 0;
x < width;
x++) {
5502 for(
y = 0;
y < height;
y++) {
5510 if (
FLT_EQ(
r, newnodatavalue)) {
5511 if (cbinfo.fn_strict) {
5512 POSTGIS_RT_DEBUG(3,
"RASTER_mapAlgebraFct: Strict callbacks cannot accept NULL arguments, skipping NODATA cell.");
5515 #if POSTGIS_PGSQL_VERSION < 120
5516 cbdata.argnull[0] =
TRUE;
5517 cbdata.arg[0] = (Datum)NULL;
5519 cbdata->args[0].isnull =
TRUE;
5520 cbdata->args[0].value = (Datum)NULL;
5524 #if POSTGIS_PGSQL_VERSION < 120
5525 cbdata.argnull[0] =
FALSE;
5526 cbdata.arg[0] = Float8GetDatum(
r);
5528 cbdata->args[0].isnull =
FALSE;
5529 cbdata->args[0].value = Float8GetDatum(
r);
5534 if (cbinfo.fn_nargs == 3) {
5538 d[0] = Int32GetDatum(
x+1);
5539 d[1] = Int32GetDatum(
y+1);
5541 a = construct_array(d, 2, INT4OID,
sizeof(
int32),
true,
'i');
5543 #if POSTGIS_PGSQL_VERSION < 120
5544 cbdata.argnull[1] =
FALSE;
5545 cbdata.arg[1] = PointerGetDatum(a);
5547 cbdata->args[1].isnull =
FALSE;
5548 cbdata->args[1].value = PointerGetDatum(a);
5555 #if POSTGIS_PGSQL_VERSION < 120
5556 tmpnewval = FunctionCallInvoke(&cbdata);
5558 if (cbdata.isnull) {
5559 newval = newnodatavalue;
5562 tmpnewval = FunctionCallInvoke(cbdata);
5566 newval = newnodatavalue;
5570 newval = DatumGetFloat8(tmpnewval);
5584 POSTGIS_RT_DEBUG(3,
"RASTER_mapAlgebraFct: raster modified, serializing it.");
5588 PG_FREE_IF_COPY(pgraster, 0);
5599 SET_VARSIZE(pgrtn, pgrtn->
size);
5600 PG_RETURN_POINTER(pgrtn);
5615 int x,
y,
nband, width, height, ngbwidth, ngbheight, winwidth, winheight, u, v, nIndex, nNullItems;
5617 double newnodatavalue = 0.0;
5618 double newinitialvalue = 0.0;
5619 double newval = 0.0;
5624 #if POSTGIS_PGSQL_VERSION < 120
5625 FunctionCallInfoData cbdata;
5627 LOCAL_FCINFO(cbdata, FUNC_MAX_ARGS);
5630 ArrayType * neighborDatum;
5631 char * strFromText = NULL;
5632 text * txtNodataMode = NULL;
5633 text * txtCallbackParam = NULL;
5635 float fltReplace = 0;
5636 bool valuereplace =
false, pixelreplace, nNodataOnly =
true, nNullSkip =
false;
5637 Datum * neighborData = NULL;
5638 bool * neighborNulls = NULL;
5639 int neighborDims[2];
5648 if (PG_ARGISNULL(0)) {
5649 elog(WARNING,
"Raster is NULL. Returning NULL");
5655 pgraster = (
rt_pgraster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
5659 PG_FREE_IF_COPY(pgraster, 0);
5660 elog(ERROR,
"RASTER_mapAlgebraFctNgb: Could not deserialize raster");
5668 if (PG_ARGISNULL(1))
5671 nband = PG_GETARG_INT32(1);
5676 POSTGIS_RT_DEBUG(3,
"RASTER_mapAlgebraFctNgb: Creating new empty raster...");
5687 if ( NULL == newrast ) {
5689 PG_FREE_IF_COPY(pgraster, 0);
5690 elog(ERROR,
"RASTER_mapAlgebraFctNgb: Could not create a new raster");
5715 elog(NOTICE,
"Raster is empty. Returning an empty raster");
5717 PG_FREE_IF_COPY(pgraster, 0);
5721 if (NULL == pgrtn) {
5722 elog(ERROR,
"RASTER_mapAlgebraFctNgb: Could not serialize raster");
5726 SET_VARSIZE(pgrtn, pgrtn->
size);
5727 PG_RETURN_POINTER(pgrtn);
5737 elog(NOTICE,
"Raster does not have the required band. Returning a raster "
5740 PG_FREE_IF_COPY(pgraster, 0);
5744 if (NULL == pgrtn) {
5745 elog(ERROR,
"RASTER_mapAlgebraFctNgb: Could not serialize raster");
5749 SET_VARSIZE(pgrtn, pgrtn->
size);
5750 PG_RETURN_POINTER(pgrtn);
5755 if ( NULL ==
band ) {
5756 elog(NOTICE,
"Could not get the required band. Returning a raster "
5759 PG_FREE_IF_COPY(pgraster, 0);
5763 if (NULL == pgrtn) {
5764 elog(ERROR,
"RASTER_mapAlgebraFctNgb: Could not serialize raster");
5768 SET_VARSIZE(pgrtn, pgrtn->
size);
5769 PG_RETURN_POINTER(pgrtn);
5775 POSTGIS_RT_DEBUG(3,
"RASTER_mapAlgebraFctNgb: Getting NODATA value for band...");
5792 newinitialvalue = newnodatavalue;
5799 if (PG_ARGISNULL(2)) {
5805 POSTGIS_RT_DEBUGF(3,
"RASTER_mapAlgebraFctNgb: Pixeltype parameter: %s", strFromText);
5808 if (newpixeltype ==
PT_END)
5812 if (newpixeltype ==
PT_END) {
5815 PG_FREE_IF_COPY(pgraster, 0);
5818 elog(ERROR,
"RASTER_mapAlgebraFctNgb: Invalid pixeltype");
5826 if (PG_ARGISNULL(5)) {
5829 PG_FREE_IF_COPY(pgraster, 0);
5832 elog(ERROR,
"RASTER_mapAlgebraFctNgb: Required function is missing");
5836 oid = PG_GETARG_OID(5);
5837 if (oid == InvalidOid) {
5840 PG_FREE_IF_COPY(pgraster, 0);
5843 elog(ERROR,
"RASTER_mapAlgebraFctNgb: Got invalid function object id");
5847 fmgr_info(oid, &cbinfo);
5850 if (cbinfo.fn_retset) {
5853 PG_FREE_IF_COPY(pgraster, 0);
5856 elog(ERROR,
"RASTER_mapAlgebraFctNgb: Function provided must return double precision not resultset");
5860 else if (cbinfo.fn_nargs != 3) {
5863 PG_FREE_IF_COPY(pgraster, 0);
5866 elog(ERROR,
"RASTER_mapAlgebraFctNgb: Function does not have three input parameters");
5870 if (func_volatile(oid) ==
'v') {
5871 elog(NOTICE,
"Function provided is VOLATILE. Unless required and for best performance, function should be IMMUTABLE or STABLE");
5875 #if POSTGIS_PGSQL_VERSION < 120
5876 InitFunctionCallInfoData(cbdata, &cbinfo, 3, InvalidOid, NULL, NULL);
5877 memset(cbdata.argnull,
FALSE,
sizeof(
bool) * 3);
5879 InitFunctionCallInfoData(*cbdata, &cbinfo, 3, InvalidOid, NULL, NULL);
5880 cbdata->args[0].isnull =
FALSE;
5881 cbdata->args[1].isnull =
FALSE;
5882 cbdata->args[2].isnull =
FALSE;
5886 if (PG_ARGISNULL(7)) {
5887 if (cbinfo.fn_strict) {
5890 PG_FREE_IF_COPY(pgraster, 0);
5893 elog(ERROR,
"RASTER_mapAlgebraFctNgb: Strict callback functions cannot have NULL parameters");
5897 #if POSTGIS_PGSQL_VERSION < 120
5898 cbdata.arg[2] = (Datum)NULL;
5899 cbdata.argnull[2] =
TRUE;
5901 cbdata->args[2].value = (Datum)NULL;
5902 cbdata->args[2].isnull =
TRUE;
5906 #if POSTGIS_PGSQL_VERSION < 120
5907 cbdata.arg[2] = PG_GETARG_DATUM(7);
5909 cbdata->args[2].value = PG_GETARG_DATUM(7);
5920 POSTGIS_RT_DEBUG(3,
"RASTER_mapAlgebraFctNgb: Band is a nodata band, returning "
5921 "a raster filled with nodata");
5924 newinitialvalue,
TRUE, newnodatavalue, 0);
5927 PG_FREE_IF_COPY(pgraster, 0);
5932 if (NULL == pgrtn) {
5933 elog(ERROR,
"RASTER_mapAlgebraFctNgb: Could not serialize raster");
5937 SET_VARSIZE(pgrtn, pgrtn->
size);
5938 PG_RETURN_POINTER(pgrtn);
5947 newinitialvalue,
TRUE, newnodatavalue, 0);
5951 if ( NULL == newband ) {
5952 elog(NOTICE,
"Could not modify band for new raster. Returning new "
5953 "raster with the original band");
5956 PG_FREE_IF_COPY(pgraster, 0);
5961 if (NULL == pgrtn) {
5962 elog(ERROR,
"RASTER_mapAlgebraFctNgb: Could not serialize raster");
5966 SET_VARSIZE(pgrtn, pgrtn->
size);
5967 PG_RETURN_POINTER(pgrtn);
5971 if (PG_ARGISNULL(3) || PG_GETARG_INT32(3) <= 0) {
5972 elog(NOTICE,
"Neighborhood width is NULL or <= 0. Returning new "
5973 "raster with the original band");
5976 PG_FREE_IF_COPY(pgraster, 0);
5981 if (NULL == pgrtn) {
5982 elog(ERROR,
"RASTER_mapAlgebraFctNgb: Could not serialize raster");
5986 SET_VARSIZE(pgrtn, pgrtn->
size);
5987 PG_RETURN_POINTER(pgrtn);
5990 ngbwidth = PG_GETARG_INT32(3);
5991 winwidth = ngbwidth * 2 + 1;
5994 if (PG_ARGISNULL(4) || PG_GETARG_INT32(4) <= 0) {
5995 elog(NOTICE,
"Neighborhood height is NULL or <= 0. Returning new "
5996 "raster with the original band");
5999 PG_FREE_IF_COPY(pgraster, 0);
6004 if (NULL == pgrtn) {
6005 elog(ERROR,
"RASTER_mapAlgebraFctNgb: Could not serialize raster");
6009 SET_VARSIZE(pgrtn, pgrtn->
size);
6010 PG_RETURN_POINTER(pgrtn);
6013 ngbheight = PG_GETARG_INT32(4);
6014 winheight = ngbheight * 2 + 1;
6017 if (PG_ARGISNULL(6)) {
6018 elog(NOTICE,
"Neighborhood NODATA behavior defaulting to 'ignore'");
6019 txtNodataMode = cstring_to_text(
"ignore");
6022 txtNodataMode = PG_GETARG_TEXT_P(6);
6025 txtCallbackParam = (text*)palloc(VARSIZE(txtNodataMode));
6026 SET_VARSIZE(txtCallbackParam, VARSIZE(txtNodataMode));
6027 memcpy((
void *)VARDATA(txtCallbackParam), (
void *)VARDATA(txtNodataMode), VARSIZE(txtNodataMode) - VARHDRSZ);
6030 #if POSTGIS_PGSQL_VERSION < 120
6031 cbdata.arg[1] = CStringGetDatum(txtCallbackParam);
6033 cbdata->args[1].value = CStringGetDatum(txtCallbackParam);
6039 if (strcmp(strFromText,
"VALUE") == 0)
6040 valuereplace =
true;
6041 else if (strcmp(strFromText,
"IGNORE") != 0 && strcmp(strFromText,
"NULL") != 0) {
6043 if (sscanf(strFromText,
"%d", &intReplace) <= 0 && sscanf(strFromText,
"%f", &fltReplace) <= 0) {
6045 elog(NOTICE,
"Neighborhood NODATA mode is not recognized. Must be one of 'value', 'ignore', "
6046 "'NULL', or a numeric value. Returning new raster with the original band");
6049 pfree(txtCallbackParam);
6053 PG_FREE_IF_COPY(pgraster, 0);
6058 if (NULL == pgrtn) {
6059 elog(ERROR,
"RASTER_mapAlgebraFctNgb: Could not serialize raster");
6063 SET_VARSIZE(pgrtn, pgrtn->
size);
6064 PG_RETURN_POINTER(pgrtn);
6067 else if (strcmp(strFromText,
"NULL") == 0) {
6076 neighborData = (Datum *)palloc(winwidth * winheight *
sizeof(Datum));
6077 neighborNulls = (
bool *)palloc(winwidth * winheight *
sizeof(
bool));
6080 neighborDims[0] = winwidth;
6081 neighborDims[1] = winheight;
6088 get_typlenbyvalalign(FLOAT8OID, &typlen, &typbyval, &typalign);
6090 for (
x = 0 + ngbwidth;
x < width - ngbwidth;
x++) {
6091 for(
y = 0 + ngbheight;
y < height - ngbheight;
y++) {
6096 pixelreplace =
false;
6100 pixelreplace =
true;
6103 for (u =
x - ngbwidth; u <=
x + ngbwidth; u++) {
6104 for (v =
y - ngbheight; v <=
y + ngbheight; v++) {
6109 neighborData[nIndex] = Float8GetDatum((
double)
r);
6110 neighborNulls[nIndex] =
false;
6111 nNodataOnly =
false;
6115 if (valuereplace && pixelreplace) {
6117 neighborData[nIndex] = Float8GetDatum((
double)rpix);
6118 neighborNulls[nIndex] =
false;
6123 neighborData[nIndex] = PointerGetDatum(NULL);
6124 neighborNulls[nIndex] =
true;
6131 neighborData[nIndex] = PointerGetDatum(NULL);
6132 neighborNulls[nIndex] =
true;
6144 if (!(nNodataOnly ||
6145 (nNullSkip && nNullItems > 0) ||
6146 (valuereplace && nNullItems > 0))) {
6148 x,
y, winwidth, winheight);
6150 neighborDatum = construct_md_array((
void *)neighborData, neighborNulls, 2, neighborDims, neighborLbs,
6151 FLOAT8OID, typlen, typbyval, typalign);
6153 #if POSTGIS_PGSQL_VERSION < 120
6155 cbdata.arg[0] = PointerGetDatum(neighborDatum);
6158 tmpnewval = FunctionCallInvoke(&cbdata);
6161 if (cbdata.isnull) {
6162 newval = newnodatavalue;
6166 cbdata->args[0].value = PointerGetDatum(neighborDatum);
6169 tmpnewval = FunctionCallInvoke(cbdata);
6174 newval = newnodatavalue;
6178 newval = DatumGetFloat8(tmpnewval);
6194 pfree(neighborNulls);
6195 pfree(neighborData);
6197 pfree(txtCallbackParam);
6200 PG_FREE_IF_COPY(pgraster, 0);
6204 POSTGIS_RT_DEBUG(3,
"RASTER_mapAlgebraFctNgb: raster modified, serializing it.");
6215 SET_VARSIZE(pgrtn, pgrtn->
size);
6216 PG_RETURN_POINTER(pgrtn);
6219 #define ARGKWCOUNT 8
6229 int pgrastpos[2] = {-1, -1};
6232 int _isempty[2] = {0};
6236 int _hasnodata[2] = {0};
6237 double _nodataval[2] = {0};
6238 double _offset[4] = {0.};
6239 double _rastoffset[2][4] = {{0.}};
6240 int _haspixel[2] = {0};
6241 double _pixel[2] = {0};
6242 int _pos[2][2] = {{0}};
6243 uint16_t _dim[2][2] = {{0}};
6245 char *pixtypename = NULL;
6247 char *extenttypename = NULL;
6252 uint16_t dim[2] = {0};
6255 double nodataval = 0;
6256 double gt[6] = {0.};
6258 Oid calltype = InvalidOid;
6261 uint16_t spi_exprpos[3] = {4, 7, 8};
6265 SPIPlanPtr spi_plan[3] = {NULL};
6266 uint16_t spi_empty = 0;
6267 Oid *argtype = NULL;
6269 char *argkw[] = {
"[rast1.x]",
"[rast1.y]",
"[rast1.val]",
"[rast1]",
"[rast2.x]",
"[rast2.y]",
"[rast2.val]",
"[rast2]"};
6273 SPITupleTable *tuptable = NULL;
6276 bool isnull =
FALSE;
6277 int hasargval[3] = {0};
6278 double argval[3] = {0.};
6279 int hasnodatanodataval = 0;
6280 double nodatanodataval = 0;
6283 Oid ufc_noid = InvalidOid;
6285 #if POSTGIS_PGSQL_VERSION < 120
6286 FunctionCallInfoData ufc_info;
6288 LOCAL_FCINFO(ufc_info, FUNC_MAX_ARGS);
6290 int ufc_nullcount = 0;
6306 for (i = 0, j = 0; i < set_count; i++) {
6307 if (!PG_ARGISNULL(j)) {
6308 pgrast[i] = (
rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(j));
6315 for (k = 0; k <= i; k++) {
6316 if (k < i &&
rast[k] != NULL)
6318 if (pgrastpos[k] != -1)
6319 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6321 elog(ERROR,
"RASTER_mapAlgebra2: Could not deserialize the %s raster", i < 1 ?
"first" :
"second");
6329 if (!PG_ARGISNULL(j)) {
6330 bandindex[i] = PG_GETARG_INT32(j);
6343 if (
rast[0] == NULL &&
rast[1] == NULL) {
6344 elog(NOTICE,
"The two rasters provided are NULL. Returning NULL");
6345 for (k = 0; k < set_count; k++) {
6346 if (pgrastpos[k] != -1)
6347 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6353 if (_isempty[0] && _isempty[1]) {
6354 elog(NOTICE,
"The two rasters provided are empty. Returning empty raster");
6358 for (k = 0; k < set_count; k++) {
6359 if (
rast[k] != NULL)
6361 if (pgrastpos[k] != -1)
6362 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6364 elog(ERROR,
"RASTER_mapAlgebra2: Could not create empty raster");
6374 SET_VARSIZE(pgrtn, pgrtn->
size);
6375 PG_RETURN_POINTER(pgrtn);
6380 (
rast[0] == NULL || _isempty[0]) ||
6381 (
rast[1] == NULL || _isempty[1])
6384 if (
rast[0] == NULL || _isempty[0]) {
6397 if (_rast[i] != NULL)
6409 if (_rast[i] == NULL) {
6411 for (k = 0; k < set_count; k++) {
6412 if (pgrastpos[k] != -1)
6413 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6415 elog(ERROR,
"RASTER_mapAlgebra2: Could not create NODATA raster");
6449 for (k = 0; k < set_count; k++) {
6450 if (_rast[k] != NULL)
6452 if (pgrastpos[k] != -1)
6453 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6455 elog(ERROR,
"RASTER_mapAlgebra2: Could not test for alignment on the two rasters");
6459 elog(NOTICE,
"The two rasters provided do not have the same alignment. Returning NULL");
6460 for (k = 0; k < set_count; k++) {
6461 if (_rast[k] != NULL)
6463 if (pgrastpos[k] != -1)
6464 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6470 if (!PG_ARGISNULL(5)) {
6474 if (pixtype ==
PT_END ) {
6475 for (k = 0; k < set_count; k++) {
6476 if (_rast[k] != NULL)
6478 if (pgrastpos[k] != -1)
6479 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6481 elog(ERROR,
"RASTER_mapAlgebra2: Invalid pixel type: %s", pixtypename);
6487 if (!PG_ARGISNULL(6)) {
6500 for (k = 0; k < set_count; k++) {
6501 if (_rast[k] != NULL)
6503 if (pgrastpos[k] != -1)
6504 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6506 elog(ERROR,
"RASTER_mapAlgebra2: Could not get output raster of correct extent");
6511 _rastoffset[0][0] = _offset[0];
6512 _rastoffset[0][1] = _offset[1];
6513 _rastoffset[1][0] = _offset[2];
6514 _rastoffset[1][1] = _offset[3];
6522 switch (extenttype) {
6532 (extenttype ==
ET_FIRST && i == 0) ||
6536 elog(NOTICE,
"The %s raster is NULL. Returning NULL", (i != 1 ?
"FIRST" :
"SECOND"));
6537 for (k = 0; k < set_count; k++) {
6538 if (_rast[k] != NULL)
6540 if (pgrastpos[k] != -1)
6541 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6549 elog(NOTICE,
"The %s raster does not have the band at index %d. Returning no band raster of correct extent",
6550 (i != 1 ?
"FIRST" :
"SECOND"), bandindex[i]
6553 for (k = 0; k < set_count; k++) {
6554 if (_rast[k] != NULL)
6556 if (pgrastpos[k] != -1)
6557 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6562 if (!pgrtn) PG_RETURN_NULL();
6564 SET_VARSIZE(pgrtn, pgrtn->
size);
6565 PG_RETURN_POINTER(pgrtn);
6573 _isempty[0] || _isempty[1] ||
6576 elog(NOTICE,
"The two rasters provided have no intersection. Returning no band raster");
6579 if (dim[0] || dim[1]) {
6584 for (k = 0; k < set_count; k++) {
6585 if (_rast[k] != NULL)
6587 if (pgrastpos[k] != -1)
6588 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6590 elog(ERROR,
"RASTER_mapAlgebra2: Could not create no band raster");
6598 for (k = 0; k < set_count; k++) {
6599 if (_rast[k] != NULL)
6601 if (pgrastpos[k] != -1)
6602 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6607 if (!pgrtn) PG_RETURN_NULL();
6609 SET_VARSIZE(pgrtn, pgrtn->
size);
6610 PG_RETURN_POINTER(pgrtn);
6615 for (k = 0; k < set_count; k++) {
6616 if (_rast[k] != NULL)
6618 if (pgrastpos[k] != -1)
6619 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6621 elog(ERROR,
"RASTER_mapAlgebra2: ET_LAST and ET_CUSTOM are not implemented");
6631 elog(NOTICE,
"The two rasters provided do not have the respectively specified band indices. Returning no band raster of correct extent");
6633 for (k = 0; k < set_count; k++) {
6634 if (_rast[k] != NULL)
6636 if (pgrastpos[k] != -1)
6637 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6642 if (!pgrtn) PG_RETURN_NULL();
6644 SET_VARSIZE(pgrtn, pgrtn->
size);
6645 PG_RETURN_POINTER(pgrtn);
6649 for (i = 0; i < set_count; i++) {
6658 if (_band[i] == NULL) {
6659 for (k = 0; k < set_count; k++) {
6660 if (_rast[k] != NULL)
6662 if (pgrastpos[k] != -1)
6663 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6666 elog(ERROR,
"RASTER_mapAlgebra2: Could not get band %d of the %s raster",
6668 (i < 1 ?
"FIRST" :
"SECOND")
6680 if ((extenttype ==
ET_SECOND && !_isempty[1]) || _isempty[0])
6687 if (extenttype ==
ET_SECOND && !_isempty[1] && _hasnodata[1]) {
6688 nodataval = _nodataval[1];
6690 else if (!_isempty[0] && _hasnodata[0]) {
6691 nodataval = _nodataval[0];
6693 else if (!_isempty[1] && _hasnodata[1]) {
6694 nodataval = _nodataval[1];
6697 elog(NOTICE,
"Neither raster provided has a NODATA value for the specified band indices. NODATA value set to minimum possible for %s",
rt_pixtype_name(pixtype));
6709 for (k = 0; k < set_count; k++) {
6710 if (_rast[k] != NULL)
6712 if (pgrastpos[k] != -1)
6713 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6716 elog(ERROR,
"RASTER_mapAlgebra2: Could not add new band to output raster");
6723 for (k = 0; k < set_count; k++) {
6724 if (_rast[k] != NULL)
6726 if (pgrastpos[k] != -1)
6727 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6730 elog(ERROR,
"RASTER_mapAlgebra2: Could not get newly added band of output raster");
6735 (
int) _rastoffset[0][0],
6736 (
int) _rastoffset[0][1],
6737 (
int) _rastoffset[1][0],
6738 (
int) _rastoffset[1][1]
6741 POSTGIS_RT_DEBUGF(4,
"metadata = (%f, %f, %d, %d, %f, %f, %f, %f, %d)",
6758 calltype = get_fn_expr_argtype(fcinfo->flinfo, 4);
6765 if (SPI_connect() != SPI_OK_CONNECT) {
6766 for (k = 0; k < set_count; k++) {
6767 if (_rast[k] != NULL)
6769 if (pgrastpos[k] != -1)
6770 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6773 elog(ERROR,
"RASTER_mapAlgebra2: Could not connect to the SPI manager");
6778 memset(hasargval, 0,
sizeof(
int) * spi_count);
6788 for (i = 0; i < spi_count; i++) {
6789 if (!PG_ARGISNULL(spi_exprpos[i])) {
6791 char place[5] =
"$1";
6806 sprintf(place,
"$%d", k);
6812 len = strlen(
"SELECT (") + strlen(expr) + strlen(
")::double precision");
6813 sql = (
char *) palloc(len + 1);
6816 for (k = 0; k < spi_count; k++) SPI_freeplan(spi_plan[k]);
6819 for (k = 0; k < set_count; k++) {
6820 if (_rast[k] != NULL)
6822 if (pgrastpos[k] != -1)
6823 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6827 elog(ERROR,
"RASTER_mapAlgebra2: Could not allocate memory for expression parameter %d", spi_exprpos[i]);
6831 strncpy(
sql,
"SELECT (", strlen(
"SELECT ("));
6832 strncpy(
sql + strlen(
"SELECT ("), expr, strlen(expr));
6833 strncpy(
sql + strlen(
"SELECT (") + strlen(expr),
")::double precision", strlen(
")::double precision"));
6839 if (spi_argcount[i]) {
6840 argtype = (Oid *) palloc(spi_argcount[i] *
sizeof(Oid));
6841 if (argtype == NULL) {
6844 for (k = 0; k < spi_count; k++) SPI_freeplan(spi_plan[k]);
6847 for (k = 0; k < set_count; k++) {
6848 if (_rast[k] != NULL)
6850 if (pgrastpos[k] != -1)
6851 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6855 elog(ERROR,
"RASTER_mapAlgebra2: Could not allocate memory for prepared plan argtypes of expression parameter %d", spi_exprpos[i]);
6861 if (argpos[i][j] < 1)
continue;
6865 (strstr(argkw[j],
"[rast1.x]") != NULL) ||
6866 (strstr(argkw[j],
"[rast1.y]") != NULL) ||
6867 (strstr(argkw[j],
"[rast2.x]") != NULL) ||
6868 (strstr(argkw[j],
"[rast2.y]") != NULL)
6870 argtype[k] = INT4OID;
6874 argtype[k] = FLOAT8OID;
6880 spi_plan[i] = SPI_prepare(
sql, spi_argcount[i], argtype);
6883 if (spi_plan[i] == NULL) {
6886 for (k = 0; k < spi_count; k++) SPI_freeplan(spi_plan[k]);
6889 for (k = 0; k < set_count; k++) {
6890 if (_rast[k] != NULL)
6892 if (pgrastpos[k] != -1)
6893 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6897 elog(ERROR,
"RASTER_mapAlgebra2: Could not create prepared plan of expression parameter %d", spi_exprpos[i]);
6903 err = SPI_execute(
sql,
TRUE, 0);
6904 if (err != SPI_OK_SELECT || SPI_tuptable == NULL || SPI_processed != 1) {
6907 for (k = 0; k < spi_count; k++) SPI_freeplan(spi_plan[k]);
6910 for (k = 0; k < set_count; k++) {
6911 if (_rast[k] != NULL)
6913 if (pgrastpos[k] != -1)
6914 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6918 elog(ERROR,
"RASTER_mapAlgebra2: Could not evaluate expression parameter %d", spi_exprpos[i]);
6923 tupdesc = SPI_tuptable->tupdesc;
6924 tuptable = SPI_tuptable;
6925 tuple = tuptable->vals[0];
6927 datum = SPI_getbinval(tuple, tupdesc, 1, &isnull);
6928 if (SPI_result == SPI_ERROR_NOATTRIBUTE) {
6931 if (SPI_tuptable) SPI_freetuptable(tuptable);
6932 for (k = 0; k < spi_count; k++) SPI_freeplan(spi_plan[k]);
6935 for (k = 0; k < set_count; k++) {
6936 if (_rast[k] != NULL)
6938 if (pgrastpos[k] != -1)
6939 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6943 elog(ERROR,
"RASTER_mapAlgebra2: Could not get result of expression parameter %d", spi_exprpos[i]);
6949 argval[i] = DatumGetFloat8(datum);
6952 if (SPI_tuptable) SPI_freetuptable(tuptable);
6962 if (!PG_ARGISNULL(9)) {
6963 hasnodatanodataval = 1;
6964 nodatanodataval = PG_GETARG_FLOAT8(9);
6967 hasnodatanodataval = 0;
6970 case REGPROCEDUREOID: {
6972 if (!PG_ARGISNULL(4)) {
6975 ufc_noid = PG_GETARG_OID(4);
6978 fmgr_info(ufc_noid, &ufl_info);
6982 if (ufl_info.fn_retset) {
6986 else if (ufl_info.fn_nargs < 3 || ufl_info.fn_nargs > 4) {
6996 for (k = 0; k < set_count; k++) {
6997 if (_rast[k] != NULL)
6999 if (pgrastpos[k] != -1)
7000 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
7005 elog(ERROR,
"RASTER_mapAlgebra2: Function provided must have three or four input parameters");
7007 elog(ERROR,
"RASTER_mapAlgebra2: Function provided must return double precision not resultset");
7011 if (func_volatile(ufc_noid) ==
'v') {
7012 elog(NOTICE,
"Function provided is VOLATILE. Unless required and for best performance, function should be IMMUTABLE or STABLE");
7016 #if POSTGIS_PGSQL_VERSION < 120
7017 InitFunctionCallInfoData(ufc_info, &ufl_info, ufl_info.fn_nargs, InvalidOid, NULL, NULL);
7018 memset(ufc_info.argnull,
FALSE,
sizeof(
bool) * ufl_info.fn_nargs);
7020 InitFunctionCallInfoData(
7021 *ufc_info, &ufl_info, ufl_info.fn_nargs, InvalidOid, NULL, NULL);
7022 ufc_info->args[0].isnull =
FALSE;
7023 ufc_info->args[1].isnull =
FALSE;
7024 ufc_info->args[2].isnull =
FALSE;
7025 if (ufl_info.fn_nargs == 4)
7026 ufc_info->args[3].isnull =
FALSE;
7029 if (ufl_info.fn_nargs != 4)
7033 #if POSTGIS_PGSQL_VERSION < 120
7034 if (!PG_ARGISNULL(7)) {
7035 ufc_info.arg[k] = PG_GETARG_DATUM(7);
7038 ufc_info.arg[k] = (Datum) NULL;
7039 ufc_info.argnull[k] =
TRUE;
7043 if (!PG_ARGISNULL(7))
7045 ufc_info->args[k].value = PG_GETARG_DATUM(7);
7049 ufc_info->args[k].value = (Datum)NULL;
7050 ufc_info->args[k].isnull =
TRUE;
7058 for (k = 0; k < set_count; k++) {
7059 if (_rast[k] != NULL)
7061 if (pgrastpos[k] != -1)
7062 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
7065 elog(ERROR,
"RASTER_mapAlgebra2: Invalid data type for expression or userfunction");
7073 (calltype == TEXTOID) && (
7074 (spi_empty != spi_count) || hasnodatanodataval
7077 (calltype == REGPROCEDUREOID) && (ufc_noid != InvalidOid)
7079 for (
x = 0;
x < dim[0];
x++) {
7080 for (
y = 0;
y < dim[1];
y++) {
7083 for (i = 0; i < set_count; i++) {
7088 _x =
x - (int) _rastoffset[i][0];
7089 _y =
y - (int) _rastoffset[i][1];
7092 _pos[i][0] = _x + 1;
7093 _pos[i][1] = _y + 1;
7096 if (_band[i] == NULL) {
7097 if (!_hasnodata[i]) {
7099 _pixel[i] = _nodataval[i];
7104 (_x >= 0 && _x < _dim[i][0]) &&
7105 (_y >= 0 && _y < _dim[i][1])
7110 if (calltype == TEXTOID) {
7111 for (k = 0; k < spi_count; k++) SPI_freeplan(spi_plan[k]);
7115 for (k = 0; k < set_count; k++) {
7116 if (_rast[k] != NULL)
7118 if (pgrastpos[k] != -1)
7119 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
7123 elog(ERROR,
"RASTER_mapAlgebra2: Could not get pixel of %s raster", (i < 1 ?
"FIRST" :
"SECOND"));
7127 if (!_hasnodata[i] || !isnodata)
7146 if (!_haspixel[0] && !_haspixel[1])
7150 else if (_haspixel[0] && !_haspixel[1])
7154 else if (!_haspixel[0] && _haspixel[1])
7163 if (hasnodatanodataval) {
7165 pixel = nodatanodataval;
7169 else if (hasargval[i]) {
7174 else if (spi_plan[i] != NULL) {
7178 memset(values, (Datum) NULL,
sizeof(Datum) *
ARGKWCOUNT);
7183 if (spi_argcount[i]) {
7187 if (idx < 1)
continue;
7190 if (strstr(argkw[j],
"[rast1.x]") != NULL) {
7191 values[idx] = _pos[0][0];
7193 else if (strstr(argkw[j],
"[rast1.y]") != NULL) {
7194 values[idx] = _pos[0][1];
7197 (strstr(argkw[j],
"[rast1.val]") != NULL) ||
7198 (strstr(argkw[j],
"[rast1]") != NULL)
7200 if (_isempty[0] || !_haspixel[0])
7203 values[idx] = Float8GetDatum(_pixel[0]);
7205 else if (strstr(argkw[j],
"[rast2.x]") != NULL) {
7206 values[idx] = _pos[1][0];
7208 else if (strstr(argkw[j],
"[rast2.y]") != NULL) {
7209 values[idx] = _pos[1][1];
7212 (strstr(argkw[j],
"[rast2.val]") != NULL) ||
7213 (strstr(argkw[j],
"[rast2]") != NULL)
7215 if (_isempty[1] || !_haspixel[1])
7218 values[idx] = Float8GetDatum(_pixel[1]);
7224 err = SPI_execute_plan(spi_plan[i], values, nulls,
TRUE, 1);
7225 if (err != SPI_OK_SELECT || SPI_tuptable == NULL || SPI_processed != 1) {
7227 for (k = 0; k < spi_count; k++) SPI_freeplan(spi_plan[k]);
7230 for (k = 0; k < set_count; k++) {
7231 if (_rast[k] != NULL)
7233 if (pgrastpos[k] != -1)
7234 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
7238 elog(ERROR,
"RASTER_mapAlgebra2: Unexpected error when running prepared statement %d", i);
7243 tupdesc = SPI_tuptable->tupdesc;
7244 tuptable = SPI_tuptable;
7245 tuple = tuptable->vals[0];
7247 datum = SPI_getbinval(tuple, tupdesc, 1, &isnull);
7248 if (SPI_result == SPI_ERROR_NOATTRIBUTE) {
7250 if (SPI_tuptable) SPI_freetuptable(tuptable);
7251 for (k = 0; k < spi_count; k++) SPI_freeplan(spi_plan[k]);
7254 for (k = 0; k < set_count; k++) {
7255 if (_rast[k] != NULL)
7257 if (pgrastpos[k] != -1)
7258 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
7262 elog(ERROR,
"RASTER_mapAlgebra2: Could not get result of prepared statement %d", i);
7268 pixel = DatumGetFloat8(datum);
7271 if (SPI_tuptable) SPI_freetuptable(tuptable);
7274 case REGPROCEDUREOID: {
7279 for (i = 0; i < set_count; i++) {
7280 #if POSTGIS_PGSQL_VERSION < 120
7281 ufc_info.arg[i] = Float8GetDatum(_pixel[i]);
7283 ufc_info->args[i].value = Float8GetDatum(_pixel[i]);
7287 #if POSTGIS_PGSQL_VERSION < 120
7288 ufc_info.argnull[i] =
FALSE;
7290 ufc_info->args[i].isnull =
FALSE;
7295 #if POSTGIS_PGSQL_VERSION < 120
7296 ufc_info.argnull[i] =
TRUE;
7298 ufc_info->args[i].isnull =
TRUE;
7306 if (ufl_info.fn_strict && ufc_nullcount)
7310 if (ufl_info.fn_nargs == 4) {
7313 for (i = 0; i < set_count; i++) {
7315 d[0] = Int32GetDatum(_pos[i][0]);
7316 d[1] = Int32GetDatum(_pos[i][1]);
7319 d[2] = Int32GetDatum(_pos[i][0]);
7320 d[3] = Int32GetDatum(_pos[i][1]);
7324 a = construct_array(d, 4, INT4OID,
sizeof(
int32),
true,
'i');
7325 #if POSTGIS_PGSQL_VERSION < 120
7326 ufc_info.arg[2] = PointerGetDatum(a);
7327 ufc_info.argnull[2] =
FALSE;
7329 ufc_info->args[2].value = PointerGetDatum(a);
7330 ufc_info->args[2].isnull =
FALSE;
7334 #if POSTGIS_PGSQL_VERSION < 120
7335 datum = FunctionCallInvoke(&ufc_info);
7338 if (!ufc_info.isnull) {
7340 pixel = DatumGetFloat8(datum);
7343 datum = FunctionCallInvoke(ufc_info);
7346 if (!ufc_info->isnull)
7349 pixel = DatumGetFloat8(datum);
7359 if (calltype == TEXTOID) {
7360 for (k = 0; k < spi_count; k++) SPI_freeplan(spi_plan[k]);
7364 for (k = 0; k < set_count; k++) {
7365 if (_rast[k] != NULL)
7367 if (pgrastpos[k] != -1)
7368 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
7372 elog(ERROR,
"RASTER_mapAlgebra2: Could not set pixel value of output raster");
7384 if (calltype == TEXTOID) {
7385 for (i = 0; i < spi_count; i++) {
7386 if (spi_plan[i] != NULL) SPI_freeplan(spi_plan[i]);
7391 for (k = 0; k < set_count; k++) {
7392 if (_rast[k] != NULL)
7394 if (pgrastpos[k] != -1)
7395 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
7400 if (!pgrtn) PG_RETURN_NULL();
7404 SET_VARSIZE(pgrtn, pgrtn->
size);
7405 PG_RETURN_POINTER(pgrtn);
int32_t gserialized_get_srid(const GSERIALIZED *s)
Extract the SRID from the serialized form (it is packed into three bytes so this is a handy function)...
LWGEOM * lwgeom_from_gserialized(const GSERIALIZED *g)
Allocate a new LWGEOM from a GSERIALIZED.
int lwgeom_ndims(const LWGEOM *geom)
Return the number of dimensions (2, 3, 4) in a geometry.
void lwgeom_free(LWGEOM *geom)
LWGEOM * lwgeom_intersection(const LWGEOM *geom1, const LWGEOM *geom2)
uint8_t * lwgeom_to_wkb(const LWGEOM *geom, uint8_t variant, size_t *size_out)
Convert LWGEOM to a char* in WKB format.
int clamp_srid(int srid)
Return a valid SRID from an arbitrary integer Raises a notice if what comes out is different from wha...
int lwgeom_is_empty(const LWGEOM *geom)
Return true or false depending on whether a geometry is an "empty" geometry (no vertices members)
#define SRID_UNKNOWN
Unknown SRID value.
LWGEOM * lwgeom_force_2d(const LWGEOM *geom)
Strip out the Z/M components of an LWGEOM.
rt_band rt_band_reclass(rt_band srcband, rt_pixtype pixtype, uint32_t hasnodata, double nodataval, rt_reclassexpr *exprset, int exprcount)
Returns new band with values reclassified.
#define RASTER_DEBUG(level, msg)
int32_t rt_raster_get_srid(rt_raster raster)
Get raster's SRID.
double rt_raster_get_x_skew(rt_raster raster)
Get skew about the X axis.
double rt_raster_get_x_offset(rt_raster raster)
Get raster x offset, in projection units.
int rt_raster_generate_new_band(rt_raster raster, rt_pixtype pixtype, double initialvalue, uint32_t hasnodata, double nodatavalue, int index)
Generate a new inline band and add it to a raster.
void rt_raster_set_geotransform_matrix(rt_raster raster, double *gt)
Set raster's geotransform using 6-element array.
rt_raster rt_raster_colormap(rt_raster raster, int nband, rt_colormap colormap)
Returns a new raster with up to four 8BUI bands (RGBA) from applying a colormap to the user-specified...
void rt_raster_set_scale(rt_raster raster, double scaleX, double scaleY)
Set scale in projection units.
int rt_raster_add_band(rt_raster raster, rt_band band, int index)
Add band data to a raster.
int rt_band_get_hasnodata_flag(rt_band band)
Get hasnodata flag value.
rt_pixtype rt_pixtype_index_from_name(const char *pixname)
rt_errorstate rt_band_get_pixel(rt_band band, int x, int y, double *value, int *nodata)
Get pixel value.
void rt_raster_destroy(rt_raster raster)
Release memory associated to a raster.
void rt_raster_set_skews(rt_raster raster, double skewX, double skewY)
Set skews about the X and Y axis.
int rt_band_get_isnodata_flag(rt_band band)
Get isnodata flag value.
rt_raster rt_raster_new(uint32_t width, uint32_t height)
Construct a raster with given dimensions.
rt_extenttype rt_util_extent_type(const char *name)
double rt_pixtype_get_min_value(rt_pixtype pixtype)
Return minimum value possible for pixel type.
int rt_raster_has_band(rt_raster raster, int nband)
Return TRUE if the raster has a band of this number.
void * rt_raster_serialize(rt_raster raster)
Return this raster in serialized form.
double rt_raster_get_x_scale(rt_raster raster)
Get scale X in projection units.
rt_band rt_raster_replace_band(rt_raster raster, rt_band band, int index)
Replace band at provided index with new band.
double rt_band_get_min_value(rt_band band)
Returns the minimal possible value for the band according to the pixel type.
rt_errorstate rt_band_set_pixel(rt_band band, int x, int y, double val, int *converted)
Set single pixel's value.
rt_errorstate rt_band_set_pixel_line(rt_band band, int x, int y, void *vals, uint32_t len)
Set values of multiple pixels.
rt_raster rt_raster_gdal_rasterize(const unsigned char *wkb, uint32_t wkb_len, const char *srs, uint32_t num_bands, rt_pixtype *pixtype, double *init, double *value, double *nodata, uint8_t *hasnodata, int *width, int *height, double *scale_x, double *scale_y, double *ul_xw, double *ul_yw, double *grid_xw, double *grid_yw, double *skew_x, double *skew_y, char **options)
Return a raster of the provided geometry.
void rt_band_destroy(rt_band band)
Destroy a raster band.
uint16_t rt_raster_get_num_bands(rt_raster raster)
rt_raster rt_raster_clone(rt_raster raster, uint8_t deep)
Clone an existing raster.
rt_errorstate rt_raster_get_convex_hull(rt_raster raster, LWGEOM **hull)
Get raster's convex hull.
uint16_t rt_raster_get_height(rt_raster raster)
rt_errorstate rt_raster_iterator(rt_iterator itrset, uint16_t itrcount, rt_extenttype extenttype, rt_raster customextent, rt_pixtype pixtype, uint8_t hasnodata, double nodataval, uint16_t distancex, uint16_t distancey, rt_mask mask, void *userarg, int(*callback)(rt_iterator_arg arg, void *userarg, double *value, int *nodata), rt_raster *rtnraster)
n-raster iterator.
void rt_raster_set_srid(rt_raster raster, int32_t srid)
Set raster's SRID.
rt_raster rt_raster_from_band(rt_raster raster, uint32_t *bandNums, int count)
Construct a new rt_raster from an existing rt_raster and an array of band numbers.
rt_errorstate rt_band_get_nodata(rt_band band, double *nodata)
Get NODATA value.
rt_pixtype rt_band_get_pixtype(rt_band band)
Return pixeltype of this band.
const char * rt_pixtype_name(rt_pixtype pixtype)
uint16_t rt_raster_get_width(rt_raster raster)
rt_errorstate rt_raster_from_two_rasters(rt_raster rast1, rt_raster rast2, rt_extenttype extenttype, rt_raster *rtnraster, double *offset)
int rt_util_same_geotransform_matrix(double *gt1, double *gt2)
int rt_raster_copy_band(rt_raster torast, rt_raster fromrast, int fromindex, int toindex)
Copy one band from one raster to another.
double rt_raster_get_y_scale(rt_raster raster)
Get scale Y in projection units.
rt_errorstate rt_raster_same_alignment(rt_raster rast1, rt_raster rast2, int *aligned, char **reason)
double rt_raster_get_y_skew(rt_raster raster)
Get skew about the Y axis.
void rt_raster_get_geotransform_matrix(rt_raster raster, double *gt)
Get 6-element array of raster geotransform matrix.
void rt_raster_set_offsets(rt_raster raster, double x, double y)
Set insertion points in projection units.
rt_bandstats rt_band_get_summary_stats(rt_band band, int exclude_nodata_value, double sample, int inc_vals, uint64_t *cK, double *cM, double *cQ)
Compute summary statistics for a band.
rt_raster rt_raster_deserialize(void *serialized, int header_only)
Return a raster from a serialized form.
rt_errorstate rt_band_get_pixel_line(rt_band band, int x, int y, uint16_t len, void **vals, uint16_t *nvals)
Get values of multiple pixels.
double rt_raster_get_y_offset(rt_raster raster)
Get raster y offset, in projection units.
int rt_raster_is_empty(rt_raster raster)
Return TRUE if the raster is empty.
rt_band rt_raster_get_band(rt_raster raster, int bandNum)
Return Nth band, or NULL if unavailable.
raster
Be careful!! Zeros function's input parameter can be a (height x width) array, not (width x height): ...
char * text_to_cstring(const text *textptr)
char * rtpg_removespaces(char *str)
char * rtpg_strtoupper(char *str)
char * rtpg_chartrim(const char *input, char *remove)
char * rtpg_strrstr(const char *s1, const char *s2)
char * rtpg_trim(const char *input)
char ** rtpg_strsplit(const char *str, const char *delimiter, uint32_t *n)
char * rtpg_strreplace(const char *str, const char *oldstr, const char *newstr, int *count)
Datum RASTER_nMapAlgebra(PG_FUNCTION_ARGS)
struct rtpg_colormap_arg_t * rtpg_colormap_arg
static rtpg_nmapalgebra_arg rtpg_nmapalgebra_arg_init()
static int rtpg_union_mean_callback(rt_iterator_arg arg, void *userarg, double *value, int *nodata)
Datum RASTER_mapAlgebraFct(PG_FUNCTION_ARGS)
static void rtpg_nmapalgebra_arg_destroy(rtpg_nmapalgebra_arg arg)
static void rtpg_nmapalgebraexpr_arg_destroy(rtpg_nmapalgebraexpr_arg arg)
static rtpg_nmapalgebraexpr_arg rtpg_nmapalgebraexpr_arg_init(int cnt, char **kw)
struct rtpg_clip_arg_t * rtpg_clip_arg
static void rtpg_colormap_arg_destroy(rtpg_colormap_arg arg)
Datum RASTER_colorMap(PG_FUNCTION_ARGS)
static rtpg_union_type rtpg_uniontype_index_from_name(const char *cutype)
struct rtpg_clip_band_t * rtpg_clip_band
static int rtpg_union_range_callback(rt_iterator_arg arg, void *userarg, double *value, int *nodata)
PG_FUNCTION_INFO_V1(RASTER_nMapAlgebra)
static int rtpg_nmapalgebra_rastbandarg_process(rtpg_nmapalgebra_arg arg, ArrayType *array, int *allnull, int *allempty, int *noband)
Datum RASTER_mapAlgebra2(PG_FUNCTION_ARGS)
static int rtpg_nmapalgebraexpr_callback(rt_iterator_arg arg, void *userarg, double *value, int *nodata)
static void rtpg_union_arg_destroy(rtpg_union_arg arg)
struct rtpg_union_arg_t * rtpg_union_arg
Datum RASTER_mapAlgebraExpr(PG_FUNCTION_ARGS)
static int rtpg_union_unionarg_process(rtpg_union_arg arg, ArrayType *array)
struct rtpg_nmapalgebra_arg_t * rtpg_nmapalgebra_arg
static rtpg_clip_arg rtpg_clip_arg_init()
Datum RASTER_mapAlgebraFctNgb(PG_FUNCTION_ARGS)
static int rtpg_union_callback(rt_iterator_arg arg, void *userarg, double *value, int *nodata)
static int rtpg_nmapalgebra_callback(rt_iterator_arg arg, void *userarg, double *value, int *nodata)
Datum RASTER_union_transfn(PG_FUNCTION_ARGS)
Datum RASTER_reclass(PG_FUNCTION_ARGS)
Datum RASTER_nMapAlgebraExpr(PG_FUNCTION_ARGS)
static void rtpg_clip_arg_destroy(rtpg_clip_arg arg)
static int rtpg_clip_callback(rt_iterator_arg arg, void *userarg, double *value, int *nodata)
struct rtpg_union_band_arg_t * rtpg_union_band_arg
struct rtpg_nmapalgebraexpr_arg_t * rtpg_nmapalgebraexpr_arg
static rtpg_colormap_arg rtpg_colormap_arg_init()
Datum RASTER_union_finalfn(PG_FUNCTION_ARGS)
Datum RASTER_clip(PG_FUNCTION_ARGS)
static int rtpg_union_noarg(rtpg_union_arg arg, rt_raster raster)
#define POSTGIS_RT_DEBUG(level, msg)
#define POSTGIS_RT_DEBUGF(level, msg,...)
enum rt_colormap_t::@9 method
struct rt_reclassexpr_t::rt_reclassrange src
struct rt_reclassexpr_t::rt_reclassrange dst
rtpg_nmapalgebra_callback_arg callback
FunctionCallInfoData ufc_info
rtpg_nmapalgebraexpr_callback_arg callback
rtpg_nmapalgebra_arg bandarg
struct rtpg_nmapalgebraexpr_callback_arg::@18 expr[3]
struct rtpg_nmapalgebraexpr_callback_arg::@20 kw
struct rtpg_nmapalgebraexpr_callback_arg::@19 nodatanodata
rtpg_union_band_arg bandarg
rtpg_union_type uniontype