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"
92 char fcinfo_data[SizeForFunctionCallInfo(FUNC_MAX_ARGS)];
97 #if defined(__clang__)
98 # pragma clang diagnostic pop
130 elog(ERROR,
"rtpg_nmapalgebra_arg_init: Could not allocate memory for arguments");
165 if (arg->
raster != NULL) {
182 if( arg->
mask != NULL )
205 if (arg == NULL || array == NULL) {
206 elog(ERROR,
"rtpg_nmapalgebra_rastbandarg_process: NULL values not permitted for parameters");
210 etype = ARR_ELEMTYPE(array);
211 get_typlenbyvalalign(etype, &typlen, &typbyval, &typalign);
216 typlen, typbyval, typalign,
221 elog(ERROR,
"rtpg_nmapalgebra_rastbandarg_process: Invalid argument for rastbandarg");
239 arg->
nband == NULL ||
243 elog(ERROR,
"rtpg_nmapalgebra_rastbandarg_process: Could not allocate memory for processing rastbandarg");
252 for (i = 0; i < n; i++) {
267 tup = (HeapTupleHeader) DatumGetPointer(e[i]);
269 elog(ERROR,
"rtpg_nmapalgebra_rastbandarg_process: Invalid argument for rastbandarg at index %d", i);
275 tupv = GetAttributeByName(tup,
"rast", &isnull);
277 elog(NOTICE,
"First argument (nband) of rastbandarg at index %d is NULL. Assuming NULL raster", i);
291 for (j = 0; j < i; j++) {
303 if (arg->
raster[i] == NULL) {
304 elog(ERROR,
"rtpg_nmapalgebra_rastbandarg_process: Could not deserialize raster at index %d", i);
320 tupv = GetAttributeByName(tup,
"nband", &isnull);
323 elog(NOTICE,
"First argument (nband) of rastbandarg at index %d is NULL. Assuming nband = %d", i,
nband);
326 nband = DatumGetInt32(tupv);
329 elog(ERROR,
"rtpg_nmapalgebra_rastbandarg_process: Band number provided for rastbandarg at index %d must be greater than zero (1-based)", i);
353 arg->
nband == NULL ||
356 elog(ERROR,
"rtpg_nmapalgebra_rastbandarg_process: Could not reallocate memory for processed rastbandarg");
371 double *
value,
int *nodata
379 ArrayType *mdValues = NULL;
380 Datum *_values = NULL;
381 bool *_nodata = NULL;
383 ArrayType *mdPos = NULL;
392 int lbound[3] = {1, 1, 1};
393 Datum datum = (Datum) NULL;
407 if (_values == NULL || _nodata == NULL) {
408 elog(ERROR,
"rtpg_nmapalgebra_callback: Could not allocate memory for values array");
415 for (z = 0; z < arg->
rasters; z++) {
417 for (
y = 0;
y < arg->
rows;
y++) {
423 _nodata[i] = (bool) arg->
nodata[z][
y][
x];
425 _values[i] = Float8GetDatum(arg->
values[z][
y][
x]);
427 _values[i] = (Datum) NULL;
435 get_typlenbyvalalign(FLOAT8OID, &typlen, &typbyval, &typalign);
438 mdValues = construct_md_array(
442 typlen, typbyval, typalign
447 _pos = palloc(
sizeof(Datum) * (arg->
rasters + 1) * 2);
448 _null = palloc(
sizeof(
bool) * (arg->
rasters + 1) * 2);
449 if (_pos == NULL || _null == NULL) {
451 elog(ERROR,
"rtpg_nmapalgebra_callback: Could not allocate memory for position array");
454 memset(_null, 0,
sizeof(
bool) * (arg->
rasters + 1) * 2);
463 for (z = 0; z < arg->
rasters; z++) {
464 _pos[i] = (Datum)arg->
src_pixel[z][0] + 1;
467 _pos[i] = (Datum)arg->
src_pixel[z][1] + 1;
472 get_typlenbyvalalign(INT4OID, &typlen, &typbyval, &typalign);
480 mdPos = construct_md_array(
484 typlen, typbyval, typalign
489 callback->
ufc_info->args[0].value = PointerGetDatum(mdValues);
490 callback->
ufc_info->args[1].value = PointerGetDatum(mdPos);
493 datum = FunctionCallInvoke(callback->
ufc_info);
502 *
value = DatumGetFloat8(datum);
505 *
value = (double) DatumGetFloat4(datum);
508 *
value = (double) DatumGetInt32(datum);
511 *
value = (double) DatumGetInt16(datum);
529 ArrayType *maskArray;
560 elog(ERROR,
"RASTER_nMapAlgebra: Could not initialize argument structure");
567 elog(ERROR,
"RASTER_nMapAlgebra: Could not process rastbandarg");
571 POSTGIS_RT_DEBUGF(4,
"allnull, allempty, noband = %d, %d, %d", allnull, allempty, noband);
575 elog(NOTICE,
"All input rasters are NULL. Returning NULL");
581 if (!PG_ARGISNULL(2)) {
582 char *pixtypename = text_to_cstring(PG_GETARG_TEXT_P(2));
588 elog(ERROR,
"RASTER_nMapAlgebra: Invalid pixel type: %s", pixtypename);
594 if (!PG_ARGISNULL(3)){
595 arg->
distance[0] = PG_GETARG_INT32(3);
600 if (!PG_ARGISNULL(4)){
601 arg->
distance[1] = PG_GETARG_INT32(4);
607 elog(ERROR,
"RASTER_nMapAlgebra: Distance for X and Y axis must be greater than or equal to zero");
612 if (!PG_ARGISNULL(5)) {
620 if (PG_ARGISNULL(6)) {
621 elog(NOTICE,
"Custom extent is NULL. Returning NULL");
632 elog(ERROR,
"RASTER_nMapAlgebra: Could not deserialize custom extent");
636 elog(NOTICE,
"Custom extent is an empty raster. Returning empty raster");
641 elog(ERROR,
"RASTER_nMapAlgebra: Could not create empty raster");
647 if (!pgraster) PG_RETURN_NULL();
649 SET_VARSIZE(pgraster, pgraster->
size);
650 PG_RETURN_POINTER(pgraster);
656 if( PG_ARGISNULL(7) ){
661 maskArray = PG_GETARG_ARRAYTYPE_P(7);
662 etype = ARR_ELEMTYPE(maskArray);
663 get_typlenbyvalalign(etype,&typlen,&typbyval,&typalign);
671 elog(ERROR,
"RASTER_nMapAlgebra: Mask data type must be FLOAT8 or FLOAT4");
675 ndims = ARR_NDIM(maskArray);
678 elog(ERROR,
"RASTER_nMapAlgebra: Mask Must be a 2D array");
683 maskDims = ARR_DIMS(maskArray);
685 if (maskDims[0] % 2 == 0 || maskDims[1] % 2 == 0) {
686 elog(ERROR,
"RASTER_nMapAlgebra: Mask dimensions must be odd");
694 typlen, typbyval,typalign,
695 &maskElements,&maskNulls,&num
698 if (num < 1 || num != (maskDims[0] * maskDims[1])) {
703 elog(ERROR,
"RASTER_nMapAlgebra: Could not deconstruct new values array");
709 arg->
mask->
values = palloc(
sizeof(
double*)* maskDims[0]);
710 arg->
mask->
nodata = palloc(
sizeof(
int*)*maskDims[0]);
711 for (i = 0; i < maskDims[0]; i++) {
712 arg->
mask->
values[i] = (
double*) palloc(
sizeof(
double) * maskDims[1]);
713 arg->
mask->
nodata[i] = (
int*) palloc(
sizeof(
int) * maskDims[1]);
718 for (
y = 0;
y < maskDims[0];
y++) {
719 for (
x = 0;
x < maskDims[1];
x++) {
727 arg->
mask->
values[
y][
x] = (double) DatumGetFloat4(maskElements[i]);
731 arg->
mask->
values[
y][
x] = (double) DatumGetFloat8(maskElements[i]);
742 if (maskDims[0] == 1 && maskDims[1] == 1) {
753 if (PG_ARGISNULL(8) || !PG_GETARG_BOOL(8)) {
754 if (arg->
mask != NULL)
765 elog(NOTICE,
"All input rasters are empty. Returning empty raster");
770 elog(NOTICE,
"All input rasters do not have bands at indicated indexes. Returning empty raster");
778 elog(ERROR,
"RASTER_nMapAlgebra: Could not create empty raster");
784 if (!pgraster) PG_RETURN_NULL();
786 SET_VARSIZE(pgraster, pgraster->
size);
787 PG_RETURN_POINTER(pgraster);
791 if (!PG_ARGISNULL(1) || get_fn_expr_argtype(fcinfo->flinfo, 1) == REGPROCEDUREOID) {
810 get_func_result_type(
837 elog(ERROR,
"RASTER_nMapAlgebra: Function provided must return a double precision, float, int or smallint");
840 elog(ERROR,
"RASTER_nMapAlgebra: Function provided must return scalar (double precision, float, int, smallint)");
843 elog(ERROR,
"RASTER_nMapAlgebra: Function provided must have three input parameters");
846 elog(ERROR,
"RASTER_nMapAlgebra: Function provided must return double precision, not resultset");
853 elog(NOTICE,
"Function provided is VOLATILE. Unless required and for best performance, function should be IMMUTABLE or STABLE");
867 if (!PG_ARGISNULL(9))
873 arg->
callback.
ufc_info->args[2].value = PointerGetDatum(construct_empty_array(TEXTOID));
884 elog(ERROR,
"RASTER_nMapAlgebra: callbackfunc must be provided");
927 if (itrset == NULL) {
929 elog(ERROR,
"RASTER_nMapAlgebra: Could not allocate memory for iterator arguments");
958 elog(ERROR,
"RASTER_nMapAlgebra: Could not run raster iterator function");
972 SET_VARSIZE(pgraster, pgraster->
size);
973 PG_RETURN_POINTER(pgraster);
1017 elog(ERROR,
"rtpg_nmapalgebraexpr_arg_init: Could not allocate memory for arguments");
1023 elog(ERROR,
"rtpg_nmapalgebraexpr_arg_init: Could not allocate memory for arg->bandarg");
1036 elog(ERROR,
"rtpg_nmapalgebraexpr_arg_init: Could not allocate memory for spi_argpos");
1067 double *
value,
int *nodata
1070 SPIPlanPtr plan = NULL;
1090 else if (arg->
nodata[0][0][0] && !arg->
nodata[1][0][0]) {
1100 else if (!arg->
nodata[0][0][0] && arg->
nodata[1][0][0]) {
1127 if (arg->
nodata[0][0][0]) {
1163 SPITupleTable *tuptable = NULL;
1166 bool isnull =
FALSE;
1171 memset(values, (Datum) NULL,
sizeof(Datum) * callback->
kw.
count);
1172 memset(nulls,
FALSE,
sizeof(
char) * callback->
kw.
count);
1177 for (i = 0; i < callback->
kw.
count; i++) {
1179 if (idx < 1)
continue;
1185 values[idx] = Int32GetDatum(arg->
src_pixel[0][0] + 1);
1189 values[idx] = Int32GetDatum(arg->
src_pixel[0][1] + 1);
1195 if (!arg->
nodata[0][0][0])
1196 values[idx] = Float8GetDatum(arg->
values[0][0][0]);
1203 values[idx] = Int32GetDatum(arg->
src_pixel[0][0] + 1);
1207 values[idx] = Int32GetDatum(arg->
src_pixel[0][1] + 1);
1213 if (!arg->
nodata[0][0][0])
1214 values[idx] = Float8GetDatum(arg->
values[0][0][0]);
1221 values[idx] = Int32GetDatum(arg->
src_pixel[1][0] + 1);
1225 values[idx] = Int32GetDatum(arg->
src_pixel[1][1] + 1);
1231 if (!arg->
nodata[1][0][0])
1232 values[idx] = Float8GetDatum(arg->
values[1][0][0]);
1242 err = SPI_execute_plan(plan, values, nulls,
TRUE, 1);
1243 if (err != SPI_OK_SELECT || SPI_tuptable == NULL || SPI_processed != 1) {
1244 elog(ERROR,
"rtpg_nmapalgebraexpr_callback: Unexpected error when running prepared statement %d",
id);
1249 tupdesc = SPI_tuptable->tupdesc;
1250 tuptable = SPI_tuptable;
1251 tuple = tuptable->vals[0];
1253 datum = SPI_getbinval(tuple, tupdesc, 1, &isnull);
1254 if (SPI_result == SPI_ERROR_NOATTRIBUTE) {
1255 if (SPI_tuptable) SPI_freetuptable(tuptable);
1256 elog(ERROR,
"rtpg_nmapalgebraexpr_callback: Could not get result of prepared statement %d",
id);
1261 *
value = DatumGetFloat8(datum);
1281 if (SPI_tuptable) SPI_freetuptable(tuptable);
1291 MemoryContext mainMemCtx = CurrentMemoryContext;
1294 uint16_t exprpos[3] = {1, 4, 5};
1308 SPITupleTable *tuptable = NULL;
1311 bool isnull =
FALSE;
1317 const int argkwcount = 12;
1333 if (PG_ARGISNULL(0))
1339 elog(ERROR,
"RASTER_nMapAlgebraExpr: Could not initialize argument structure");
1346 elog(ERROR,
"RASTER_nMapAlgebra: Could not process rastbandarg");
1350 POSTGIS_RT_DEBUGF(4,
"allnull, allempty, noband = %d, %d, %d", allnull, allempty, noband);
1354 elog(NOTICE,
"All input rasters are NULL. Returning NULL");
1366 if (!PG_ARGISNULL(2)) {
1367 char *pixtypename = text_to_cstring(PG_GETARG_TEXT_P(2));
1373 elog(ERROR,
"RASTER_nMapAlgebraExpr: Invalid pixel type: %s", pixtypename);
1380 if (!PG_ARGISNULL(3)) {
1386 if (numraster < 2) {
1387 elog(NOTICE,
"CUSTOM extent type not supported. Defaulting to FIRST");
1391 elog(NOTICE,
"CUSTOM extent type not supported. Defaulting to INTERSECTION");
1395 else if (numraster < 2)
1401 if (!PG_ARGISNULL(6)) {
1409 elog(NOTICE,
"All input rasters are empty. Returning empty raster");
1414 elog(NOTICE,
"All input rasters do not have bands at indicated indexes. Returning empty raster");
1422 elog(ERROR,
"RASTER_nMapAlgebraExpr: Could not create empty raster");
1428 if (!pgraster) PG_RETURN_NULL();
1430 SET_VARSIZE(pgraster, pgraster->
size);
1431 PG_RETURN_POINTER(pgraster);
1435 if (SPI_connect() != SPI_OK_CONNECT) {
1437 elog(ERROR,
"RASTER_nMapAlgebraExpr: Could not connect to the SPI manager");
1453 char place[12] =
"$1";
1455 if (PG_ARGISNULL(exprpos[i]))
1458 expr = text_to_cstring(PG_GETARG_TEXT_P(exprpos[i]));
1461 for (j = 0, k = 1; j < argkwcount; j++) {
1473 sprintf(place,
"$%d", k);
1479 len = strlen(
"SELECT (") + strlen(expr) + strlen(
")::double precision");
1480 sql = (
char *) palloc(len + 1);
1484 elog(ERROR,
"RASTER_nMapAlgebraExpr: Could not allocate memory for expression parameter %d", exprpos[i]);
1488 memcpy(
sql,
"SELECT (", strlen(
"SELECT ("));
1489 memcpy(
sql + strlen(
"SELECT ("), expr, strlen(expr));
1490 memcpy(
sql + strlen(
"SELECT (") + strlen(expr),
")::double precision", strlen(
")::double precision"));
1499 if (argtype == NULL) {
1503 elog(ERROR,
"RASTER_nMapAlgebraExpr: Could not allocate memory for prepared plan argtypes of expression parameter %d", exprpos[i]);
1508 for (j = 0, k = 0; j < argkwcount; j++) {
1513 (strstr(argkw[j],
"[rast.x]") != NULL) ||
1514 (strstr(argkw[j],
"[rast.y]") != NULL) ||
1515 (strstr(argkw[j],
"[rast1.x]") != NULL) ||
1516 (strstr(argkw[j],
"[rast1.y]") != NULL) ||
1517 (strstr(argkw[j],
"[rast2.x]") != NULL) ||
1518 (strstr(argkw[j],
"[rast2.y]") != NULL)
1520 argtype[k] = INT4OID;
1523 argtype[k] = FLOAT8OID;
1535 elog(ERROR,
"RASTER_nMapAlgebraExpr: Could not create prepared plan of expression parameter %d", exprpos[i]);
1541 POSTGIS_RT_DEBUGF(3,
"expression parameter %d has no args, simply executing", exprpos[i]);
1542 err = SPI_execute(
sql,
TRUE, 0);
1545 if (err != SPI_OK_SELECT || SPI_tuptable == NULL || SPI_processed != 1) {
1548 elog(ERROR,
"RASTER_nMapAlgebraExpr: Could not evaluate expression parameter %d", exprpos[i]);
1553 tupdesc = SPI_tuptable->tupdesc;
1554 tuptable = SPI_tuptable;
1555 tuple = tuptable->vals[0];
1557 datum = SPI_getbinval(tuple, tupdesc, 1, &isnull);
1558 if (SPI_result == SPI_ERROR_NOATTRIBUTE) {
1559 if (SPI_tuptable) SPI_freetuptable(tuptable);
1562 elog(ERROR,
"RASTER_nMapAlgebraExpr: Could not get result of expression parameter %d", exprpos[i]);
1571 if (SPI_tuptable) SPI_freetuptable(tuptable);
1591 for (i = 0; i < numraster; i++) {
1615 if (itrset == NULL) {
1618 elog(ERROR,
"RASTER_nMapAlgebra: Could not allocate memory for iterator arguments");
1623 for (i = 0; i < numraster; i++) {
1647 elog(ERROR,
"RASTER_nMapAlgebraExpr: Could not run raster iterator function");
1650 else if (
raster == NULL) {
1656 MemoryContextSwitchTo(mainMemCtx);
1667 SET_VARSIZE(pgraster, pgraster->
size);
1668 PG_RETURN_POINTER(pgraster);
1688 assert(cutype && strlen(cutype) > 0);
1690 if (strcmp(cutype,
"LAST") == 0)
1692 else if (strcmp(cutype,
"FIRST") == 0)
1694 else if (strcmp(cutype,
"MIN") == 0)
1696 else if (strcmp(cutype,
"MAX") == 0)
1698 else if (strcmp(cutype,
"COUNT") == 0)
1700 else if (strcmp(cutype,
"SUM") == 0)
1702 else if (strcmp(cutype,
"MEAN") == 0)
1704 else if (strcmp(cutype,
"RANGE") == 0)
1731 for (i = 0; i < arg->
numband; i++) {
1755 double *
value,
int *nodata
1767 elog(ERROR,
"rtpg_union_callback: Invalid arguments passed to callback");
1783 else if (!arg->
nodata[0][0][0] && arg->
nodata[1][0][0]) {
1789 else if (arg->
nodata[0][0][0] && !arg->
nodata[1][0][0]) {
1817 else if (!arg->
nodata[0][0][0] && arg->
nodata[1][0][0])
1820 else if (arg->
nodata[0][0][0] && !arg->
nodata[1][0][0])
1846 double *
value,
int *nodata
1856 elog(ERROR,
"rtpg_union_mean_callback: Invalid arguments passed to callback");
1879 double *
value,
int *nodata
1889 elog(ERROR,
"rtpg_union_range_callback: Invalid arguments passed to callback");
1922 HeapTupleHeader tup;
1928 char *utypename = NULL;
1931 etype = ARR_ELEMTYPE(array);
1932 get_typlenbyvalalign(etype, &typlen, &typbyval, &typalign);
1937 typlen, typbyval, typalign,
1942 elog(ERROR,
"rtpg_union_unionarg_process: Invalid argument for unionarg");
1950 elog(ERROR,
"rtpg_union_unionarg_process: Could not allocate memory for band information");
1955 for (i = 0; i < n; i++) {
1964 tup = (HeapTupleHeader) DatumGetPointer(e[i]);
1966 elog(ERROR,
"rtpg_union_unionarg_process: Invalid argument for unionarg");
1971 tupv = GetAttributeByName(tup,
"nband", &isnull);
1974 elog(NOTICE,
"First argument (nband) of unionarg is NULL. Assuming nband = %d",
nband);
1977 nband = DatumGetInt32(tupv);
1980 elog(ERROR,
"rtpg_union_unionarg_process: Band number must be greater than zero (1-based)");
1985 tupv = GetAttributeByName(tup,
"uniontype", &isnull);
1987 elog(NOTICE,
"Second argument (uniontype) of unionarg is NULL. Assuming uniontype = LAST");
1991 utypename = text_to_cstring((text *) DatumGetPointer(tupv));
2012 elog(ERROR,
"rtpg_union_unionarg_process: Could not reallocate memory for band information");
2029 if (numbands <= arg->numband)
2039 elog(ERROR,
"rtpg_union_noarg: Could not reallocate memory for band information");
2045 for (; i < arg->
numband; i++) {
2053 elog(ERROR,
"rtpg_union_noarg: Could not allocate memory for working rasters");
2062 elog(ERROR,
"rtpg_union_noarg: Could not create working raster");
2075 MemoryContext aggcontext;
2076 MemoryContext oldcontext;
2086 int isempty[2] = {0};
2087 int hasband[2] = {0};
2089 double _offset[4] = {0.};
2097 char *utypename = NULL;
2101 double nodataval = 0;
2107 uint16_t _dim[2] = {0};
2114 if (!AggCheckCallContext(fcinfo, &aggcontext)) {
2115 elog(ERROR,
"RASTER_union_transfn: Cannot be called in a non-aggregate context");
2120 oldcontext = MemoryContextSwitchTo(aggcontext);
2122 if (PG_ARGISNULL(0)) {
2127 MemoryContextSwitchTo(oldcontext);
2128 elog(ERROR,
"RASTER_union_transfn: Could not allocate memory for state variable");
2144 if (!PG_ARGISNULL(1)) {
2146 pgraster = (
rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
2153 PG_FREE_IF_COPY(pgraster, 1);
2155 MemoryContextSwitchTo(oldcontext);
2156 elog(ERROR,
"RASTER_union_transfn: Could not deserialize raster");
2169 if (!PG_ARGISNULL(2)) {
2170 Oid calltype = get_fn_expr_argtype(fcinfo->flinfo, 2);
2181 utypename = text_to_cstring(PG_GETARG_TEXT_P(2));
2205 if (numband > idx) {
2223 PG_FREE_IF_COPY(pgraster, 1);
2226 MemoryContextSwitchTo(oldcontext);
2227 elog(ERROR,
"RASTER_union_transfn: Could not allocate memory for band information");
2232 for (i = idx; i < iwr->
numband; i++) {
2256 nband = PG_GETARG_INT32(2);
2262 PG_FREE_IF_COPY(pgraster, 1);
2265 MemoryContextSwitchTo(oldcontext);
2266 elog(ERROR,
"RASTER_union_transfn: Band number must be greater than zero (1-based)");
2277 PG_FREE_IF_COPY(pgraster, 1);
2280 MemoryContextSwitchTo(oldcontext);
2281 elog(ERROR,
"RASTER_union_transfn: Could not allocate memory for band information");
2302 PG_FREE_IF_COPY(pgraster, 1);
2305 MemoryContextSwitchTo(oldcontext);
2306 elog(ERROR,
"RASTER_union_transfn: Could not process unionarg");
2315 if (nargs > 3 && !PG_ARGISNULL(3)) {
2316 utypename = text_to_cstring(PG_GETARG_TEXT_P(3));
2330 for (i = 0; i < iwr->
numband; i++) {
2345 PG_FREE_IF_COPY(pgraster, 1);
2348 MemoryContextSwitchTo(oldcontext);
2349 elog(ERROR,
"RASTER_union_transfn: Could not allocate memory for working raster(s)");
2364 PG_FREE_IF_COPY(pgraster, 1);
2367 MemoryContextSwitchTo(oldcontext);
2368 elog(ERROR,
"RASTER_union_transfn: Could not create working raster");
2385 PG_FREE_IF_COPY(pgraster, 1);
2388 MemoryContextSwitchTo(oldcontext);
2389 elog(ERROR,
"RASTER_union_transfn: Could not check and balance number of bands");
2396 if (itrset == NULL) {
2401 PG_FREE_IF_COPY(pgraster, 1);
2404 MemoryContextSwitchTo(oldcontext);
2405 elog(ERROR,
"RASTER_union_transfn: Could not allocate memory for iterator arguments");
2410 for (i = 0; i < iwr->
numband; i++) {
2430 if (!isempty[0] && hasband[0])
2432 else if (!isempty[1] && hasband[1])
2439 if (_band != NULL) {
2477 itrset[0].
nband = 0;
2493 if (!isempty[0] && !isempty[1] && hasband[0] && hasband[1]) {
2507 PG_FREE_IF_COPY(pgraster, 1);
2510 MemoryContextSwitchTo(oldcontext);
2511 elog(ERROR,
"RASTER_union_transfn: Could not create internal raster");
2515 _offset[0], _offset[1], _offset[2], _offset[3]);
2522 double igt[6] = {0};
2538 hasnodata, nodataval,
2547 PG_FREE_IF_COPY(pgraster, 1);
2550 MemoryContextSwitchTo(oldcontext);
2551 elog(ERROR,
"RASTER_union_transfn: Could not add new band to internal raster");
2559 for (
y = 0;
y < _dim[1];
y++) {
2560 POSTGIS_RT_DEBUGF(4,
"Getting pixel line of working raster at (x, y, length) = (0, %d, %d)",
y, _dim[0]);
2574 PG_FREE_IF_COPY(pgraster, 1);
2577 MemoryContextSwitchTo(oldcontext);
2578 elog(ERROR,
"RASTER_union_transfn: Could not get pixel line from band of working raster");
2582 POSTGIS_RT_DEBUGF(4,
"Setting pixel line at (x, y, length) = (%d, %d, %d)", (
int) _offset[0], (
int) _offset[1] +
y, nvals);
2585 (
int) _offset[0], (
int) _offset[1] +
y,
2595 PG_FREE_IF_COPY(pgraster, 1);
2598 MemoryContextSwitchTo(oldcontext);
2599 elog(ERROR,
"RASTER_union_transfn: Could not set pixel line to band of internal raster");
2615 hasnodata, nodataval,
2632 PG_FREE_IF_COPY(pgraster, 1);
2635 MemoryContextSwitchTo(oldcontext);
2636 elog(ERROR,
"RASTER_union_transfn: Could not run raster iterator function");
2644 for (
y = 0;
y < _dim[1];
y++) {
2645 POSTGIS_RT_DEBUGF(4,
"Getting pixel line of iterator raster at (x, y, length) = (0, %d, %d)",
y, _dim[0]);
2661 PG_FREE_IF_COPY(pgraster, 1);
2664 MemoryContextSwitchTo(oldcontext);
2665 elog(ERROR,
"RASTER_union_transfn: Could not get pixel line from band of working raster");
2669 POSTGIS_RT_DEBUGF(4,
"Setting pixel line at (x, y, length) = (%d, %d, %d)", (
int) _offset[2], (
int) _offset[3] +
y, nvals);
2672 (
int) _offset[2], (
int) _offset[3] +
y,
2684 PG_FREE_IF_COPY(pgraster, 1);
2687 MemoryContextSwitchTo(oldcontext);
2688 elog(ERROR,
"RASTER_union_transfn: Could not set pixel line to band of internal raster");
2708 hasnodata, nodataval,
2722 PG_FREE_IF_COPY(pgraster, 1);
2725 MemoryContextSwitchTo(oldcontext);
2726 elog(ERROR,
"RASTER_union_transfn: Could not run raster iterator function");
2745 PG_FREE_IF_COPY(pgraster, 1);
2749 MemoryContextSwitchTo(oldcontext);
2753 PG_RETURN_POINTER(iwr);
2773 double nodataval = 0;
2778 if (!AggCheckCallContext(fcinfo, NULL)) {
2779 elog(ERROR,
"RASTER_union_finalfn: Cannot be called in a non-aggregate context");
2784 if (PG_ARGISNULL(0))
2791 if (itrset == NULL) {
2793 elog(ERROR,
"RASTER_union_finalfn: Could not allocate memory for iterator arguments");
2797 for (i = 0; i < iwr->
numband; i++) {
2812 itrset[0].
nband = 0;
2814 itrset[1].
nband = 0;
2822 hasnodata, nodataval,
2835 hasnodata, nodataval,
2849 elog(ERROR,
"RASTER_union_finalfn: Could not run raster iterator function");
2855 if (_raster == NULL)
2861 uint32_t bandNums[1] = {0};
2863 status = (_rtn == NULL) ? -1 : 0;
2888 elog(ERROR,
"RASTER_union_finalfn: Could not add band to final raster");
2901 if (!_rtn) PG_RETURN_NULL();
2911 SET_VARSIZE(pgraster, pgraster->
size);
2912 PG_RETURN_POINTER(pgraster);
2941 elog(ERROR,
"rtpg_clip_arg_init: Could not allocate memory for function arguments");
2955 if (arg->
band != NULL)
2960 if (arg->
mask != NULL)
3000 double offset[4] = {0.};
3010 int mask_height = 0;
3020 if (PG_ARGISNULL(0) || PG_ARGISNULL(2))
3026 elog(ERROR,
"RASTER_clip: Could not initialize argument structure");
3031 pgraster = (
rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
3035 if (arg->
raster == NULL) {
3037 PG_FREE_IF_COPY(pgraster, 0);
3038 elog(ERROR,
"RASTER_clip: Could not deserialize raster");
3044 elog(NOTICE,
"Input raster is empty or has no bands. Returning empty raster");
3047 PG_FREE_IF_COPY(pgraster, 0);
3051 elog(ERROR,
"RASTER_clip: Could not create empty raster");
3060 SET_VARSIZE(pgrtn, pgrtn->
size);
3061 PG_RETURN_POINTER(pgrtn);
3069 gser = PG_GETARG_GSERIALIZED_P(2);
3081 elog(NOTICE,
"Geometry provided does not have the same SRID as the raster. Returning NULL");
3084 PG_FREE_IF_COPY(pgraster, 0);
3086 PG_FREE_IF_COPY(gser, 2);
3092 if (!PG_ARGISNULL(4) && !PG_GETARG_BOOL(4))
3099 PG_FREE_IF_COPY(pgraster, 0);
3101 PG_FREE_IF_COPY(gser, 2);
3103 elog(ERROR,
"RASTER_clip: Could not get convex hull of raster");
3110 PG_FREE_IF_COPY(gser, 2);
3115 elog(NOTICE,
"The input raster and input geometry do not intersect. Returning empty raster");
3118 PG_FREE_IF_COPY(pgraster, 0);
3123 elog(ERROR,
"RASTER_clip: Could not create empty raster");
3132 SET_VARSIZE(pgrtn, pgrtn->
size);
3133 PG_RETURN_POINTER(pgrtn);
3137 if (!PG_ARGISNULL(1)) {
3138 array = PG_GETARG_ARRAYTYPE_P(1);
3139 etype = ARR_ELEMTYPE(array);
3140 get_typlenbyvalalign(etype, &typlen, &typbyval, &typalign);
3148 PG_FREE_IF_COPY(pgraster, 0);
3150 elog(ERROR,
"RASTER_clip: Invalid data type for band indexes");
3157 typlen, typbyval, typalign,
3162 if (arg->
band == NULL) {
3164 PG_FREE_IF_COPY(pgraster, 0);
3166 elog(ERROR,
"RASTER_clip: Could not allocate memory for band arguments");
3170 for (i = 0, j = 0; i < arg->
numbands; i++) {
3171 if (nulls[i])
continue;
3175 arg->
band[j].
nband = DatumGetInt16(e[i]) - 1;
3178 arg->
band[j].
nband = DatumGetInt32(e[i]) - 1;
3185 if (j < arg->numbands) {
3187 if (arg->
band == NULL) {
3189 PG_FREE_IF_COPY(pgraster, 0);
3191 elog(ERROR,
"RASTER_clip: Could not reallocate memory for band arguments");
3199 for (i = 0; i < arg->
numbands; i++) {
3201 elog(NOTICE,
"Band at index %d not found in raster", arg->
band[i].
nband + 1);
3203 PG_FREE_IF_COPY(pgraster, 0);
3218 if (arg->
band == NULL) {
3221 PG_FREE_IF_COPY(pgraster, 0);
3224 elog(ERROR,
"RASTER_clip: Could not allocate memory for band arguments");
3228 for (i = 0; i < arg->
numbands; i++) {
3237 if (!PG_ARGISNULL(3)) {
3238 array = PG_GETARG_ARRAYTYPE_P(3);
3239 etype = ARR_ELEMTYPE(array);
3240 get_typlenbyvalalign(etype, &typlen, &typbyval, &typalign);
3248 PG_FREE_IF_COPY(pgraster, 0);
3250 elog(ERROR,
"RASTER_clip: Invalid data type for NODATA values");
3257 typlen, typbyval, typalign,
3262 for (i = 0, j = 0; i < arg->
numbands; i++, j++) {
3310 if (arg->
mask == NULL) {
3312 PG_FREE_IF_COPY(pgraster, 0);
3313 elog(ERROR,
"RASTER_clip: Could not rasterize intersection geometry");
3325 PG_FREE_IF_COPY(pgraster, 0);
3326 elog(ERROR,
"RASTER_clip: Could not compute extent of rasters");
3335 for (i = 0; i < arg->
numbands; i++) {
3356 PG_FREE_IF_COPY(pgraster, 0);
3357 elog(ERROR,
"RASTER_clip: Could not add new band in output raster");
3371 for (
y = 0;
y < height;
y++) {
3372 for (
x = 0;
x < width;
x++) {
3373 mask_x =
x - (int)offset[2];
3374 mask_y =
y - (int)offset[3];
3378 mask_x < mask_width &&
3380 mask_y < mask_height
3387 PG_FREE_IF_COPY(pgraster, 0);
3388 elog(ERROR,
"RASTER_clip: Could not get pixel value");
3396 input_x =
x - (int)offset[0];
3397 input_y =
y - (int)offset[1];
3401 PG_FREE_IF_COPY(pgraster, 0);
3402 elog(ERROR,
"RASTER_clip: Could not get pixel value");
3412 PG_FREE_IF_COPY(pgraster, 0);
3413 elog(ERROR,
"RASTER_clip: Could not set pixel value");
3421 PG_FREE_IF_COPY(pgraster, 0);
3431 SET_VARSIZE(pgrtn, pgrtn->
size);
3432 PG_RETURN_POINTER(pgrtn);
3445 uint32_t numBands = 0;
3465 HeapTupleHeader tup;
3470 text *exprtext = NULL;
3475 char *pixeltype = NULL;
3476 text *pixeltypetext = NULL;
3478 double nodataval = 0;
3479 bool hasnodata =
FALSE;
3481 char **comma_set = NULL;
3482 uint32_t comma_n = 0;
3483 char **colon_set = NULL;
3484 uint32_t colon_n = 0;
3485 char **dash_set = NULL;
3486 uint32_t dash_n = 0;
3491 if (PG_ARGISNULL(0))
3493 pgraster = (
rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
3498 PG_FREE_IF_COPY(pgraster, 0);
3499 elog(ERROR,
"RASTER_reclass: Could not deserialize raster");
3503 POSTGIS_RT_DEBUGF(3,
"RASTER_reclass: %d possible bands to be reclassified", numBands);
3507 array = PG_GETARG_ARRAYTYPE_P(1);
3508 etype = ARR_ELEMTYPE(array);
3509 get_typlenbyvalalign(etype, &typlen, &typbyval, &typalign);
3511 deconstruct_array(array, etype, typlen, typbyval, typalign, &e,
3515 elog(NOTICE,
"Invalid argument for reclassargset. Returning original raster");
3519 PG_FREE_IF_COPY(pgraster, 0);
3523 SET_VARSIZE(pgrtn, pgrtn->
size);
3524 PG_RETURN_POINTER(pgrtn);
3532 for (i = 0; i < n; i++) {
3533 if (nulls[i])
continue;
3536 tup = (HeapTupleHeader) DatumGetPointer(e[i]);
3538 elog(NOTICE,
"Invalid argument for reclassargset. Returning original raster");
3542 PG_FREE_IF_COPY(pgraster, 0);
3546 SET_VARSIZE(pgrtn, pgrtn->
size);
3547 PG_RETURN_POINTER(pgrtn);
3551 tupv = GetAttributeByName(tup,
"nband", &isnull);
3553 elog(NOTICE,
"Invalid argument for reclassargset. Missing value of nband for reclassarg of index %d . Returning original raster", i);
3557 PG_FREE_IF_COPY(pgraster, 0);
3561 SET_VARSIZE(pgrtn, pgrtn->
size);
3562 PG_RETURN_POINTER(pgrtn);
3564 nband = DatumGetInt32(tupv);
3568 if (nband < 1 || nband > numBands) {
3569 elog(NOTICE,
"Invalid argument for reclassargset. Invalid band index (must use 1-based) for reclassarg of index %d . Returning original raster", i);
3573 PG_FREE_IF_COPY(pgraster, 0);
3577 SET_VARSIZE(pgrtn, pgrtn->
size);
3578 PG_RETURN_POINTER(pgrtn);
3582 tupv = GetAttributeByName(tup,
"reclassexpr", &isnull);
3584 elog(NOTICE,
"Invalid argument for reclassargset. Missing value of reclassexpr for reclassarg of index %d . Returning original raster", i);
3588 PG_FREE_IF_COPY(pgraster, 0);
3592 SET_VARSIZE(pgrtn, pgrtn->
size);
3593 PG_RETURN_POINTER(pgrtn);
3595 exprtext = (text *) DatumGetPointer(tupv);
3596 if (NULL == exprtext) {
3597 elog(NOTICE,
"Invalid argument for reclassargset. Missing value of reclassexpr for reclassarg of index %d . Returning original raster", i);
3601 PG_FREE_IF_COPY(pgraster, 0);
3605 SET_VARSIZE(pgrtn, pgrtn->
size);
3606 PG_RETURN_POINTER(pgrtn);
3608 expr = text_to_cstring(exprtext);
3617 elog(NOTICE,
"Invalid argument for reclassargset. Invalid expression of reclassexpr for reclassarg of index %d . Returning original raster", i);
3621 PG_FREE_IF_COPY(pgraster, 0);
3625 SET_VARSIZE(pgrtn, pgrtn->
size);
3626 PG_RETURN_POINTER(pgrtn);
3633 for (a = 0, j = 0; a < comma_n; a++) {
3639 elog(NOTICE,
"Invalid argument for reclassargset. Invalid expression of reclassexpr for reclassarg of index %d . Returning original raster", i);
3640 for (k = 0; k < j; k++) pfree(exprset[k]);
3645 PG_FREE_IF_COPY(pgraster, 0);
3649 SET_VARSIZE(pgrtn, pgrtn->
size);
3650 PG_RETURN_POINTER(pgrtn);
3656 for (b = 0; b < colon_n; b++) {
3661 if (dash_n < 1 || dash_n > 3) {
3662 elog(NOTICE,
"Invalid argument for reclassargset. Invalid expression of reclassexpr for reclassarg of index %d . Returning original raster", i);
3663 for (k = 0; k < j; k++) pfree(exprset[k]);
3668 PG_FREE_IF_COPY(pgraster, 0);
3672 SET_VARSIZE(pgrtn, pgrtn->
size);
3673 PG_RETURN_POINTER(pgrtn);
3676 for (c = 0; c < dash_n; c++) {
3680 strlen(dash_set[c]) == 1 && (
3681 strchr(dash_set[c],
'(') != NULL ||
3682 strchr(dash_set[c],
'[') != NULL ||
3683 strchr(dash_set[c],
')') != NULL ||
3684 strchr(dash_set[c],
']') != NULL
3688 junk = palloc(
sizeof(
char) * (strlen(dash_set[c + 1]) + 2));
3690 for (k = 0; k <= j; k++) pfree(exprset[k]);
3693 PG_FREE_IF_COPY(pgraster, 0);
3695 elog(ERROR,
"RASTER_reclass: Could not allocate memory");
3699 sprintf(junk,
"%s%s", dash_set[c], dash_set[c + 1]);
3701 dash_set[c] = repalloc(dash_set[c],
sizeof(
char) * (strlen(junk) + 1));
3702 strcpy(dash_set[c], junk);
3706 for (dash_it = 1; dash_it < dash_n; dash_it++) {
3707 dash_set[dash_it - 1] = repalloc(dash_set[dash_it - 1], (strlen(dash_set[dash_it]) + 1) *
sizeof(
char));
3708 strcpy(dash_set[dash_it - 1], dash_set[dash_it]);
3712 pfree(dash_set[dash_n]);
3713 dash_set = repalloc(dash_set,
sizeof(
char *) * dash_n);
3717 if (c < 1 && dash_n > 2) {
3718 elog(NOTICE,
"Invalid argument for reclassargset. Invalid expression of reclassexpr for reclassarg of index %d . Returning original raster", i);
3719 for (k = 0; k < j; k++) pfree(exprset[k]);
3724 PG_FREE_IF_COPY(pgraster, 0);
3728 SET_VARSIZE(pgrtn, pgrtn->
size);
3729 PG_RETURN_POINTER(pgrtn);
3740 strchr(dash_set[c],
')') != NULL ||
3741 strchr(dash_set[c],
']') != NULL
3746 else if (strchr(dash_set[c],
'(') != NULL){
3756 strrchr(dash_set[c],
'(') != NULL ||
3757 strrchr(dash_set[c],
'[') != NULL
3762 else if (strrchr(dash_set[c],
']') != NULL) {
3770 POSTGIS_RT_DEBUGF(4,
"RASTER_reclass: exc_val %d inc_val %d", exc_val, inc_val);
3778 val = strtod(dash_set[c], &junk);
3779 if (errno != 0 || dash_set[c] == junk) {
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);
3797 junk = strstr(colon_set[b], dash_set[c]);
3802 "(colon_set[%d], dash_set[%d], junk) = (%s, %s, %s)",
3803 b, c, colon_set[b], dash_set[c], junk
3806 if (junk != colon_set[b]) {
3808 if (*(junk - 1) ==
'-') {
3811 ((junk - 1) == colon_set[b]) ||
3812 (*(junk - 2) ==
'-') ||
3813 (*(junk - 2) ==
'[') ||
3814 (*(junk - 2) ==
'(')
3834 exprset[j]->
src.
min = val;
3840 exprset[j]->
src.
max = val;
3850 exprset[j]->
dst.
min = val;
3853 exprset[j]->
dst.
max = val;
3861 , exprset[j]->src.min
3871 tupv = GetAttributeByName(tup,
"pixeltype", &isnull);
3873 elog(NOTICE,
"Invalid argument for reclassargset. Missing value of pixeltype for reclassarg of index %d . Returning original raster", i);
3877 PG_FREE_IF_COPY(pgraster, 0);
3881 SET_VARSIZE(pgrtn, pgrtn->
size);
3882 PG_RETURN_POINTER(pgrtn);
3884 pixeltypetext = (text *) DatumGetPointer(tupv);
3885 if (NULL == pixeltypetext) {
3886 elog(NOTICE,
"Invalid argument for reclassargset. Missing value of pixeltype for reclassarg of index %d . Returning original raster", i);
3890 PG_FREE_IF_COPY(pgraster, 0);
3894 SET_VARSIZE(pgrtn, pgrtn->
size);
3895 PG_RETURN_POINTER(pgrtn);
3897 pixeltype = text_to_cstring(pixeltypetext);
3902 tupv = GetAttributeByName(tup,
"nodataval", &isnull);
3908 nodataval = DatumGetFloat8(tupv);
3917 elog(NOTICE,
"Could not find raster band of index %d. Returning original raster",
nband);
3918 for (k = 0; k < j; k++) pfree(exprset[k]);
3923 PG_FREE_IF_COPY(pgraster, 0);
3927 SET_VARSIZE(pgrtn, pgrtn->
size);
3928 PG_RETURN_POINTER(pgrtn);
3932 for (k = 0; k < j; k++) pfree(exprset[k]);
3936 PG_FREE_IF_COPY(pgraster, 0);
3938 elog(ERROR,
"RASTER_reclass: Could not reclassify raster band of index %d",
nband);
3944 for (k = 0; k < j; k++) pfree(exprset[k]);
3949 PG_FREE_IF_COPY(pgraster, 0);
3951 elog(ERROR,
"RASTER_reclass: Could not replace raster band of index %d with reclassified band",
nband);
3959 for (k = 0; k < j; k++) pfree(exprset[k]);
3965 PG_FREE_IF_COPY(pgraster, 0);
3971 SET_VARSIZE(pgrtn, pgrtn->
size);
3972 PG_RETURN_POINTER(pgrtn);
4001 elog(ERROR,
"rtpg_colormap_arg: Could not allocate memory for function arguments");
4012 elog(ERROR,
"rtpg_colormap_arg: Could not allocate memory for function arguments");
4045 for (i = 0; i < arg->
nentry; i++) {
4046 if (arg->
entry[i] != NULL)
4047 pfree(arg->
entry[i]);
4053 for (i = 0; i < arg->
nelement; i++)
4073 if (PG_ARGISNULL(0))
4079 elog(ERROR,
"RASTER_colorMap: Could not initialize argument structure");
4084 pgraster = (
rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
4090 PG_FREE_IF_COPY(pgraster, 0);
4091 elog(ERROR,
"RASTER_colorMap: Could not deserialize raster");
4096 if (!PG_ARGISNULL(1))
4097 arg->
nband = PG_GETARG_INT32(1);
4102 elog(NOTICE,
"Raster does not have band at index %d. Returning empty raster", arg->
nband);
4107 PG_FREE_IF_COPY(pgraster, 0);
4108 elog(ERROR,
"RASTER_colorMap: Could not create empty raster");
4113 PG_FREE_IF_COPY(pgraster, 0);
4117 if (pgraster == NULL)
4120 SET_VARSIZE(pgraster, ((
rt_pgraster*) pgraster)->size);
4121 PG_RETURN_POINTER(pgraster);
4126 if (arg->
band == NULL) {
4129 PG_FREE_IF_COPY(pgraster, 0);
4130 elog(ERROR,
"RASTER_colorMap: Could not get band at index %d",
nband);
4135 if (!PG_ARGISNULL(3)) {
4136 char *method = NULL;
4137 char *tmp = text_to_cstring(PG_GETARG_TEXT_P(3));
4144 if (strcmp(method,
"INTERPOLATE") == 0)
4146 else if (strcmp(method,
"EXACT") == 0)
4148 else if (strcmp(method,
"NEAREST") == 0)
4151 elog(NOTICE,
"Unknown value provided for method. Defaulting to INTERPOLATE");
4161 if (PG_ARGISNULL(2)) {
4163 PG_FREE_IF_COPY(pgraster, 0);
4164 elog(ERROR,
"RASTER_colorMap: Value must be provided for colormap");
4169 char *colormap = text_to_cstring(PG_GETARG_TEXT_P(2));
4178 if (!strlen(colormap)) {
4180 PG_FREE_IF_COPY(pgraster, 0);
4181 elog(ERROR,
"RASTER_colorMap: Value must be provided for colormap");
4189 PG_FREE_IF_COPY(pgraster, 0);
4190 elog(ERROR,
"RASTER_colorMap: Could not process the value provided for colormap");
4198 PG_FREE_IF_COPY(pgraster, 0);
4199 elog(ERROR,
"RASTER_colorMap: Could not allocate memory for colormap entries");
4205 for (i = 0; i < arg->
nentry; i++) {
4219 if (!strlen(_entry)) {
4229 PG_FREE_IF_COPY(pgraster, 0);
4230 elog(ERROR,
"RASTER_colorMap: Could not process colormap entry %d", i + 1);
4234 elog(NOTICE,
"More than five elements in colormap entry %d. Using at most five elements", i + 1);
4243 for (j = 0; j < arg->
nelement; j++) {
4252 char *percent = NULL;
4256 strcmp(_element,
"NV") == 0 ||
4257 strcmp(_element,
"NULL") == 0 ||
4258 strcmp(_element,
"NODATA") == 0
4263 elog(NOTICE,
"More than one NODATA entry found. Using only the first one");
4272 else if ((percent = strchr(_element,
'%')) != NULL) {
4284 PG_FREE_IF_COPY(pgraster, 0);
4285 elog(ERROR,
"RASTER_colorMap: Could not get band's summary stats to process percentages");
4291 tmp = palloc(
sizeof(
char) * (percent - _element + 1));
4295 PG_FREE_IF_COPY(pgraster, 0);
4296 elog(ERROR,
"RASTER_colorMap: Could not allocate memory for value of percentage");
4300 memcpy(tmp, _element, percent - _element);
4301 tmp[percent - _element] =
'\0';
4306 value = strtod(tmp, NULL);
4308 if (errno != 0 || _element == junk) {
4311 PG_FREE_IF_COPY(pgraster, 0);
4312 elog(ERROR,
"RASTER_colorMap: Could not process percent string to value");
4318 elog(NOTICE,
"Percentage values cannot be less than zero. Defaulting to zero");
4321 else if (
value > 100.) {
4322 elog(NOTICE,
"Percentage values cannot be greater than 100. Defaulting to 100");
4334 if (errno != 0 || _element == junk) {
4337 PG_FREE_IF_COPY(pgraster, 0);
4338 elog(ERROR,
"RASTER_colorMap: Could not process string to value");
4349 value = (int) strtod(_element, &junk);
4350 if (errno != 0 || _element == junk) {
4353 PG_FREE_IF_COPY(pgraster, 0);
4354 elog(ERROR,
"RASTER_colorMap: Could not process string to value");
4359 elog(NOTICE,
"RGBA value cannot be greater than 255. Defaulting to 255");
4362 else if (
value < 0) {
4363 elog(NOTICE,
"RGBA value cannot be less than zero. Defaulting to zero");
4372 POSTGIS_RT_DEBUGF(4,
"colormap->entry[%d] (isnodata, value, R, G, B, A) = (%d, %f, %d, %d, %d, %d)",
4393 PG_FREE_IF_COPY(pgraster, 0);
4394 elog(ERROR,
"RASTER_colorMap: Could not create new raster with applied colormap");
4399 PG_FREE_IF_COPY(pgraster, 0);
4405 if (pgraster == NULL)
4408 SET_VARSIZE(pgraster, ((
rt_pgraster*) pgraster)->size);
4409 PG_RETURN_POINTER(pgraster);
4421 int x,
y,
nband, width, height;
4423 double newnodatavalue = 0.0;
4424 double newinitialvalue = 0.0;
4425 double newval = 0.0;
4426 char *newexpr = NULL;
4427 char *initexpr = NULL;
4428 char *expression = NULL;
4429 int hasnodataval = 0;
4430 double nodataval = 0.;
4432 int skipcomputation = 0;
4434 const int argkwcount = 3;
4435 enum KEYWORDS { kVAL=0, kX=1, kY=2 };
4436 char *argkw[] = {
"[rast]",
"[rast.x]",
"[rast.y]"};
4437 Oid argkwtypes[] = { FLOAT8OID, INT4OID, INT4OID };
4439 Oid argtype[] = { FLOAT8OID, INT4OID, INT4OID };
4440 uint8_t argpos[3] = {0};
4445 SPIPlanPtr spi_plan = NULL;
4446 SPITupleTable * tuptable = NULL;
4448 char * strFromText = NULL;
4449 Datum *values = NULL;
4450 Datum datum = (Datum)NULL;
4452 bool isnull =
FALSE;
4459 if (PG_ARGISNULL(0)) {
4460 elog(NOTICE,
"Raster is NULL. Returning NULL");
4466 pgraster = (
rt_pgraster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
4469 PG_FREE_IF_COPY(pgraster, 0);
4470 elog(ERROR,
"RASTER_mapAlgebraExpr: Could not deserialize raster");
4476 if (PG_ARGISNULL(1))
4479 nband = PG_GETARG_INT32(1);
4485 POSTGIS_RT_DEBUG(3,
"RASTER_mapAlgebraExpr: Creating new empty raster...");
4496 if ( NULL == newrast ) {
4497 PG_FREE_IF_COPY(pgraster, 0);
4498 elog(ERROR,
"RASTER_mapAlgebraExpr: Could not create a new raster");
4523 elog(NOTICE,
"Raster is empty. Returning an empty raster");
4525 PG_FREE_IF_COPY(pgraster, 0);
4529 if (NULL == pgrtn) {
4531 elog(ERROR,
"RASTER_mapAlgebraExpr: Could not serialize raster");
4535 SET_VARSIZE(pgrtn, pgrtn->
size);
4536 PG_RETURN_POINTER(pgrtn);
4547 elog(NOTICE,
"Raster does not have the required band. Returning a raster "
4550 PG_FREE_IF_COPY(pgraster, 0);
4554 if (NULL == pgrtn) {
4555 elog(ERROR,
"RASTER_mapAlgebraExpr: Could not serialize raster");
4559 SET_VARSIZE(pgrtn, pgrtn->
size);
4560 PG_RETURN_POINTER(pgrtn);
4565 if ( NULL ==
band ) {
4566 elog(NOTICE,
"Could not get the required band. Returning a raster "
4569 PG_FREE_IF_COPY(pgraster, 0);
4573 if (NULL == pgrtn) {
4574 elog(ERROR,
"RASTER_mapAlgebraExpr: Could not serialize raster");
4578 SET_VARSIZE(pgrtn, pgrtn->
size);
4579 PG_RETURN_POINTER(pgrtn);
4585 POSTGIS_RT_DEBUG(3,
"RASTER_mapAlgebraExpr: Getting NODATA value for band...");
4603 newinitialvalue = newnodatavalue;
4610 if (PG_ARGISNULL(2)) {
4615 strFromText = text_to_cstring(PG_GETARG_TEXT_P(2));
4618 if (newpixeltype ==
PT_END)
4622 if (newpixeltype ==
PT_END) {
4623 PG_FREE_IF_COPY(pgraster, 0);
4624 elog(ERROR,
"RASTER_mapAlgebraExpr: Invalid pixeltype");
4633 if (!PG_ARGISNULL(3)) {
4634 expression = text_to_cstring(PG_GETARG_TEXT_P(3));
4635 len = strlen(
"SELECT (") + strlen(expression) + strlen(
")::double precision");
4636 initexpr = (
char *)palloc(len + 1);
4638 memcpy(initexpr,
"SELECT (", strlen(
"SELECT ("));
4639 memcpy(initexpr + strlen(
"SELECT ("), expression, strlen(expression));
4640 memcpy(initexpr + strlen(
"SELECT (") + strlen(expression),
")::double precision", strlen(
")::double precision"));
4641 initexpr[len] =
'\0';
4659 if (!PG_ARGISNULL(4)) {
4661 nodataval = PG_GETARG_FLOAT8(4);
4662 newinitialvalue = nodataval;
4679 POSTGIS_RT_DEBUG(3,
"RASTER_mapAlgebraExpr: Band is a nodata band, returning "
4680 "a raster filled with nodata");
4683 newinitialvalue,
TRUE, newnodatavalue, 0);
4689 PG_FREE_IF_COPY(pgraster, 0);
4694 if (NULL == pgrtn) {
4695 elog(ERROR,
"RASTER_mapAlgebraExpr: Could not serialize raster");
4699 SET_VARSIZE(pgrtn, pgrtn->
size);
4700 PG_RETURN_POINTER(pgrtn);
4708 if (initexpr != NULL && ( !strcmp(initexpr,
"SELECT [rast]") || !strcmp(initexpr,
"SELECT [rast.val]") ) && !hasnodataval) {
4711 "Returning raster with band %d from original raster",
nband);
4723 PG_FREE_IF_COPY(pgraster, 0);
4728 if (NULL == pgrtn) {
4729 elog(ERROR,
"RASTER_mapAlgebraExpr: Could not serialize raster");
4733 SET_VARSIZE(pgrtn, pgrtn->
size);
4734 PG_RETURN_POINTER(pgrtn);
4741 if (initexpr != NULL && strstr(initexpr,
"[rast") == NULL) {
4742 ret = SPI_connect();
4743 if (ret != SPI_OK_CONNECT) {
4744 PG_FREE_IF_COPY(pgraster, 0);
4745 elog(ERROR,
"RASTER_mapAlgebraExpr: Could not connect to the SPI manager");
4750 ret = SPI_execute(initexpr,
FALSE, 0);
4752 if (ret != SPI_OK_SELECT || SPI_tuptable == NULL || SPI_processed != 1) {
4756 SPI_freetuptable(tuptable);
4757 PG_FREE_IF_COPY(pgraster, 0);
4760 elog(ERROR,
"RASTER_mapAlgebraExpr: Invalid construction for expression");
4764 tupdesc = SPI_tuptable->tupdesc;
4765 tuptable = SPI_tuptable;
4767 tuple = tuptable->vals[0];
4768 newexpr = SPI_getvalue(tuple, tupdesc, 1);
4770 POSTGIS_RT_DEBUG(3,
"Constant expression evaluated to NULL, keeping initvalue");
4771 newval = newinitialvalue;
4773 newval = atof(newexpr);
4776 SPI_freetuptable(tuptable);
4783 skipcomputation = 1;
4789 if (!hasnodataval) {
4790 newinitialvalue = newval;
4791 skipcomputation = 2;
4795 else if (
FLT_NEQ(newval, newinitialvalue)) {
4796 skipcomputation = 2;
4805 newinitialvalue,
TRUE, newnodatavalue, 0);
4812 if (expression == NULL || skipcomputation == 2) {
4818 PG_FREE_IF_COPY(pgraster, 0);
4823 if (NULL == pgrtn) {
4824 elog(ERROR,
"RASTER_mapAlgebraExpr: Could not serialize raster");
4828 SET_VARSIZE(pgrtn, pgrtn->
size);
4829 PG_RETURN_POINTER(pgrtn);
4832 RASTER_DEBUG(3,
"RASTER_mapAlgebraExpr: Creating new raster band...");
4836 if ( NULL == newband ) {
4837 elog(NOTICE,
"Could not modify band for new raster. Returning new "
4838 "raster with the original band");
4843 PG_FREE_IF_COPY(pgraster, 0);
4848 if (NULL == pgrtn) {
4849 elog(ERROR,
"RASTER_mapAlgebraExpr: Could not serialize raster");
4853 SET_VARSIZE(pgrtn, pgrtn->
size);
4854 PG_RETURN_POINTER(pgrtn);
4860 if (initexpr != NULL) {
4863 pfree(initexpr); initexpr=newexpr;
4865 sprintf(place,
"$1");
4866 for (i = 0, j = 1; i < argkwcount; i++) {
4869 pfree(initexpr); initexpr=newexpr;
4871 argtype[argcount] = argkwtypes[i];
4875 sprintf(place,
"$%d", j);
4885 values = (Datum *) palloc(
sizeof(Datum) * argcount);
4886 if (values == NULL) {
4891 PG_FREE_IF_COPY(pgraster, 0);
4894 elog(ERROR,
"RASTER_mapAlgebraExpr: Could not allocate memory for value parameters of prepared statement");
4899 nulls = (
char *)palloc(argcount);
4900 if (nulls == NULL) {
4905 PG_FREE_IF_COPY(pgraster, 0);
4908 elog(ERROR,
"RASTER_mapAlgebraExpr: Could not allocate memory for null parameters of prepared statement");
4913 ret = SPI_connect();
4914 if (ret != SPI_OK_CONNECT) {
4919 PG_FREE_IF_COPY(pgraster, 0);
4922 elog(ERROR,
"RASTER_mapAlgebraExpr: Could not connect to the SPI manager");
4927 spi_plan = SPI_prepare(initexpr, argcount, argtype);
4929 if (spi_plan == NULL) {
4932 PG_FREE_IF_COPY(pgraster, 0);
4939 elog(ERROR,
"RASTER_mapAlgebraExpr: Could not prepare expression");
4944 for (
x = 0;
x < width;
x++) {
4945 for(
y = 0;
y < height;
y++) {
4953 if (skipcomputation == 0) {
4954 if (initexpr != NULL) {
4956 memset(nulls,
'n', argcount);
4958 for (i = 0; i < argkwcount; i++) {
4960 if (idx < 1)
continue;
4965 values[idx] = Int32GetDatum(
x+1);
4970 values[idx] = Int32GetDatum(
y+1);
4973 else if (i == kVAL ) {
4974 values[idx] = Float8GetDatum(
r);
4980 ret = SPI_execute_plan(spi_plan, values, nulls,
FALSE, 0);
4981 if (ret != SPI_OK_SELECT || SPI_tuptable == NULL ||
4982 SPI_processed != 1) {
4984 SPI_freetuptable(tuptable);
4986 SPI_freeplan(spi_plan);
4994 PG_FREE_IF_COPY(pgraster, 0);
4997 elog(ERROR,
"RASTER_mapAlgebraExpr: Error executing prepared plan");
5002 tupdesc = SPI_tuptable->tupdesc;
5003 tuptable = SPI_tuptable;
5005 tuple = tuptable->vals[0];
5006 datum = SPI_getbinval(tuple, tupdesc, 1, &isnull);
5007 if ( SPI_result == SPI_ERROR_NOATTRIBUTE ) {
5008 POSTGIS_RT_DEBUGF(3,
"Expression for pixel %d,%d (value %g) errored, skip setting",
x+1,
y+1,
r);
5009 newval = newinitialvalue;
5011 else if ( isnull ) {
5012 POSTGIS_RT_DEBUGF(3,
"Expression for pixel %d,%d (value %g) evaluated to NULL, skip setting",
x+1,
y+1,
r);
5013 newval = newinitialvalue;
5015 newval = DatumGetFloat8(datum);
5018 SPI_freetuptable(tuptable);
5022 newval = newinitialvalue;
5035 if (initexpr != NULL) {
5036 SPI_freeplan(spi_plan);
5050 POSTGIS_RT_DEBUG(3,
"RASTER_mapAlgebraExpr: raster modified, serializing it.");
5054 PG_FREE_IF_COPY(pgraster, 0);
5061 SET_VARSIZE(pgrtn, pgrtn->
size);
5069 PG_RETURN_POINTER(pgrtn);
5084 int x,
y,
nband, width, height;
5086 double newnodatavalue = 0.0;
5087 double newinitialvalue = 0.0;
5088 double newval = 0.0;
5093 LOCAL_FCINFO(cbdata, FUNC_MAX_ARGS);
5096 char * strFromText = NULL;
5102 if (PG_ARGISNULL(0)) {
5103 elog(WARNING,
"Raster is NULL. Returning NULL");
5109 pgraster = (
rt_pgraster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
5112 PG_FREE_IF_COPY(pgraster, 0);
5113 elog(ERROR,
"RASTER_mapAlgebraFct: Could not deserialize raster");
5121 if (PG_ARGISNULL(1))
5124 nband = PG_GETARG_INT32(1);
5140 if ( NULL == newrast ) {
5143 PG_FREE_IF_COPY(pgraster, 0);
5145 elog(ERROR,
"RASTER_mapAlgebraFct: Could not create a new raster");
5170 elog(NOTICE,
"Raster is empty. Returning an empty raster");
5172 PG_FREE_IF_COPY(pgraster, 0);
5176 if (NULL == pgrtn) {
5177 elog(ERROR,
"RASTER_mapAlgebraFct: Could not serialize raster");
5181 SET_VARSIZE(pgrtn, pgrtn->
size);
5182 PG_RETURN_POINTER(pgrtn);
5192 elog(NOTICE,
"Raster does not have the required band. Returning a raster "
5195 PG_FREE_IF_COPY(pgraster, 0);
5199 if (NULL == pgrtn) {
5200 elog(ERROR,
"RASTER_mapAlgebraFct: Could not serialize raster");
5204 SET_VARSIZE(pgrtn, pgrtn->
size);
5205 PG_RETURN_POINTER(pgrtn);
5210 if ( NULL ==
band ) {
5211 elog(NOTICE,
"Could not get the required band. Returning a raster "
5214 PG_FREE_IF_COPY(pgraster, 0);
5218 if (NULL == pgrtn) {
5219 elog(ERROR,
"RASTER_mapAlgebraFct: Could not serialize raster");
5223 SET_VARSIZE(pgrtn, pgrtn->
size);
5224 PG_RETURN_POINTER(pgrtn);
5230 POSTGIS_RT_DEBUG(3,
"RASTER_mapAlgebraFct: Getting NODATA value for band...");
5247 newinitialvalue = newnodatavalue;
5254 if (PG_ARGISNULL(2)) {
5259 strFromText = text_to_cstring(PG_GETARG_TEXT_P(2));
5262 if (newpixeltype ==
PT_END)
5266 if (newpixeltype ==
PT_END) {
5269 PG_FREE_IF_COPY(pgraster, 0);
5272 elog(ERROR,
"RASTER_mapAlgebraFct: Invalid pixeltype");
5280 if (PG_ARGISNULL(3)) {
5283 PG_FREE_IF_COPY(pgraster, 0);
5286 elog(ERROR,
"RASTER_mapAlgebraFct: Required function is missing. Returning NULL");
5290 oid = PG_GETARG_OID(3);
5291 if (oid == InvalidOid) {
5294 PG_FREE_IF_COPY(pgraster, 0);
5297 elog(ERROR,
"RASTER_mapAlgebraFct: Got invalid function object id. Returning NULL");
5301 fmgr_info(oid, &cbinfo);
5304 if (cbinfo.fn_retset) {
5307 PG_FREE_IF_COPY(pgraster, 0);
5310 elog(ERROR,
"RASTER_mapAlgebraFct: Function provided must return double precision not resultset");
5314 else if (cbinfo.fn_nargs < 2 || cbinfo.fn_nargs > 3) {
5317 PG_FREE_IF_COPY(pgraster, 0);
5320 elog(ERROR,
"RASTER_mapAlgebraFct: Function does not have two or three input parameters");
5324 if (cbinfo.fn_nargs == 2)
5329 if (func_volatile(oid) ==
'v') {
5330 elog(NOTICE,
"Function provided is VOLATILE. Unless required and for best performance, function should be IMMUTABLE or STABLE");
5334 InitFunctionCallInfoData(*cbdata, &cbinfo, 2, InvalidOid, NULL, NULL);
5336 cbdata->args[0].isnull =
FALSE;
5337 cbdata->args[1].isnull =
FALSE;
5338 cbdata->args[2].isnull =
FALSE;
5341 if (PG_ARGISNULL(4)) {
5342 if (cbinfo.fn_strict) {
5345 PG_FREE_IF_COPY(pgraster, 0);
5348 elog(ERROR,
"RASTER_mapAlgebraFct: Strict callback functions cannot have null parameters");
5352 cbdata->args[k].value = (Datum)NULL;
5353 cbdata->args[k].isnull =
TRUE;
5356 cbdata->args[k].value = PG_GETARG_DATUM(4);
5366 POSTGIS_RT_DEBUG(3,
"RASTER_mapAlgebraFct: Band is a nodata band, returning "
5367 "a raster filled with nodata");
5370 newinitialvalue,
TRUE, newnodatavalue, 0);
5373 PG_FREE_IF_COPY(pgraster, 0);
5378 if (NULL == pgrtn) {
5379 elog(ERROR,
"RASTER_mapAlgebraFct: Could not serialize raster");
5383 SET_VARSIZE(pgrtn, pgrtn->
size);
5384 PG_RETURN_POINTER(pgrtn);
5393 newinitialvalue,
TRUE, newnodatavalue, 0);
5397 if ( NULL == newband ) {
5398 elog(NOTICE,
"Could not modify band for new raster. Returning new "
5399 "raster with the original band");
5402 PG_FREE_IF_COPY(pgraster, 0);
5407 if (NULL == pgrtn) {
5408 elog(ERROR,
"RASTER_mapAlgebraFct: Could not serialize raster");
5412 SET_VARSIZE(pgrtn, pgrtn->
size);
5413 PG_RETURN_POINTER(pgrtn);
5419 for (
x = 0;
x < width;
x++) {
5420 for(
y = 0;
y < height;
y++) {
5428 if (
FLT_EQ(
r, newnodatavalue)) {
5429 if (cbinfo.fn_strict) {
5430 POSTGIS_RT_DEBUG(3,
"RASTER_mapAlgebraFct: Strict callbacks cannot accept NULL arguments, skipping NODATA cell.");
5433 cbdata->args[0].isnull =
TRUE;
5434 cbdata->args[0].value = (Datum)NULL;
5437 cbdata->args[0].isnull =
FALSE;
5438 cbdata->args[0].value = Float8GetDatum(
r);
5442 if (cbinfo.fn_nargs == 3) {
5446 d[0] = Int32GetDatum(
x+1);
5447 d[1] = Int32GetDatum(
y+1);
5449 a = construct_array(d, 2, INT4OID,
sizeof(
int32),
true,
'i');
5451 cbdata->args[1].isnull =
FALSE;
5452 cbdata->args[1].value = PointerGetDatum(a);
5458 tmpnewval = FunctionCallInvoke(cbdata);
5462 newval = newnodatavalue;
5465 newval = DatumGetFloat8(tmpnewval);
5479 POSTGIS_RT_DEBUG(3,
"RASTER_mapAlgebraFct: raster modified, serializing it.");
5483 PG_FREE_IF_COPY(pgraster, 0);
5494 SET_VARSIZE(pgrtn, pgrtn->
size);
5495 PG_RETURN_POINTER(pgrtn);
5510 int x,
y,
nband, width, height, ngbwidth, ngbheight, winwidth, winheight, u, v, nIndex, nNullItems;
5512 double newnodatavalue = 0.0;
5513 double newinitialvalue = 0.0;
5514 double newval = 0.0;
5519 LOCAL_FCINFO(cbdata, FUNC_MAX_ARGS);
5521 ArrayType * neighborDatum;
5522 char * strFromText = NULL;
5523 text * txtNodataMode = NULL;
5524 text * txtCallbackParam = NULL;
5526 float fltReplace = 0;
5527 bool valuereplace =
false, pixelreplace, nNodataOnly =
true, nNullSkip =
false;
5528 Datum * neighborData = NULL;
5529 bool * neighborNulls = NULL;
5530 int neighborDims[2];
5539 if (PG_ARGISNULL(0)) {
5540 elog(WARNING,
"Raster is NULL. Returning NULL");
5546 pgraster = (
rt_pgraster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
5550 PG_FREE_IF_COPY(pgraster, 0);
5551 elog(ERROR,
"RASTER_mapAlgebraFctNgb: Could not deserialize raster");
5559 if (PG_ARGISNULL(1))
5562 nband = PG_GETARG_INT32(1);
5567 POSTGIS_RT_DEBUG(3,
"RASTER_mapAlgebraFctNgb: Creating new empty raster...");
5578 if ( NULL == newrast ) {
5580 PG_FREE_IF_COPY(pgraster, 0);
5581 elog(ERROR,
"RASTER_mapAlgebraFctNgb: Could not create a new raster");
5606 elog(NOTICE,
"Raster is empty. Returning an empty raster");
5608 PG_FREE_IF_COPY(pgraster, 0);
5612 if (NULL == pgrtn) {
5613 elog(ERROR,
"RASTER_mapAlgebraFctNgb: Could not serialize raster");
5617 SET_VARSIZE(pgrtn, pgrtn->
size);
5618 PG_RETURN_POINTER(pgrtn);
5628 elog(NOTICE,
"Raster does not have the required band. Returning a raster "
5631 PG_FREE_IF_COPY(pgraster, 0);
5635 if (NULL == pgrtn) {
5636 elog(ERROR,
"RASTER_mapAlgebraFctNgb: Could not serialize raster");
5640 SET_VARSIZE(pgrtn, pgrtn->
size);
5641 PG_RETURN_POINTER(pgrtn);
5646 if ( NULL ==
band ) {
5647 elog(NOTICE,
"Could not get the required band. Returning a raster "
5650 PG_FREE_IF_COPY(pgraster, 0);
5654 if (NULL == pgrtn) {
5655 elog(ERROR,
"RASTER_mapAlgebraFctNgb: Could not serialize raster");
5659 SET_VARSIZE(pgrtn, pgrtn->
size);
5660 PG_RETURN_POINTER(pgrtn);
5666 POSTGIS_RT_DEBUG(3,
"RASTER_mapAlgebraFctNgb: Getting NODATA value for band...");
5683 newinitialvalue = newnodatavalue;
5690 if (PG_ARGISNULL(2)) {
5695 strFromText = text_to_cstring(PG_GETARG_TEXT_P(2));
5696 POSTGIS_RT_DEBUGF(3,
"RASTER_mapAlgebraFctNgb: Pixeltype parameter: %s", strFromText);
5699 if (newpixeltype ==
PT_END)
5703 if (newpixeltype ==
PT_END) {
5706 PG_FREE_IF_COPY(pgraster, 0);
5709 elog(ERROR,
"RASTER_mapAlgebraFctNgb: Invalid pixeltype");
5717 if (PG_ARGISNULL(5)) {
5720 PG_FREE_IF_COPY(pgraster, 0);
5723 elog(ERROR,
"RASTER_mapAlgebraFctNgb: Required function is missing");
5727 oid = PG_GETARG_OID(5);
5728 if (oid == InvalidOid) {
5731 PG_FREE_IF_COPY(pgraster, 0);
5734 elog(ERROR,
"RASTER_mapAlgebraFctNgb: Got invalid function object id");
5738 fmgr_info(oid, &cbinfo);
5741 if (cbinfo.fn_retset) {
5744 PG_FREE_IF_COPY(pgraster, 0);
5747 elog(ERROR,
"RASTER_mapAlgebraFctNgb: Function provided must return double precision not resultset");
5751 else if (cbinfo.fn_nargs != 3) {
5754 PG_FREE_IF_COPY(pgraster, 0);
5757 elog(ERROR,
"RASTER_mapAlgebraFctNgb: Function does not have three input parameters");
5761 if (func_volatile(oid) ==
'v') {
5762 elog(NOTICE,
"Function provided is VOLATILE. Unless required and for best performance, function should be IMMUTABLE or STABLE");
5766 InitFunctionCallInfoData(*cbdata, &cbinfo, 3, InvalidOid, NULL, NULL);
5767 cbdata->args[0].isnull =
FALSE;
5768 cbdata->args[1].isnull =
FALSE;
5769 cbdata->args[2].isnull =
FALSE;
5772 if (PG_ARGISNULL(7)) {
5773 if (cbinfo.fn_strict) {
5776 PG_FREE_IF_COPY(pgraster, 0);
5779 elog(ERROR,
"RASTER_mapAlgebraFctNgb: Strict callback functions cannot have NULL parameters");
5783 cbdata->args[2].value = (Datum)NULL;
5784 cbdata->args[2].isnull =
TRUE;
5787 cbdata->args[2].value = PG_GETARG_DATUM(7);
5797 POSTGIS_RT_DEBUG(3,
"RASTER_mapAlgebraFctNgb: Band is a nodata band, returning "
5798 "a raster filled with nodata");
5801 newinitialvalue,
TRUE, newnodatavalue, 0);
5804 PG_FREE_IF_COPY(pgraster, 0);
5809 if (NULL == pgrtn) {
5810 elog(ERROR,
"RASTER_mapAlgebraFctNgb: Could not serialize raster");
5814 SET_VARSIZE(pgrtn, pgrtn->
size);
5815 PG_RETURN_POINTER(pgrtn);
5824 newinitialvalue,
TRUE, newnodatavalue, 0);
5828 if ( NULL == newband ) {
5829 elog(NOTICE,
"Could not modify band for new raster. Returning new "
5830 "raster with the original band");
5833 PG_FREE_IF_COPY(pgraster, 0);
5838 if (NULL == pgrtn) {
5839 elog(ERROR,
"RASTER_mapAlgebraFctNgb: Could not serialize raster");
5843 SET_VARSIZE(pgrtn, pgrtn->
size);
5844 PG_RETURN_POINTER(pgrtn);
5848 if (PG_ARGISNULL(3) || PG_GETARG_INT32(3) <= 0) {
5849 elog(NOTICE,
"Neighborhood width is NULL or <= 0. Returning new "
5850 "raster with the original band");
5853 PG_FREE_IF_COPY(pgraster, 0);
5858 if (NULL == pgrtn) {
5859 elog(ERROR,
"RASTER_mapAlgebraFctNgb: Could not serialize raster");
5863 SET_VARSIZE(pgrtn, pgrtn->
size);
5864 PG_RETURN_POINTER(pgrtn);
5867 ngbwidth = PG_GETARG_INT32(3);
5868 winwidth = ngbwidth * 2 + 1;
5871 if (PG_ARGISNULL(4) || PG_GETARG_INT32(4) <= 0) {
5872 elog(NOTICE,
"Neighborhood height is NULL or <= 0. Returning new "
5873 "raster with the original band");
5876 PG_FREE_IF_COPY(pgraster, 0);
5881 if (NULL == pgrtn) {
5882 elog(ERROR,
"RASTER_mapAlgebraFctNgb: Could not serialize raster");
5886 SET_VARSIZE(pgrtn, pgrtn->
size);
5887 PG_RETURN_POINTER(pgrtn);
5890 ngbheight = PG_GETARG_INT32(4);
5891 winheight = ngbheight * 2 + 1;
5894 if (PG_ARGISNULL(6)) {
5895 elog(NOTICE,
"Neighborhood NODATA behavior defaulting to 'ignore'");
5896 txtNodataMode = cstring_to_text(
"ignore");
5899 txtNodataMode = PG_GETARG_TEXT_P(6);
5902 txtCallbackParam = (text*)palloc(VARSIZE(txtNodataMode));
5903 SET_VARSIZE(txtCallbackParam, VARSIZE(txtNodataMode));
5904 memcpy((
void *)VARDATA(txtCallbackParam), (
void *)VARDATA(txtNodataMode), VARSIZE(txtNodataMode) - VARHDRSZ);
5907 cbdata->args[1].value = PointerGetDatum(txtCallbackParam);
5909 strFromText = text_to_cstring(txtNodataMode);
5912 if (strcmp(strFromText,
"VALUE") == 0)
5913 valuereplace =
true;
5914 else if (strcmp(strFromText,
"IGNORE") != 0 && strcmp(strFromText,
"NULL") != 0) {
5916 if (sscanf(strFromText,
"%d", &intReplace) <= 0 && sscanf(strFromText,
"%f", &fltReplace) <= 0) {
5918 elog(NOTICE,
"Neighborhood NODATA mode is not recognized. Must be one of 'value', 'ignore', "
5919 "'NULL', or a numeric value. Returning new raster with the original band");
5922 pfree(txtCallbackParam);
5926 PG_FREE_IF_COPY(pgraster, 0);
5931 if (NULL == pgrtn) {
5932 elog(ERROR,
"RASTER_mapAlgebraFctNgb: Could not serialize raster");
5936 SET_VARSIZE(pgrtn, pgrtn->
size);
5937 PG_RETURN_POINTER(pgrtn);
5940 else if (strcmp(strFromText,
"NULL") == 0) {
5949 neighborData = (Datum *)palloc(winwidth * winheight *
sizeof(Datum));
5950 neighborNulls = (
bool *)palloc(winwidth * winheight *
sizeof(
bool));
5953 neighborDims[0] = winwidth;
5954 neighborDims[1] = winheight;
5961 get_typlenbyvalalign(FLOAT8OID, &typlen, &typbyval, &typalign);
5963 for (
x = 0 + ngbwidth;
x < width - ngbwidth;
x++) {
5964 for(
y = 0 + ngbheight;
y < height - ngbheight;
y++) {
5969 pixelreplace =
false;
5973 pixelreplace =
true;
5976 for (u =
x - ngbwidth; u <=
x + ngbwidth; u++) {
5977 for (v =
y - ngbheight; v <=
y + ngbheight; v++) {
5982 neighborData[nIndex] = Float8GetDatum((
double)
r);
5983 neighborNulls[nIndex] =
false;
5984 nNodataOnly =
false;
5988 if (valuereplace && pixelreplace) {
5990 neighborData[nIndex] = Float8GetDatum((
double)rpix);
5991 neighborNulls[nIndex] =
false;
5996 neighborData[nIndex] = PointerGetDatum(NULL);
5997 neighborNulls[nIndex] =
true;
6004 neighborData[nIndex] = PointerGetDatum(NULL);
6005 neighborNulls[nIndex] =
true;
6017 if (!(nNodataOnly ||
6018 (nNullSkip && nNullItems > 0) ||
6019 (valuereplace && nNullItems > 0))) {
6021 x,
y, winwidth, winheight);
6023 neighborDatum = construct_md_array((
void *)neighborData, neighborNulls, 2, neighborDims, neighborLbs,
6024 FLOAT8OID, typlen, typbyval, typalign);
6027 cbdata->args[0].value = PointerGetDatum(neighborDatum);
6030 tmpnewval = FunctionCallInvoke(cbdata);
6035 newval = newnodatavalue;
6038 newval = DatumGetFloat8(tmpnewval);
6054 pfree(neighborNulls);
6055 pfree(neighborData);
6057 pfree(txtCallbackParam);
6060 PG_FREE_IF_COPY(pgraster, 0);
6064 POSTGIS_RT_DEBUG(3,
"RASTER_mapAlgebraFctNgb: raster modified, serializing it.");
6075 SET_VARSIZE(pgrtn, pgrtn->
size);
6076 PG_RETURN_POINTER(pgrtn);
6079 #define ARGKWCOUNT 8
6087 const uint32_t set_count = 2;
6089 int pgrastpos[2] = {-1, -1};
6092 int _isempty[2] = {0};
6093 uint32_t bandindex[2] = {0};
6096 int _hasnodata[2] = {0};
6097 double _nodataval[2] = {0};
6098 double _offset[4] = {0.};
6099 double _rastoffset[2][4] = {{0.}};
6100 int _haspixel[2] = {0};
6101 double _pixel[2] = {0};
6102 int _pos[2][2] = {{0}};
6103 uint16_t _dim[2][2] = {{0}};
6105 char *pixtypename = NULL;
6107 char *extenttypename = NULL;
6112 uint16_t dim[2] = {0};
6115 double nodataval = 0;
6116 double gt[6] = {0.};
6118 Oid calltype = InvalidOid;
6120 const uint32_t spi_count = 3;
6121 uint16_t spi_exprpos[3] = {4, 7, 8};
6122 uint32_t spi_argcount[3] = {0};
6125 SPIPlanPtr spi_plan[3] = {NULL};
6126 uint16_t spi_empty = 0;
6127 Oid *argtype = NULL;
6128 uint8_t argpos[3][8] = {{0}};
6129 char *argkw[] = {
"[rast1.x]",
"[rast1.y]",
"[rast1.val]",
"[rast1]",
"[rast2.x]",
"[rast2.y]",
"[rast2.val]",
"[rast2]"};
6133 SPITupleTable *tuptable = NULL;
6136 bool isnull =
FALSE;
6137 int hasargval[3] = {0};
6138 double argval[3] = {0.};
6139 int hasnodatanodataval = 0;
6140 double nodatanodataval = 0;
6143 Oid ufc_noid = InvalidOid;
6145 LOCAL_FCINFO(ufc_info, FUNC_MAX_ARGS);
6147 int ufc_nullcount = 0;
6163 for (i = 0, j = 0; i < set_count; i++) {
6164 if (!PG_ARGISNULL(j)) {
6165 pgrast[i] = (
rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(j));
6172 for (k = 0; k <= i; k++) {
6173 if (k < i &&
rast[k] != NULL)
6175 if (pgrastpos[k] != -1)
6176 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6178 elog(ERROR,
"RASTER_mapAlgebra2: Could not deserialize the %s raster", i < 1 ?
"first" :
"second");
6186 if (!PG_ARGISNULL(j)) {
6187 bandindex[i] = PG_GETARG_INT32(j);
6200 if (
rast[0] == NULL &&
rast[1] == NULL) {
6201 elog(NOTICE,
"The two rasters provided are NULL. Returning NULL");
6202 for (k = 0; k < set_count; k++) {
6203 if (pgrastpos[k] != -1)
6204 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6210 if (_isempty[0] && _isempty[1]) {
6211 elog(NOTICE,
"The two rasters provided are empty. Returning empty raster");
6215 for (k = 0; k < set_count; k++) {
6216 if (
rast[k] != NULL)
6218 if (pgrastpos[k] != -1)
6219 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6221 elog(ERROR,
"RASTER_mapAlgebra2: Could not create empty raster");
6231 SET_VARSIZE(pgrtn, pgrtn->
size);
6232 PG_RETURN_POINTER(pgrtn);
6237 (
rast[0] == NULL || _isempty[0]) ||
6238 (
rast[1] == NULL || _isempty[1])
6241 if (
rast[0] == NULL || _isempty[0]) {
6254 if (_rast[i] != NULL)
6266 if (_rast[i] == NULL) {
6268 for (k = 0; k < set_count; k++) {
6269 if (pgrastpos[k] != -1)
6270 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6272 elog(ERROR,
"RASTER_mapAlgebra2: Could not create NODATA raster");
6306 for (k = 0; k < set_count; k++) {
6307 if (_rast[k] != NULL)
6309 if (pgrastpos[k] != -1)
6310 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6312 elog(ERROR,
"RASTER_mapAlgebra2: Could not test for alignment on the two rasters");
6316 elog(NOTICE,
"The two rasters provided do not have the same alignment. Returning NULL");
6317 for (k = 0; k < set_count; k++) {
6318 if (_rast[k] != NULL)
6320 if (pgrastpos[k] != -1)
6321 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6327 if (!PG_ARGISNULL(5)) {
6328 pixtypename = text_to_cstring(PG_GETARG_TEXT_P(5));
6331 if (pixtype ==
PT_END ) {
6332 for (k = 0; k < set_count; k++) {
6333 if (_rast[k] != NULL)
6335 if (pgrastpos[k] != -1)
6336 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6338 elog(ERROR,
"RASTER_mapAlgebra2: Invalid pixel type: %s", pixtypename);
6344 if (!PG_ARGISNULL(6)) {
6357 for (k = 0; k < set_count; k++) {
6358 if (_rast[k] != NULL)
6360 if (pgrastpos[k] != -1)
6361 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6363 elog(ERROR,
"RASTER_mapAlgebra2: Could not get output raster of correct extent");
6368 _rastoffset[0][0] = _offset[0];
6369 _rastoffset[0][1] = _offset[1];
6370 _rastoffset[1][0] = _offset[2];
6371 _rastoffset[1][1] = _offset[3];
6379 switch (extenttype) {
6389 (extenttype ==
ET_FIRST && i == 0) ||
6393 elog(NOTICE,
"The %s raster is NULL. Returning NULL", (i != 1 ?
"FIRST" :
"SECOND"));
6394 for (k = 0; k < set_count; k++) {
6395 if (_rast[k] != NULL)
6397 if (pgrastpos[k] != -1)
6398 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6406 elog(NOTICE,
"The %s raster does not have the band at index %d. Returning no band raster of correct extent",
6407 (i != 1 ?
"FIRST" :
"SECOND"), bandindex[i]
6410 for (k = 0; k < set_count; k++) {
6411 if (_rast[k] != NULL)
6413 if (pgrastpos[k] != -1)
6414 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6419 if (!pgrtn) PG_RETURN_NULL();
6421 SET_VARSIZE(pgrtn, pgrtn->
size);
6422 PG_RETURN_POINTER(pgrtn);
6430 _isempty[0] || _isempty[1] ||
6433 elog(NOTICE,
"The two rasters provided have no intersection. Returning no band raster");
6436 if (dim[0] || dim[1]) {
6441 for (k = 0; k < set_count; k++) {
6442 if (_rast[k] != NULL)
6444 if (pgrastpos[k] != -1)
6445 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6447 elog(ERROR,
"RASTER_mapAlgebra2: Could not create no band raster");
6455 for (k = 0; k < set_count; k++) {
6456 if (_rast[k] != NULL)
6458 if (pgrastpos[k] != -1)
6459 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6464 if (!pgrtn) PG_RETURN_NULL();
6466 SET_VARSIZE(pgrtn, pgrtn->
size);
6467 PG_RETURN_POINTER(pgrtn);
6472 for (k = 0; k < set_count; k++) {
6473 if (_rast[k] != NULL)
6475 if (pgrastpos[k] != -1)
6476 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6478 elog(ERROR,
"RASTER_mapAlgebra2: ET_LAST and ET_CUSTOM are not implemented");
6488 elog(NOTICE,
"The two rasters provided do not have the respectively specified band indices. Returning no band raster of correct extent");
6490 for (k = 0; k < set_count; k++) {
6491 if (_rast[k] != NULL)
6493 if (pgrastpos[k] != -1)
6494 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6499 if (!pgrtn) PG_RETURN_NULL();
6501 SET_VARSIZE(pgrtn, pgrtn->
size);
6502 PG_RETURN_POINTER(pgrtn);
6506 for (i = 0; i < set_count; i++) {
6515 if (_band[i] == NULL) {
6516 for (k = 0; k < set_count; k++) {
6517 if (_rast[k] != NULL)
6519 if (pgrastpos[k] != -1)
6520 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6523 elog(ERROR,
"RASTER_mapAlgebra2: Could not get band %d of the %s raster",
6525 (i < 1 ?
"FIRST" :
"SECOND")
6537 if ((extenttype ==
ET_SECOND && !_isempty[1]) || _isempty[0])
6544 if (extenttype ==
ET_SECOND && !_isempty[1] && _hasnodata[1]) {
6545 nodataval = _nodataval[1];
6547 else if (!_isempty[0] && _hasnodata[0]) {
6548 nodataval = _nodataval[0];
6550 else if (!_isempty[1] && _hasnodata[1]) {
6551 nodataval = _nodataval[1];
6554 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));
6566 for (k = 0; k < set_count; k++) {
6567 if (_rast[k] != NULL)
6569 if (pgrastpos[k] != -1)
6570 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6573 elog(ERROR,
"RASTER_mapAlgebra2: Could not add new band to output raster");
6580 for (k = 0; k < set_count; k++) {
6581 if (_rast[k] != NULL)
6583 if (pgrastpos[k] != -1)
6584 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6587 elog(ERROR,
"RASTER_mapAlgebra2: Could not get newly added band of output raster");
6592 (
int) _rastoffset[0][0],
6593 (
int) _rastoffset[0][1],
6594 (
int) _rastoffset[1][0],
6595 (
int) _rastoffset[1][1]
6598 POSTGIS_RT_DEBUGF(4,
"metadata = (%f, %f, %d, %d, %f, %f, %f, %f, %d)",
6615 calltype = get_fn_expr_argtype(fcinfo->flinfo, 4);
6622 if (SPI_connect() != SPI_OK_CONNECT) {
6623 for (k = 0; k < set_count; k++) {
6624 if (_rast[k] != NULL)
6626 if (pgrastpos[k] != -1)
6627 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6630 elog(ERROR,
"RASTER_mapAlgebra2: Could not connect to the SPI manager");
6635 memset(hasargval, 0,
sizeof(
int) * spi_count);
6645 for (i = 0; i < spi_count; i++) {
6646 if (!PG_ARGISNULL(spi_exprpos[i])) {
6648 char place[5] =
"$1";
6649 expr = text_to_cstring(PG_GETARG_TEXT_P(spi_exprpos[i]));
6663 sprintf(place,
"$%d", k);
6669 len = strlen(
"SELECT (") + strlen(expr) + strlen(
")::double precision");
6670 sql = (
char *) palloc(len + 1);
6673 for (k = 0; k < spi_count; k++) SPI_freeplan(spi_plan[k]);
6676 for (k = 0; k < set_count; k++) {
6677 if (_rast[k] != NULL)
6679 if (pgrastpos[k] != -1)
6680 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6684 elog(ERROR,
"RASTER_mapAlgebra2: Could not allocate memory for expression parameter %d", spi_exprpos[i]);
6688 memcpy(
sql,
"SELECT (", strlen(
"SELECT ("));
6689 memcpy(
sql + strlen(
"SELECT ("), expr, strlen(expr));
6690 memcpy(
sql + strlen(
"SELECT (") + strlen(expr),
")::double precision", strlen(
")::double precision"));
6696 if (spi_argcount[i]) {
6697 argtype = (Oid *) palloc(spi_argcount[i] *
sizeof(Oid));
6698 if (argtype == NULL) {
6701 for (k = 0; k < spi_count; k++) SPI_freeplan(spi_plan[k]);
6704 for (k = 0; k < set_count; k++) {
6705 if (_rast[k] != NULL)
6707 if (pgrastpos[k] != -1)
6708 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6712 elog(ERROR,
"RASTER_mapAlgebra2: Could not allocate memory for prepared plan argtypes of expression parameter %d", spi_exprpos[i]);
6718 if (argpos[i][j] < 1)
continue;
6722 (strstr(argkw[j],
"[rast1.x]") != NULL) ||
6723 (strstr(argkw[j],
"[rast1.y]") != NULL) ||
6724 (strstr(argkw[j],
"[rast2.x]") != NULL) ||
6725 (strstr(argkw[j],
"[rast2.y]") != NULL)
6727 argtype[k] = INT4OID;
6731 argtype[k] = FLOAT8OID;
6737 spi_plan[i] = SPI_prepare(
sql, spi_argcount[i], argtype);
6740 if (spi_plan[i] == NULL) {
6743 for (k = 0; k < spi_count; k++) SPI_freeplan(spi_plan[k]);
6746 for (k = 0; k < set_count; k++) {
6747 if (_rast[k] != NULL)
6749 if (pgrastpos[k] != -1)
6750 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6754 elog(ERROR,
"RASTER_mapAlgebra2: Could not create prepared plan of expression parameter %d", spi_exprpos[i]);
6760 err = SPI_execute(
sql,
TRUE, 0);
6761 if (err != SPI_OK_SELECT || SPI_tuptable == NULL || SPI_processed != 1) {
6764 for (k = 0; k < spi_count; k++) SPI_freeplan(spi_plan[k]);
6767 for (k = 0; k < set_count; k++) {
6768 if (_rast[k] != NULL)
6770 if (pgrastpos[k] != -1)
6771 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6775 elog(ERROR,
"RASTER_mapAlgebra2: Could not evaluate expression parameter %d", spi_exprpos[i]);
6780 tupdesc = SPI_tuptable->tupdesc;
6781 tuptable = SPI_tuptable;
6782 tuple = tuptable->vals[0];
6784 datum = SPI_getbinval(tuple, tupdesc, 1, &isnull);
6785 if (SPI_result == SPI_ERROR_NOATTRIBUTE) {
6788 if (SPI_tuptable) SPI_freetuptable(tuptable);
6789 for (k = 0; k < spi_count; k++) SPI_freeplan(spi_plan[k]);
6792 for (k = 0; k < set_count; k++) {
6793 if (_rast[k] != NULL)
6795 if (pgrastpos[k] != -1)
6796 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6800 elog(ERROR,
"RASTER_mapAlgebra2: Could not get result of expression parameter %d", spi_exprpos[i]);
6806 argval[i] = DatumGetFloat8(datum);
6809 if (SPI_tuptable) SPI_freetuptable(tuptable);
6819 if (!PG_ARGISNULL(9)) {
6820 hasnodatanodataval = 1;
6821 nodatanodataval = PG_GETARG_FLOAT8(9);
6824 hasnodatanodataval = 0;
6827 case REGPROCEDUREOID: {
6829 if (!PG_ARGISNULL(4)) {
6832 ufc_noid = PG_GETARG_OID(4);
6835 fmgr_info(ufc_noid, &ufl_info);
6839 if (ufl_info.fn_retset) {
6843 else if (ufl_info.fn_nargs < 3 || ufl_info.fn_nargs > 4) {
6853 for (k = 0; k < set_count; k++) {
6854 if (_rast[k] != NULL)
6856 if (pgrastpos[k] != -1)
6857 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6862 elog(ERROR,
"RASTER_mapAlgebra2: Function provided must have three or four input parameters");
6864 elog(ERROR,
"RASTER_mapAlgebra2: Function provided must return double precision not resultset");
6868 if (func_volatile(ufc_noid) ==
'v') {
6869 elog(NOTICE,
"Function provided is VOLATILE. Unless required and for best performance, function should be IMMUTABLE or STABLE");
6873 InitFunctionCallInfoData(
6874 *ufc_info, &ufl_info, ufl_info.fn_nargs, InvalidOid, NULL, NULL);
6875 ufc_info->args[0].isnull =
FALSE;
6876 ufc_info->args[1].isnull =
FALSE;
6877 ufc_info->args[2].isnull =
FALSE;
6878 if (ufl_info.fn_nargs == 4)
6879 ufc_info->args[3].isnull =
FALSE;
6881 if (ufl_info.fn_nargs != 4)
6885 if (!PG_ARGISNULL(7))
6887 ufc_info->args[k].value = PG_GETARG_DATUM(7);
6891 ufc_info->args[k].value = (Datum)NULL;
6892 ufc_info->args[k].isnull =
TRUE;
6899 for (k = 0; k < set_count; k++) {
6900 if (_rast[k] != NULL)
6902 if (pgrastpos[k] != -1)
6903 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6906 elog(ERROR,
"RASTER_mapAlgebra2: Invalid data type for expression or userfunction");
6914 (calltype == TEXTOID) && (
6915 (spi_empty != spi_count) || hasnodatanodataval
6918 (calltype == REGPROCEDUREOID) && (ufc_noid != InvalidOid)
6920 for (
x = 0;
x < dim[0];
x++) {
6921 for (
y = 0;
y < dim[1];
y++) {
6924 for (i = 0; i < set_count; i++) {
6929 _x = (int)
x - (
int)_rastoffset[i][0];
6930 _y = (int)
y - (
int)_rastoffset[i][1];
6933 _pos[i][0] = _x + 1;
6934 _pos[i][1] = _y + 1;
6937 if (_band[i] == NULL) {
6938 if (!_hasnodata[i]) {
6940 _pixel[i] = _nodataval[i];
6945 (_x >= 0 && _x < _dim[i][0]) &&
6946 (_y >= 0 && _y < _dim[i][1])
6951 if (calltype == TEXTOID) {
6952 for (k = 0; k < spi_count; k++) SPI_freeplan(spi_plan[k]);
6956 for (k = 0; k < set_count; k++) {
6957 if (_rast[k] != NULL)
6959 if (pgrastpos[k] != -1)
6960 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6964 elog(ERROR,
"RASTER_mapAlgebra2: Could not get pixel of %s raster", (i < 1 ?
"FIRST" :
"SECOND"));
6968 if (!_hasnodata[i] || !isnodata)
6987 if (!_haspixel[0] && !_haspixel[1])
6991 else if (_haspixel[0] && !_haspixel[1])
6995 else if (!_haspixel[0] && _haspixel[1])
7004 if (hasnodatanodataval) {
7006 pixel = nodatanodataval;
7010 else if (hasargval[i]) {
7015 else if (spi_plan[i] != NULL) {
7019 memset(values, (Datum) NULL,
sizeof(Datum) *
ARGKWCOUNT);
7024 if (spi_argcount[i]) {
7028 if (idx < 1)
continue;
7031 if (strstr(argkw[j],
"[rast1.x]") != NULL) {
7032 values[idx] = _pos[0][0];
7034 else if (strstr(argkw[j],
"[rast1.y]") != NULL) {
7035 values[idx] = _pos[0][1];
7038 (strstr(argkw[j],
"[rast1.val]") != NULL) ||
7039 (strstr(argkw[j],
"[rast1]") != NULL)
7041 if (_isempty[0] || !_haspixel[0])
7044 values[idx] = Float8GetDatum(_pixel[0]);
7046 else if (strstr(argkw[j],
"[rast2.x]") != NULL) {
7047 values[idx] = _pos[1][0];
7049 else if (strstr(argkw[j],
"[rast2.y]") != NULL) {
7050 values[idx] = _pos[1][1];
7053 (strstr(argkw[j],
"[rast2.val]") != NULL) ||
7054 (strstr(argkw[j],
"[rast2]") != NULL)
7056 if (_isempty[1] || !_haspixel[1])
7059 values[idx] = Float8GetDatum(_pixel[1]);
7065 err = SPI_execute_plan(spi_plan[i], values, nulls,
TRUE, 1);
7066 if (err != SPI_OK_SELECT || SPI_tuptable == NULL || SPI_processed != 1) {
7068 for (k = 0; k < spi_count; k++) SPI_freeplan(spi_plan[k]);
7071 for (k = 0; k < set_count; k++) {
7072 if (_rast[k] != NULL)
7074 if (pgrastpos[k] != -1)
7075 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
7079 elog(ERROR,
"RASTER_mapAlgebra2: Unexpected error when running prepared statement %d", i);
7084 tupdesc = SPI_tuptable->tupdesc;
7085 tuptable = SPI_tuptable;
7086 tuple = tuptable->vals[0];
7088 datum = SPI_getbinval(tuple, tupdesc, 1, &isnull);
7089 if (SPI_result == SPI_ERROR_NOATTRIBUTE) {
7091 if (SPI_tuptable) SPI_freetuptable(tuptable);
7092 for (k = 0; k < spi_count; k++) SPI_freeplan(spi_plan[k]);
7095 for (k = 0; k < set_count; k++) {
7096 if (_rast[k] != NULL)
7098 if (pgrastpos[k] != -1)
7099 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
7103 elog(ERROR,
"RASTER_mapAlgebra2: Could not get result of prepared statement %d", i);
7109 pixel = DatumGetFloat8(datum);
7112 if (SPI_tuptable) SPI_freetuptable(tuptable);
7115 case REGPROCEDUREOID: {
7120 for (i = 0; i < set_count; i++) {
7121 ufc_info->args[i].value = Float8GetDatum(_pixel[i]);
7123 ufc_info->args[i].isnull =
FALSE;
7127 ufc_info->args[i].isnull =
TRUE;
7134 if (ufl_info.fn_strict && ufc_nullcount)
7138 if (ufl_info.fn_nargs == 4) {
7141 for (i = 0; i < set_count; i++) {
7143 d[0] = Int32GetDatum(_pos[i][0]);
7144 d[1] = Int32GetDatum(_pos[i][1]);
7147 d[2] = Int32GetDatum(_pos[i][0]);
7148 d[3] = Int32GetDatum(_pos[i][1]);
7152 a = construct_array(d, 4, INT4OID,
sizeof(
int32),
true,
'i');
7153 ufc_info->args[2].value = PointerGetDatum(a);
7154 ufc_info->args[2].isnull =
FALSE;
7157 datum = FunctionCallInvoke(ufc_info);
7160 if (!ufc_info->isnull)
7163 pixel = DatumGetFloat8(datum);
7172 if (calltype == TEXTOID) {
7173 for (k = 0; k < spi_count; k++) SPI_freeplan(spi_plan[k]);
7177 for (k = 0; k < set_count; k++) {
7178 if (_rast[k] != NULL)
7180 if (pgrastpos[k] != -1)
7181 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
7185 elog(ERROR,
"RASTER_mapAlgebra2: Could not set pixel value of output raster");
7197 if (calltype == TEXTOID) {
7198 for (i = 0; i < spi_count; i++) {
7199 if (spi_plan[i] != NULL) SPI_freeplan(spi_plan[i]);
7204 for (k = 0; k < set_count; k++) {
7205 if (_rast[k] != NULL)
7207 if (pgrastpos[k] != -1)
7208 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
7213 if (!pgrtn) PG_RETURN_NULL();
7217 SET_VARSIZE(pgrtn, pgrtn->
size);
7218 PG_RETURN_POINTER(pgrtn);
int32_t gserialized_get_srid(const GSERIALIZED *g)
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)
#define LWSIZE_GET(varsize)
Macro for reading the size from the GSERIALIZED size attribute.
#define SRID_UNKNOWN
Unknown SRID value.
lwvarlena_t * lwgeom_to_wkb_varlena(const LWGEOM *geom, uint8_t variant)
int32_t clamp_srid(int32_t srid)
Return a valid SRID from an arbitrary integer Raises a notice if what comes out is different from wha...
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_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.
static int lwgeom_is_empty(const LWGEOM *geom)
Return true or false depending on whether a geometry is an "empty" geometry (no vertices members)
raster
Be careful!! Zeros function's input parameter can be a (height x width) array, not (width x height): ...
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)
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::@12 method
struct rt_reclassexpr_t::rt_reclassrange src
struct rt_reclassexpr_t::rt_reclassrange dst
rtpg_nmapalgebra_callback_arg callback
union rtpg_nmapalgebra_callback_arg::@23 ufc_info_data
FunctionCallInfo ufc_info
FunctionCallInfoBaseData fcinfo
rtpg_nmapalgebraexpr_callback_arg callback
rtpg_nmapalgebra_arg bandarg
struct rtpg_nmapalgebraexpr_callback_arg::@26 kw
struct rtpg_nmapalgebraexpr_callback_arg::@24 expr[3]
struct rtpg_nmapalgebraexpr_callback_arg::@25 nodatanodata
rtpg_union_band_arg bandarg
rtpg_union_type uniontype