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++) {
469 _pos[i] = (Datum)arg->
src_pixel[z][0] + 1;
472 _pos[i] = (Datum)arg->
src_pixel[z][1] + 1;
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");
968 if (itrset == NULL) {
970 elog(ERROR,
"RASTER_nMapAlgebra: Could not allocate memory for iterator arguments");
999 elog(ERROR,
"RASTER_nMapAlgebra: Could not run raster iterator function");
1013 SET_VARSIZE(pgraster, pgraster->
size);
1014 PG_RETURN_POINTER(pgraster);
1058 elog(ERROR,
"rtpg_nmapalgebraexpr_arg_init: Could not allocate memory for arguments");
1064 elog(ERROR,
"rtpg_nmapalgebraexpr_arg_init: Could not allocate memory for arg->bandarg");
1077 elog(ERROR,
"rtpg_nmapalgebraexpr_arg_init: Could not allocate memory for spi_argpos");
1108 double *
value,
int *nodata
1111 SPIPlanPtr plan = NULL;
1131 else if (arg->
nodata[0][0][0] && !arg->
nodata[1][0][0]) {
1141 else if (!arg->
nodata[0][0][0] && arg->
nodata[1][0][0]) {
1168 if (arg->
nodata[0][0][0]) {
1204 SPITupleTable *tuptable = NULL;
1207 bool isnull =
FALSE;
1212 memset(values, (Datum) NULL,
sizeof(Datum) * callback->
kw.
count);
1213 memset(nulls,
FALSE,
sizeof(
char) * callback->
kw.
count);
1218 for (i = 0; i < callback->
kw.
count; i++) {
1220 if (idx < 1)
continue;
1226 values[idx] = Int32GetDatum(arg->
src_pixel[0][0] + 1);
1230 values[idx] = Int32GetDatum(arg->
src_pixel[0][1] + 1);
1236 if (!arg->
nodata[0][0][0])
1237 values[idx] = Float8GetDatum(arg->
values[0][0][0]);
1244 values[idx] = Int32GetDatum(arg->
src_pixel[0][0] + 1);
1248 values[idx] = Int32GetDatum(arg->
src_pixel[0][1] + 1);
1254 if (!arg->
nodata[0][0][0])
1255 values[idx] = Float8GetDatum(arg->
values[0][0][0]);
1262 values[idx] = Int32GetDatum(arg->
src_pixel[1][0] + 1);
1266 values[idx] = Int32GetDatum(arg->
src_pixel[1][1] + 1);
1272 if (!arg->
nodata[1][0][0])
1273 values[idx] = Float8GetDatum(arg->
values[1][0][0]);
1283 err = SPI_execute_plan(plan, values, nulls,
TRUE, 1);
1284 if (err != SPI_OK_SELECT || SPI_tuptable == NULL || SPI_processed != 1) {
1285 elog(ERROR,
"rtpg_nmapalgebraexpr_callback: Unexpected error when running prepared statement %d",
id);
1290 tupdesc = SPI_tuptable->tupdesc;
1291 tuptable = SPI_tuptable;
1292 tuple = tuptable->vals[0];
1294 datum = SPI_getbinval(tuple, tupdesc, 1, &isnull);
1295 if (SPI_result == SPI_ERROR_NOATTRIBUTE) {
1296 if (SPI_tuptable) SPI_freetuptable(tuptable);
1297 elog(ERROR,
"rtpg_nmapalgebraexpr_callback: Could not get result of prepared statement %d",
id);
1302 *
value = DatumGetFloat8(datum);
1322 if (SPI_tuptable) SPI_freetuptable(tuptable);
1332 MemoryContext mainMemCtx = CurrentMemoryContext;
1335 uint16_t exprpos[3] = {1, 4, 5};
1349 SPITupleTable *tuptable = NULL;
1352 bool isnull =
FALSE;
1358 const int argkwcount = 12;
1374 if (PG_ARGISNULL(0))
1380 elog(ERROR,
"RASTER_nMapAlgebraExpr: Could not initialize argument structure");
1387 elog(ERROR,
"RASTER_nMapAlgebra: Could not process rastbandarg");
1391 POSTGIS_RT_DEBUGF(4,
"allnull, allempty, noband = %d, %d, %d", allnull, allempty, noband);
1395 elog(NOTICE,
"All input rasters are NULL. Returning NULL");
1407 if (!PG_ARGISNULL(2)) {
1414 elog(ERROR,
"RASTER_nMapAlgebraExpr: Invalid pixel type: %s", pixtypename);
1421 if (!PG_ARGISNULL(3)) {
1427 if (numraster < 2) {
1428 elog(NOTICE,
"CUSTOM extent type not supported. Defaulting to FIRST");
1432 elog(NOTICE,
"CUSTOM extent type not supported. Defaulting to INTERSECTION");
1436 else if (numraster < 2)
1442 if (!PG_ARGISNULL(6)) {
1450 elog(NOTICE,
"All input rasters are empty. Returning empty raster");
1455 elog(NOTICE,
"All input rasters do not have bands at indicated indexes. Returning empty raster");
1463 elog(ERROR,
"RASTER_nMapAlgebraExpr: Could not create empty raster");
1469 if (!pgraster) PG_RETURN_NULL();
1471 SET_VARSIZE(pgraster, pgraster->
size);
1472 PG_RETURN_POINTER(pgraster);
1476 if (SPI_connect() != SPI_OK_CONNECT) {
1478 elog(ERROR,
"RASTER_nMapAlgebraExpr: Could not connect to the SPI manager");
1494 char place[12] =
"$1";
1496 if (PG_ARGISNULL(exprpos[i]))
1502 for (j = 0, k = 1; j < argkwcount; j++) {
1514 sprintf(place,
"$%d", k);
1520 len = strlen(
"SELECT (") + strlen(expr) + strlen(
")::double precision");
1521 sql = (
char *) palloc(len + 1);
1525 elog(ERROR,
"RASTER_nMapAlgebraExpr: Could not allocate memory for expression parameter %d", exprpos[i]);
1529 memcpy(
sql,
"SELECT (", strlen(
"SELECT ("));
1530 memcpy(
sql + strlen(
"SELECT ("), expr, strlen(expr));
1531 memcpy(
sql + strlen(
"SELECT (") + strlen(expr),
")::double precision", strlen(
")::double precision"));
1540 if (argtype == NULL) {
1544 elog(ERROR,
"RASTER_nMapAlgebraExpr: Could not allocate memory for prepared plan argtypes of expression parameter %d", exprpos[i]);
1549 for (j = 0, k = 0; j < argkwcount; j++) {
1554 (strstr(argkw[j],
"[rast.x]") != NULL) ||
1555 (strstr(argkw[j],
"[rast.y]") != NULL) ||
1556 (strstr(argkw[j],
"[rast1.x]") != NULL) ||
1557 (strstr(argkw[j],
"[rast1.y]") != NULL) ||
1558 (strstr(argkw[j],
"[rast2.x]") != NULL) ||
1559 (strstr(argkw[j],
"[rast2.y]") != NULL)
1561 argtype[k] = INT4OID;
1564 argtype[k] = FLOAT8OID;
1576 elog(ERROR,
"RASTER_nMapAlgebraExpr: Could not create prepared plan of expression parameter %d", exprpos[i]);
1582 POSTGIS_RT_DEBUGF(3,
"expression parameter %d has no args, simply executing", exprpos[i]);
1583 err = SPI_execute(
sql,
TRUE, 0);
1586 if (err != SPI_OK_SELECT || SPI_tuptable == NULL || SPI_processed != 1) {
1589 elog(ERROR,
"RASTER_nMapAlgebraExpr: Could not evaluate expression parameter %d", exprpos[i]);
1594 tupdesc = SPI_tuptable->tupdesc;
1595 tuptable = SPI_tuptable;
1596 tuple = tuptable->vals[0];
1598 datum = SPI_getbinval(tuple, tupdesc, 1, &isnull);
1599 if (SPI_result == SPI_ERROR_NOATTRIBUTE) {
1600 if (SPI_tuptable) SPI_freetuptable(tuptable);
1603 elog(ERROR,
"RASTER_nMapAlgebraExpr: Could not get result of expression parameter %d", exprpos[i]);
1612 if (SPI_tuptable) SPI_freetuptable(tuptable);
1632 for (i = 0; i < numraster; i++) {
1656 if (itrset == NULL) {
1659 elog(ERROR,
"RASTER_nMapAlgebra: Could not allocate memory for iterator arguments");
1664 for (i = 0; i < numraster; i++) {
1688 elog(ERROR,
"RASTER_nMapAlgebraExpr: Could not run raster iterator function");
1691 else if (
raster == NULL) {
1697 MemoryContextSwitchTo(mainMemCtx);
1708 SET_VARSIZE(pgraster, pgraster->
size);
1709 PG_RETURN_POINTER(pgraster);
1729 assert(cutype && strlen(cutype) > 0);
1731 if (strcmp(cutype,
"LAST") == 0)
1733 else if (strcmp(cutype,
"FIRST") == 0)
1735 else if (strcmp(cutype,
"MIN") == 0)
1737 else if (strcmp(cutype,
"MAX") == 0)
1739 else if (strcmp(cutype,
"COUNT") == 0)
1741 else if (strcmp(cutype,
"SUM") == 0)
1743 else if (strcmp(cutype,
"MEAN") == 0)
1745 else if (strcmp(cutype,
"RANGE") == 0)
1772 for (i = 0; i < arg->
numband; i++) {
1796 double *
value,
int *nodata
1808 elog(ERROR,
"rtpg_union_callback: Invalid arguments passed to callback");
1824 else if (!arg->
nodata[0][0][0] && arg->
nodata[1][0][0]) {
1830 else if (arg->
nodata[0][0][0] && !arg->
nodata[1][0][0]) {
1858 else if (!arg->
nodata[0][0][0] && arg->
nodata[1][0][0])
1861 else if (arg->
nodata[0][0][0] && !arg->
nodata[1][0][0])
1887 double *
value,
int *nodata
1897 elog(ERROR,
"rtpg_union_mean_callback: Invalid arguments passed to callback");
1920 double *
value,
int *nodata
1930 elog(ERROR,
"rtpg_union_range_callback: Invalid arguments passed to callback");
1963 HeapTupleHeader tup;
1969 char *utypename = NULL;
1972 etype = ARR_ELEMTYPE(array);
1973 get_typlenbyvalalign(etype, &typlen, &typbyval, &typalign);
1978 typlen, typbyval, typalign,
1983 elog(ERROR,
"rtpg_union_unionarg_process: Invalid argument for unionarg");
1991 elog(ERROR,
"rtpg_union_unionarg_process: Could not allocate memory for band information");
1996 for (i = 0; i < n; i++) {
2005 tup = (HeapTupleHeader) DatumGetPointer(e[i]);
2007 elog(ERROR,
"rtpg_union_unionarg_process: Invalid argument for unionarg");
2012 tupv = GetAttributeByName(tup,
"nband", &isnull);
2015 elog(NOTICE,
"First argument (nband) of unionarg is NULL. Assuming nband = %d",
nband);
2018 nband = DatumGetInt32(tupv);
2021 elog(ERROR,
"rtpg_union_unionarg_process: Band number must be greater than zero (1-based)");
2026 tupv = GetAttributeByName(tup,
"uniontype", &isnull);
2028 elog(NOTICE,
"Second argument (uniontype) of unionarg is NULL. Assuming uniontype = LAST");
2053 elog(ERROR,
"rtpg_union_unionarg_process: Could not reallocate memory for band information");
2070 if (numbands <= arg->numband)
2080 elog(ERROR,
"rtpg_union_noarg: Could not reallocate memory for band information");
2086 for (; i < arg->
numband; i++) {
2094 elog(ERROR,
"rtpg_union_noarg: Could not allocate memory for working rasters");
2103 elog(ERROR,
"rtpg_union_noarg: Could not create working raster");
2116 MemoryContext aggcontext;
2117 MemoryContext oldcontext;
2127 int isempty[2] = {0};
2128 int hasband[2] = {0};
2130 double _offset[4] = {0.};
2138 char *utypename = NULL;
2142 double nodataval = 0;
2148 uint16_t _dim[2] = {0};
2155 if (!AggCheckCallContext(fcinfo, &aggcontext)) {
2156 elog(ERROR,
"RASTER_union_transfn: Cannot be called in a non-aggregate context");
2161 oldcontext = MemoryContextSwitchTo(aggcontext);
2163 if (PG_ARGISNULL(0)) {
2168 MemoryContextSwitchTo(oldcontext);
2169 elog(ERROR,
"RASTER_union_transfn: Could not allocate memory for state variable");
2185 if (!PG_ARGISNULL(1)) {
2187 pgraster = (
rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
2194 PG_FREE_IF_COPY(pgraster, 1);
2196 MemoryContextSwitchTo(oldcontext);
2197 elog(ERROR,
"RASTER_union_transfn: Could not deserialize raster");
2210 if (!PG_ARGISNULL(2)) {
2211 Oid calltype = get_fn_expr_argtype(fcinfo->flinfo, 2);
2246 if (numband > idx) {
2264 PG_FREE_IF_COPY(pgraster, 1);
2267 MemoryContextSwitchTo(oldcontext);
2268 elog(ERROR,
"RASTER_union_transfn: Could not allocate memory for band information");
2273 for (i = idx; i < iwr->
numband; i++) {
2297 nband = PG_GETARG_INT32(2);
2303 PG_FREE_IF_COPY(pgraster, 1);
2306 MemoryContextSwitchTo(oldcontext);
2307 elog(ERROR,
"RASTER_union_transfn: Band number must be greater than zero (1-based)");
2318 PG_FREE_IF_COPY(pgraster, 1);
2321 MemoryContextSwitchTo(oldcontext);
2322 elog(ERROR,
"RASTER_union_transfn: Could not allocate memory for band information");
2343 PG_FREE_IF_COPY(pgraster, 1);
2346 MemoryContextSwitchTo(oldcontext);
2347 elog(ERROR,
"RASTER_union_transfn: Could not process unionarg");
2356 if (nargs > 3 && !PG_ARGISNULL(3)) {
2371 for (i = 0; i < iwr->
numband; i++) {
2386 PG_FREE_IF_COPY(pgraster, 1);
2389 MemoryContextSwitchTo(oldcontext);
2390 elog(ERROR,
"RASTER_union_transfn: Could not allocate memory for working raster(s)");
2405 PG_FREE_IF_COPY(pgraster, 1);
2408 MemoryContextSwitchTo(oldcontext);
2409 elog(ERROR,
"RASTER_union_transfn: Could not create working raster");
2426 PG_FREE_IF_COPY(pgraster, 1);
2429 MemoryContextSwitchTo(oldcontext);
2430 elog(ERROR,
"RASTER_union_transfn: Could not check and balance number of bands");
2437 if (itrset == NULL) {
2442 PG_FREE_IF_COPY(pgraster, 1);
2445 MemoryContextSwitchTo(oldcontext);
2446 elog(ERROR,
"RASTER_union_transfn: Could not allocate memory for iterator arguments");
2451 for (i = 0; i < iwr->
numband; i++) {
2471 if (!isempty[0] && hasband[0])
2473 else if (!isempty[1] && hasband[1])
2480 if (_band != NULL) {
2518 itrset[0].
nband = 0;
2534 if (!isempty[0] && !isempty[1] && hasband[0] && hasband[1]) {
2548 PG_FREE_IF_COPY(pgraster, 1);
2551 MemoryContextSwitchTo(oldcontext);
2552 elog(ERROR,
"RASTER_union_transfn: Could not create internal raster");
2556 _offset[0], _offset[1], _offset[2], _offset[3]);
2563 double igt[6] = {0};
2579 hasnodata, nodataval,
2588 PG_FREE_IF_COPY(pgraster, 1);
2591 MemoryContextSwitchTo(oldcontext);
2592 elog(ERROR,
"RASTER_union_transfn: Could not add new band to internal raster");
2600 for (
y = 0;
y < _dim[1];
y++) {
2601 POSTGIS_RT_DEBUGF(4,
"Getting pixel line of working raster at (x, y, length) = (0, %d, %d)",
y, _dim[0]);
2615 PG_FREE_IF_COPY(pgraster, 1);
2618 MemoryContextSwitchTo(oldcontext);
2619 elog(ERROR,
"RASTER_union_transfn: Could not get pixel line from band of working raster");
2623 POSTGIS_RT_DEBUGF(4,
"Setting pixel line at (x, y, length) = (%d, %d, %d)", (
int) _offset[0], (
int) _offset[1] +
y, nvals);
2626 (
int) _offset[0], (
int) _offset[1] +
y,
2636 PG_FREE_IF_COPY(pgraster, 1);
2639 MemoryContextSwitchTo(oldcontext);
2640 elog(ERROR,
"RASTER_union_transfn: Could not set pixel line to band of internal raster");
2656 hasnodata, nodataval,
2673 PG_FREE_IF_COPY(pgraster, 1);
2676 MemoryContextSwitchTo(oldcontext);
2677 elog(ERROR,
"RASTER_union_transfn: Could not run raster iterator function");
2685 for (
y = 0;
y < _dim[1];
y++) {
2686 POSTGIS_RT_DEBUGF(4,
"Getting pixel line of iterator raster at (x, y, length) = (0, %d, %d)",
y, _dim[0]);
2702 PG_FREE_IF_COPY(pgraster, 1);
2705 MemoryContextSwitchTo(oldcontext);
2706 elog(ERROR,
"RASTER_union_transfn: Could not get pixel line from band of working raster");
2710 POSTGIS_RT_DEBUGF(4,
"Setting pixel line at (x, y, length) = (%d, %d, %d)", (
int) _offset[2], (
int) _offset[3] +
y, nvals);
2713 (
int) _offset[2], (
int) _offset[3] +
y,
2725 PG_FREE_IF_COPY(pgraster, 1);
2728 MemoryContextSwitchTo(oldcontext);
2729 elog(ERROR,
"RASTER_union_transfn: Could not set pixel line to band of internal raster");
2749 hasnodata, nodataval,
2763 PG_FREE_IF_COPY(pgraster, 1);
2766 MemoryContextSwitchTo(oldcontext);
2767 elog(ERROR,
"RASTER_union_transfn: Could not run raster iterator function");
2786 PG_FREE_IF_COPY(pgraster, 1);
2790 MemoryContextSwitchTo(oldcontext);
2794 PG_RETURN_POINTER(iwr);
2814 double nodataval = 0;
2819 if (!AggCheckCallContext(fcinfo, NULL)) {
2820 elog(ERROR,
"RASTER_union_finalfn: Cannot be called in a non-aggregate context");
2825 if (PG_ARGISNULL(0))
2832 if (itrset == NULL) {
2834 elog(ERROR,
"RASTER_union_finalfn: Could not allocate memory for iterator arguments");
2838 for (i = 0; i < iwr->
numband; i++) {
2853 itrset[0].
nband = 0;
2855 itrset[1].
nband = 0;
2863 hasnodata, nodataval,
2876 hasnodata, nodataval,
2890 elog(ERROR,
"RASTER_union_finalfn: Could not run raster iterator function");
2896 if (_raster == NULL)
2902 uint32_t bandNums[1] = {0};
2904 status = (_rtn == NULL) ? -1 : 0;
2929 elog(ERROR,
"RASTER_union_finalfn: Could not add band to final raster");
2942 if (!_rtn) PG_RETURN_NULL();
2952 SET_VARSIZE(pgraster, pgraster->
size);
2953 PG_RETURN_POINTER(pgraster);
2982 elog(ERROR,
"rtpg_clip_arg_init: Could not allocate memory for function arguments");
2996 if (arg->
band != NULL)
3001 if (arg->
mask != NULL)
3009 double *
value,
int *nodata
3037 unsigned char *wkb = NULL;
3066 if (PG_ARGISNULL(0) || PG_ARGISNULL(2))
3072 elog(ERROR,
"RASTER_clip: Could not initialize argument structure");
3077 pgraster = (
rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
3081 if (arg->
raster == NULL) {
3083 PG_FREE_IF_COPY(pgraster, 0);
3084 elog(ERROR,
"RASTER_clip: Could not deserialize raster");
3090 elog(NOTICE,
"Input raster is empty or has no bands. Returning empty raster");
3093 PG_FREE_IF_COPY(pgraster, 0);
3097 elog(ERROR,
"RASTER_clip: Could not create empty raster");
3106 SET_VARSIZE(pgrtn, pgrtn->
size);
3107 PG_RETURN_POINTER(pgrtn);
3115 gser = PG_GETARG_GSERIALIZED_P(2);
3127 elog(NOTICE,
"Geometry provided does not have the same SRID as the raster. Returning NULL");
3130 PG_FREE_IF_COPY(pgraster, 0);
3132 PG_FREE_IF_COPY(gser, 2);
3138 if (!PG_ARGISNULL(4) && !PG_GETARG_BOOL(4))
3145 PG_FREE_IF_COPY(pgraster, 0);
3147 PG_FREE_IF_COPY(gser, 2);
3149 elog(ERROR,
"RASTER_clip: Could not get convex hull of raster");
3156 PG_FREE_IF_COPY(gser, 2);
3161 elog(NOTICE,
"The input raster and input geometry do not intersect. Returning empty raster");
3164 PG_FREE_IF_COPY(pgraster, 0);
3169 elog(ERROR,
"RASTER_clip: Could not create empty raster");
3178 SET_VARSIZE(pgrtn, pgrtn->
size);
3179 PG_RETURN_POINTER(pgrtn);
3183 if (!PG_ARGISNULL(1)) {
3184 array = PG_GETARG_ARRAYTYPE_P(1);
3185 etype = ARR_ELEMTYPE(array);
3186 get_typlenbyvalalign(etype, &typlen, &typbyval, &typalign);
3194 PG_FREE_IF_COPY(pgraster, 0);
3196 elog(ERROR,
"RASTER_clip: Invalid data type for band indexes");
3203 typlen, typbyval, typalign,
3208 if (arg->
band == NULL) {
3210 PG_FREE_IF_COPY(pgraster, 0);
3212 elog(ERROR,
"RASTER_clip: Could not allocate memory for band arguments");
3216 for (i = 0, j = 0; i < arg->
numbands; i++) {
3217 if (nulls[i])
continue;
3221 arg->
band[j].
nband = DatumGetInt16(e[i]) - 1;
3224 arg->
band[j].
nband = DatumGetInt32(e[i]) - 1;
3231 if (j < arg->numbands) {
3233 if (arg->
band == NULL) {
3235 PG_FREE_IF_COPY(pgraster, 0);
3237 elog(ERROR,
"RASTER_clip: Could not reallocate memory for band arguments");
3245 for (i = 0; i < arg->
numbands; i++) {
3247 elog(NOTICE,
"Band at index %d not found in raster", arg->
band[i].
nband + 1);
3249 PG_FREE_IF_COPY(pgraster, 0);
3264 if (arg->
band == NULL) {
3267 PG_FREE_IF_COPY(pgraster, 0);
3270 elog(ERROR,
"RASTER_clip: Could not allocate memory for band arguments");
3274 for (i = 0; i < arg->
numbands; i++) {
3283 if (!PG_ARGISNULL(3)) {
3284 array = PG_GETARG_ARRAYTYPE_P(3);
3285 etype = ARR_ELEMTYPE(array);
3286 get_typlenbyvalalign(etype, &typlen, &typbyval, &typalign);
3294 PG_FREE_IF_COPY(pgraster, 0);
3296 elog(ERROR,
"RASTER_clip: Invalid data type for NODATA values");
3303 typlen, typbyval, typalign,
3308 for (i = 0, j = 0; i < arg->
numbands; i++, j++) {
3349 if (arg->
mask == NULL) {
3351 PG_FREE_IF_COPY(pgraster, 0);
3352 elog(ERROR,
"RASTER_clip: Could not rasterize intersection geometry");
3363 if (itrset == NULL) {
3365 PG_FREE_IF_COPY(pgraster, 0);
3366 elog(ERROR,
"RASTER_clip: Could not allocate memory for iterator arguments");
3371 for (i = 0; i < arg->
numbands; i++) {
3372 POSTGIS_RT_DEBUGF(4,
"band arg %d (nband, hasnodata, nodataval) = (%d, %d, %f)",
3400 PG_FREE_IF_COPY(pgraster, 0);
3401 elog(ERROR,
"RASTER_clip: Could not create output raster");
3410 PG_FREE_IF_COPY(pgraster, 0);
3411 elog(ERROR,
"RASTER_clip: Could not add NODATA band to output raster");
3425 itrset[1].
nband = 0;
3433 hasnodata, nodataval,
3445 PG_FREE_IF_COPY(pgraster, 0);
3446 elog(ERROR,
"RASTER_clip: Could not run raster iterator function");
3461 PG_FREE_IF_COPY(pgraster, 0);
3462 elog(NOTICE,
"RASTER_clip: Could not get band from working raster");
3471 PG_FREE_IF_COPY(pgraster, 0);
3472 elog(ERROR,
"RASTER_clip: Could not add new band to output raster");
3482 PG_FREE_IF_COPY(pgraster, 0);
3492 SET_VARSIZE(pgrtn, pgrtn->
size);
3493 PG_RETURN_POINTER(pgrtn);
3506 uint32_t numBands = 0;
3526 HeapTupleHeader tup;
3531 text *exprtext = NULL;
3536 char *pixeltype = NULL;
3537 text *pixeltypetext = NULL;
3539 double nodataval = 0;
3540 bool hasnodata =
FALSE;
3542 char **comma_set = NULL;
3543 uint32_t comma_n = 0;
3544 char **colon_set = NULL;
3545 uint32_t colon_n = 0;
3546 char **dash_set = NULL;
3547 uint32_t dash_n = 0;
3552 if (PG_ARGISNULL(0))
3554 pgraster = (
rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
3559 PG_FREE_IF_COPY(pgraster, 0);
3560 elog(ERROR,
"RASTER_reclass: Could not deserialize raster");
3564 POSTGIS_RT_DEBUGF(3,
"RASTER_reclass: %d possible bands to be reclassified", numBands);
3568 array = PG_GETARG_ARRAYTYPE_P(1);
3569 etype = ARR_ELEMTYPE(array);
3570 get_typlenbyvalalign(etype, &typlen, &typbyval, &typalign);
3572 deconstruct_array(array, etype, typlen, typbyval, typalign, &e,
3576 elog(NOTICE,
"Invalid argument for reclassargset. Returning original raster");
3580 PG_FREE_IF_COPY(pgraster, 0);
3584 SET_VARSIZE(pgrtn, pgrtn->
size);
3585 PG_RETURN_POINTER(pgrtn);
3593 for (i = 0; i < n; i++) {
3594 if (nulls[i])
continue;
3597 tup = (HeapTupleHeader) DatumGetPointer(e[i]);
3599 elog(NOTICE,
"Invalid argument for reclassargset. Returning original raster");
3603 PG_FREE_IF_COPY(pgraster, 0);
3607 SET_VARSIZE(pgrtn, pgrtn->
size);
3608 PG_RETURN_POINTER(pgrtn);
3612 tupv = GetAttributeByName(tup,
"nband", &isnull);
3614 elog(NOTICE,
"Invalid argument for reclassargset. Missing value of nband for reclassarg of index %d . Returning original raster", i);
3618 PG_FREE_IF_COPY(pgraster, 0);
3622 SET_VARSIZE(pgrtn, pgrtn->
size);
3623 PG_RETURN_POINTER(pgrtn);
3625 nband = DatumGetInt32(tupv);
3629 if (nband < 1 || nband > numBands) {
3630 elog(NOTICE,
"Invalid argument for reclassargset. Invalid band index (must use 1-based) for reclassarg of index %d . Returning original raster", i);
3634 PG_FREE_IF_COPY(pgraster, 0);
3638 SET_VARSIZE(pgrtn, pgrtn->
size);
3639 PG_RETURN_POINTER(pgrtn);
3643 tupv = GetAttributeByName(tup,
"reclassexpr", &isnull);
3645 elog(NOTICE,
"Invalid argument for reclassargset. Missing value of reclassexpr for reclassarg of index %d . Returning original raster", i);
3649 PG_FREE_IF_COPY(pgraster, 0);
3653 SET_VARSIZE(pgrtn, pgrtn->
size);
3654 PG_RETURN_POINTER(pgrtn);
3656 exprtext = (text *) DatumGetPointer(tupv);
3657 if (NULL == exprtext) {
3658 elog(NOTICE,
"Invalid argument for reclassargset. Missing value of reclassexpr for reclassarg of index %d . Returning original raster", i);
3662 PG_FREE_IF_COPY(pgraster, 0);
3666 SET_VARSIZE(pgrtn, pgrtn->
size);
3667 PG_RETURN_POINTER(pgrtn);
3678 elog(NOTICE,
"Invalid argument for reclassargset. Invalid expression of reclassexpr for reclassarg of index %d . Returning original raster", i);
3682 PG_FREE_IF_COPY(pgraster, 0);
3686 SET_VARSIZE(pgrtn, pgrtn->
size);
3687 PG_RETURN_POINTER(pgrtn);
3694 for (a = 0, j = 0; a < comma_n; a++) {
3700 elog(NOTICE,
"Invalid argument for reclassargset. Invalid expression of reclassexpr for reclassarg of index %d . Returning original raster", i);
3701 for (k = 0; k < j; k++) pfree(exprset[k]);
3706 PG_FREE_IF_COPY(pgraster, 0);
3710 SET_VARSIZE(pgrtn, pgrtn->
size);
3711 PG_RETURN_POINTER(pgrtn);
3717 for (b = 0; b < colon_n; b++) {
3722 if (dash_n < 1 || dash_n > 3) {
3723 elog(NOTICE,
"Invalid argument for reclassargset. Invalid expression of reclassexpr for reclassarg of index %d . Returning original raster", i);
3724 for (k = 0; k < j; k++) pfree(exprset[k]);
3729 PG_FREE_IF_COPY(pgraster, 0);
3733 SET_VARSIZE(pgrtn, pgrtn->
size);
3734 PG_RETURN_POINTER(pgrtn);
3737 for (c = 0; c < dash_n; c++) {
3741 strlen(dash_set[c]) == 1 && (
3742 strchr(dash_set[c],
'(') != NULL ||
3743 strchr(dash_set[c],
'[') != NULL ||
3744 strchr(dash_set[c],
')') != NULL ||
3745 strchr(dash_set[c],
']') != NULL
3749 junk = palloc(
sizeof(
char) * (strlen(dash_set[c + 1]) + 2));
3751 for (k = 0; k <= j; k++) pfree(exprset[k]);
3754 PG_FREE_IF_COPY(pgraster, 0);
3756 elog(ERROR,
"RASTER_reclass: Could not allocate memory");
3760 sprintf(junk,
"%s%s", dash_set[c], dash_set[c + 1]);
3762 dash_set[c] = repalloc(dash_set[c],
sizeof(
char) * (strlen(junk) + 1));
3763 strcpy(dash_set[c], junk);
3767 for (dash_it = 1; dash_it < dash_n; dash_it++) {
3768 dash_set[dash_it - 1] = repalloc(dash_set[dash_it - 1], (strlen(dash_set[dash_it]) + 1) *
sizeof(
char));
3769 strcpy(dash_set[dash_it - 1], dash_set[dash_it]);
3773 pfree(dash_set[dash_n]);
3774 dash_set = repalloc(dash_set,
sizeof(
char *) * dash_n);
3778 if (c < 1 && dash_n > 2) {
3779 elog(NOTICE,
"Invalid argument for reclassargset. Invalid expression of reclassexpr for reclassarg of index %d . Returning original raster", i);
3780 for (k = 0; k < j; k++) pfree(exprset[k]);
3785 PG_FREE_IF_COPY(pgraster, 0);
3789 SET_VARSIZE(pgrtn, pgrtn->
size);
3790 PG_RETURN_POINTER(pgrtn);
3801 strchr(dash_set[c],
')') != NULL ||
3802 strchr(dash_set[c],
']') != NULL
3807 else if (strchr(dash_set[c],
'(') != NULL){
3817 strrchr(dash_set[c],
'(') != NULL ||
3818 strrchr(dash_set[c],
'[') != NULL
3823 else if (strrchr(dash_set[c],
']') != NULL) {
3831 POSTGIS_RT_DEBUGF(4,
"RASTER_reclass: exc_val %d inc_val %d", exc_val, inc_val);
3839 val = strtod(dash_set[c], &junk);
3840 if (errno != 0 || dash_set[c] == junk) {
3841 elog(NOTICE,
"Invalid argument for reclassargset. Invalid expression of reclassexpr for reclassarg of index %d . Returning original raster", i);
3842 for (k = 0; k < j; k++) pfree(exprset[k]);
3847 PG_FREE_IF_COPY(pgraster, 0);
3851 SET_VARSIZE(pgrtn, pgrtn->
size);
3852 PG_RETURN_POINTER(pgrtn);
3858 junk = strstr(colon_set[b], dash_set[c]);
3863 "(colon_set[%d], dash_set[%d], junk) = (%s, %s, %s)",
3864 b, c, colon_set[b], dash_set[c], junk
3867 if (junk != colon_set[b]) {
3869 if (*(junk - 1) ==
'-') {
3872 ((junk - 1) == colon_set[b]) ||
3873 (*(junk - 2) ==
'-') ||
3874 (*(junk - 2) ==
'[') ||
3875 (*(junk - 2) ==
'(')
3895 exprset[j]->
src.
min = val;
3901 exprset[j]->
src.
max = val;
3911 exprset[j]->
dst.
min = val;
3914 exprset[j]->
dst.
max = val;
3922 , exprset[j]->src.min
3932 tupv = GetAttributeByName(tup,
"pixeltype", &isnull);
3934 elog(NOTICE,
"Invalid argument for reclassargset. Missing value of pixeltype for reclassarg of index %d . Returning original raster", i);
3938 PG_FREE_IF_COPY(pgraster, 0);
3942 SET_VARSIZE(pgrtn, pgrtn->
size);
3943 PG_RETURN_POINTER(pgrtn);
3945 pixeltypetext = (text *) DatumGetPointer(tupv);
3946 if (NULL == pixeltypetext) {
3947 elog(NOTICE,
"Invalid argument for reclassargset. Missing value of pixeltype for reclassarg of index %d . Returning original raster", i);
3951 PG_FREE_IF_COPY(pgraster, 0);
3955 SET_VARSIZE(pgrtn, pgrtn->
size);
3956 PG_RETURN_POINTER(pgrtn);
3963 tupv = GetAttributeByName(tup,
"nodataval", &isnull);
3969 nodataval = DatumGetFloat8(tupv);
3978 elog(NOTICE,
"Could not find raster band of index %d. Returning original raster",
nband);
3979 for (k = 0; k < j; k++) pfree(exprset[k]);
3984 PG_FREE_IF_COPY(pgraster, 0);
3988 SET_VARSIZE(pgrtn, pgrtn->
size);
3989 PG_RETURN_POINTER(pgrtn);
3993 for (k = 0; k < j; k++) pfree(exprset[k]);
3997 PG_FREE_IF_COPY(pgraster, 0);
3999 elog(ERROR,
"RASTER_reclass: Could not reclassify raster band of index %d",
nband);
4005 for (k = 0; k < j; k++) pfree(exprset[k]);
4010 PG_FREE_IF_COPY(pgraster, 0);
4012 elog(ERROR,
"RASTER_reclass: Could not replace raster band of index %d with reclassified band",
nband);
4020 for (k = 0; k < j; k++) pfree(exprset[k]);
4026 PG_FREE_IF_COPY(pgraster, 0);
4032 SET_VARSIZE(pgrtn, pgrtn->
size);
4033 PG_RETURN_POINTER(pgrtn);
4062 elog(ERROR,
"rtpg_colormap_arg: Could not allocate memory for function arguments");
4073 elog(ERROR,
"rtpg_colormap_arg: Could not allocate memory for function arguments");
4106 for (i = 0; i < arg->
nentry; i++) {
4107 if (arg->
entry[i] != NULL)
4108 pfree(arg->
entry[i]);
4114 for (i = 0; i < arg->
nelement; i++)
4134 if (PG_ARGISNULL(0))
4140 elog(ERROR,
"RASTER_colorMap: Could not initialize argument structure");
4145 pgraster = (
rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
4151 PG_FREE_IF_COPY(pgraster, 0);
4152 elog(ERROR,
"RASTER_colorMap: Could not deserialize raster");
4157 if (!PG_ARGISNULL(1))
4158 arg->
nband = PG_GETARG_INT32(1);
4163 elog(NOTICE,
"Raster does not have band at index %d. Returning empty raster", arg->
nband);
4168 PG_FREE_IF_COPY(pgraster, 0);
4169 elog(ERROR,
"RASTER_colorMap: Could not create empty raster");
4174 PG_FREE_IF_COPY(pgraster, 0);
4178 if (pgraster == NULL)
4181 SET_VARSIZE(pgraster, ((
rt_pgraster*) pgraster)->size);
4182 PG_RETURN_POINTER(pgraster);
4187 if (arg->
band == NULL) {
4190 PG_FREE_IF_COPY(pgraster, 0);
4191 elog(ERROR,
"RASTER_colorMap: Could not get band at index %d",
nband);
4196 if (!PG_ARGISNULL(3)) {
4197 char *method = NULL;
4205 if (strcmp(method,
"INTERPOLATE") == 0)
4207 else if (strcmp(method,
"EXACT") == 0)
4209 else if (strcmp(method,
"NEAREST") == 0)
4212 elog(NOTICE,
"Unknown value provided for method. Defaulting to INTERPOLATE");
4222 if (PG_ARGISNULL(2)) {
4224 PG_FREE_IF_COPY(pgraster, 0);
4225 elog(ERROR,
"RASTER_colorMap: Value must be provided for colormap");
4239 if (!strlen(colormap)) {
4241 PG_FREE_IF_COPY(pgraster, 0);
4242 elog(ERROR,
"RASTER_colorMap: Value must be provided for colormap");
4250 PG_FREE_IF_COPY(pgraster, 0);
4251 elog(ERROR,
"RASTER_colorMap: Could not process the value provided for colormap");
4259 PG_FREE_IF_COPY(pgraster, 0);
4260 elog(ERROR,
"RASTER_colorMap: Could not allocate memory for colormap entries");
4266 for (i = 0; i < arg->
nentry; i++) {
4280 if (!strlen(_entry)) {
4290 PG_FREE_IF_COPY(pgraster, 0);
4291 elog(ERROR,
"RASTER_colorMap: Could not process colormap entry %d", i + 1);
4295 elog(NOTICE,
"More than five elements in colormap entry %d. Using at most five elements", i + 1);
4304 for (j = 0; j < arg->
nelement; j++) {
4313 char *percent = NULL;
4317 strcmp(_element,
"NV") == 0 ||
4318 strcmp(_element,
"NULL") == 0 ||
4319 strcmp(_element,
"NODATA") == 0
4324 elog(NOTICE,
"More than one NODATA entry found. Using only the first one");
4333 else if ((percent = strchr(_element,
'%')) != NULL) {
4345 PG_FREE_IF_COPY(pgraster, 0);
4346 elog(ERROR,
"RASTER_colorMap: Could not get band's summary stats to process percentages");
4352 tmp = palloc(
sizeof(
char) * (percent - _element + 1));
4356 PG_FREE_IF_COPY(pgraster, 0);
4357 elog(ERROR,
"RASTER_colorMap: Could not allocate memory for value of percentage");
4361 memcpy(tmp, _element, percent - _element);
4362 tmp[percent - _element] =
'\0';
4367 value = strtod(tmp, NULL);
4369 if (errno != 0 || _element == junk) {
4372 PG_FREE_IF_COPY(pgraster, 0);
4373 elog(ERROR,
"RASTER_colorMap: Could not process percent string to value");
4379 elog(NOTICE,
"Percentage values cannot be less than zero. Defaulting to zero");
4382 else if (
value > 100.) {
4383 elog(NOTICE,
"Percentage values cannot be greater than 100. Defaulting to 100");
4395 if (errno != 0 || _element == junk) {
4398 PG_FREE_IF_COPY(pgraster, 0);
4399 elog(ERROR,
"RASTER_colorMap: Could not process string to value");
4410 value = (int) strtod(_element, &junk);
4411 if (errno != 0 || _element == junk) {
4414 PG_FREE_IF_COPY(pgraster, 0);
4415 elog(ERROR,
"RASTER_colorMap: Could not process string to value");
4420 elog(NOTICE,
"RGBA value cannot be greater than 255. Defaulting to 255");
4423 else if (
value < 0) {
4424 elog(NOTICE,
"RGBA value cannot be less than zero. Defaulting to zero");
4433 POSTGIS_RT_DEBUGF(4,
"colormap->entry[%d] (isnodata, value, R, G, B, A) = (%d, %f, %d, %d, %d, %d)",
4454 PG_FREE_IF_COPY(pgraster, 0);
4455 elog(ERROR,
"RASTER_colorMap: Could not create new raster with applied colormap");
4460 PG_FREE_IF_COPY(pgraster, 0);
4466 if (pgraster == NULL)
4469 SET_VARSIZE(pgraster, ((
rt_pgraster*) pgraster)->size);
4470 PG_RETURN_POINTER(pgraster);
4482 int x,
y,
nband, width, height;
4484 double newnodatavalue = 0.0;
4485 double newinitialvalue = 0.0;
4486 double newval = 0.0;
4487 char *newexpr = NULL;
4488 char *initexpr = NULL;
4489 char *expression = NULL;
4490 int hasnodataval = 0;
4491 double nodataval = 0.;
4493 int skipcomputation = 0;
4495 const int argkwcount = 3;
4496 enum KEYWORDS { kVAL=0, kX=1, kY=2 };
4497 char *argkw[] = {
"[rast]",
"[rast.x]",
"[rast.y]"};
4498 Oid argkwtypes[] = { FLOAT8OID, INT4OID, INT4OID };
4500 Oid argtype[] = { FLOAT8OID, INT4OID, INT4OID };
4501 uint8_t argpos[3] = {0};
4506 SPIPlanPtr spi_plan = NULL;
4507 SPITupleTable * tuptable = NULL;
4509 char * strFromText = NULL;
4510 Datum *values = NULL;
4511 Datum datum = (Datum)NULL;
4513 bool isnull =
FALSE;
4520 if (PG_ARGISNULL(0)) {
4521 elog(NOTICE,
"Raster is NULL. Returning NULL");
4527 pgraster = (
rt_pgraster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
4530 PG_FREE_IF_COPY(pgraster, 0);
4531 elog(ERROR,
"RASTER_mapAlgebraExpr: Could not deserialize raster");
4537 if (PG_ARGISNULL(1))
4540 nband = PG_GETARG_INT32(1);
4546 POSTGIS_RT_DEBUG(3,
"RASTER_mapAlgebraExpr: Creating new empty raster...");
4557 if ( NULL == newrast ) {
4558 PG_FREE_IF_COPY(pgraster, 0);
4559 elog(ERROR,
"RASTER_mapAlgebraExpr: Could not create a new raster");
4584 elog(NOTICE,
"Raster is empty. Returning an empty raster");
4586 PG_FREE_IF_COPY(pgraster, 0);
4590 if (NULL == pgrtn) {
4592 elog(ERROR,
"RASTER_mapAlgebraExpr: Could not serialize raster");
4596 SET_VARSIZE(pgrtn, pgrtn->
size);
4597 PG_RETURN_POINTER(pgrtn);
4608 elog(NOTICE,
"Raster does not have the required band. Returning a raster "
4611 PG_FREE_IF_COPY(pgraster, 0);
4615 if (NULL == pgrtn) {
4616 elog(ERROR,
"RASTER_mapAlgebraExpr: Could not serialize raster");
4620 SET_VARSIZE(pgrtn, pgrtn->
size);
4621 PG_RETURN_POINTER(pgrtn);
4626 if ( NULL ==
band ) {
4627 elog(NOTICE,
"Could not get the required band. Returning a raster "
4630 PG_FREE_IF_COPY(pgraster, 0);
4634 if (NULL == pgrtn) {
4635 elog(ERROR,
"RASTER_mapAlgebraExpr: Could not serialize raster");
4639 SET_VARSIZE(pgrtn, pgrtn->
size);
4640 PG_RETURN_POINTER(pgrtn);
4646 POSTGIS_RT_DEBUG(3,
"RASTER_mapAlgebraExpr: Getting NODATA value for band...");
4664 newinitialvalue = newnodatavalue;
4671 if (PG_ARGISNULL(2)) {
4679 if (newpixeltype ==
PT_END)
4683 if (newpixeltype ==
PT_END) {
4684 PG_FREE_IF_COPY(pgraster, 0);
4685 elog(ERROR,
"RASTER_mapAlgebraExpr: Invalid pixeltype");
4694 if (!PG_ARGISNULL(3)) {
4696 len = strlen(
"SELECT (") + strlen(expression) + strlen(
")::double precision");
4697 initexpr = (
char *)palloc(len + 1);
4699 memcpy(initexpr,
"SELECT (", strlen(
"SELECT ("));
4700 memcpy(initexpr + strlen(
"SELECT ("), expression, strlen(expression));
4701 memcpy(initexpr + strlen(
"SELECT (") + strlen(expression),
")::double precision", strlen(
")::double precision"));
4702 initexpr[len] =
'\0';
4720 if (!PG_ARGISNULL(4)) {
4722 nodataval = PG_GETARG_FLOAT8(4);
4723 newinitialvalue = nodataval;
4740 POSTGIS_RT_DEBUG(3,
"RASTER_mapAlgebraExpr: Band is a nodata band, returning "
4741 "a raster filled with nodata");
4744 newinitialvalue,
TRUE, newnodatavalue, 0);
4750 PG_FREE_IF_COPY(pgraster, 0);
4755 if (NULL == pgrtn) {
4756 elog(ERROR,
"RASTER_mapAlgebraExpr: Could not serialize raster");
4760 SET_VARSIZE(pgrtn, pgrtn->
size);
4761 PG_RETURN_POINTER(pgrtn);
4769 if (initexpr != NULL && ( !strcmp(initexpr,
"SELECT [rast]") || !strcmp(initexpr,
"SELECT [rast.val]") ) && !hasnodataval) {
4772 "Returning raster with band %d from original raster",
nband);
4784 PG_FREE_IF_COPY(pgraster, 0);
4789 if (NULL == pgrtn) {
4790 elog(ERROR,
"RASTER_mapAlgebraExpr: Could not serialize raster");
4794 SET_VARSIZE(pgrtn, pgrtn->
size);
4795 PG_RETURN_POINTER(pgrtn);
4802 if (initexpr != NULL && strstr(initexpr,
"[rast") == NULL) {
4803 ret = SPI_connect();
4804 if (ret != SPI_OK_CONNECT) {
4805 PG_FREE_IF_COPY(pgraster, 0);
4806 elog(ERROR,
"RASTER_mapAlgebraExpr: Could not connect to the SPI manager");
4811 ret = SPI_execute(initexpr,
FALSE, 0);
4813 if (ret != SPI_OK_SELECT || SPI_tuptable == NULL || SPI_processed != 1) {
4817 SPI_freetuptable(tuptable);
4818 PG_FREE_IF_COPY(pgraster, 0);
4821 elog(ERROR,
"RASTER_mapAlgebraExpr: Invalid construction for expression");
4825 tupdesc = SPI_tuptable->tupdesc;
4826 tuptable = SPI_tuptable;
4828 tuple = tuptable->vals[0];
4829 newexpr = SPI_getvalue(tuple, tupdesc, 1);
4831 POSTGIS_RT_DEBUG(3,
"Constant expression evaluated to NULL, keeping initvalue");
4832 newval = newinitialvalue;
4834 newval = atof(newexpr);
4837 SPI_freetuptable(tuptable);
4844 skipcomputation = 1;
4850 if (!hasnodataval) {
4851 newinitialvalue = newval;
4852 skipcomputation = 2;
4856 else if (
FLT_NEQ(newval, newinitialvalue)) {
4857 skipcomputation = 2;
4866 newinitialvalue,
TRUE, newnodatavalue, 0);
4873 if (expression == NULL || skipcomputation == 2) {
4879 PG_FREE_IF_COPY(pgraster, 0);
4884 if (NULL == pgrtn) {
4885 elog(ERROR,
"RASTER_mapAlgebraExpr: Could not serialize raster");
4889 SET_VARSIZE(pgrtn, pgrtn->
size);
4890 PG_RETURN_POINTER(pgrtn);
4893 RASTER_DEBUG(3,
"RASTER_mapAlgebraExpr: Creating new raster band...");
4897 if ( NULL == newband ) {
4898 elog(NOTICE,
"Could not modify band for new raster. Returning new "
4899 "raster with the original band");
4904 PG_FREE_IF_COPY(pgraster, 0);
4909 if (NULL == pgrtn) {
4910 elog(ERROR,
"RASTER_mapAlgebraExpr: Could not serialize raster");
4914 SET_VARSIZE(pgrtn, pgrtn->
size);
4915 PG_RETURN_POINTER(pgrtn);
4921 if (initexpr != NULL) {
4924 pfree(initexpr); initexpr=newexpr;
4926 sprintf(place,
"$1");
4927 for (i = 0, j = 1; i < argkwcount; i++) {
4930 pfree(initexpr); initexpr=newexpr;
4932 argtype[argcount] = argkwtypes[i];
4936 sprintf(place,
"$%d", j);
4946 values = (Datum *) palloc(
sizeof(Datum) * argcount);
4947 if (values == NULL) {
4952 PG_FREE_IF_COPY(pgraster, 0);
4955 elog(ERROR,
"RASTER_mapAlgebraExpr: Could not allocate memory for value parameters of prepared statement");
4960 nulls = (
char *)palloc(argcount);
4961 if (nulls == NULL) {
4966 PG_FREE_IF_COPY(pgraster, 0);
4969 elog(ERROR,
"RASTER_mapAlgebraExpr: Could not allocate memory for null parameters of prepared statement");
4974 ret = SPI_connect();
4975 if (ret != SPI_OK_CONNECT) {
4980 PG_FREE_IF_COPY(pgraster, 0);
4983 elog(ERROR,
"RASTER_mapAlgebraExpr: Could not connect to the SPI manager");
4988 spi_plan = SPI_prepare(initexpr, argcount, argtype);
4990 if (spi_plan == NULL) {
4993 PG_FREE_IF_COPY(pgraster, 0);
5000 elog(ERROR,
"RASTER_mapAlgebraExpr: Could not prepare expression");
5005 for (
x = 0;
x < width;
x++) {
5006 for(
y = 0;
y < height;
y++) {
5014 if (skipcomputation == 0) {
5015 if (initexpr != NULL) {
5017 memset(nulls,
'n', argcount);
5019 for (i = 0; i < argkwcount; i++) {
5021 if (idx < 1)
continue;
5026 values[idx] = Int32GetDatum(
x+1);
5031 values[idx] = Int32GetDatum(
y+1);
5034 else if (i == kVAL ) {
5035 values[idx] = Float8GetDatum(
r);
5041 ret = SPI_execute_plan(spi_plan, values, nulls,
FALSE, 0);
5042 if (ret != SPI_OK_SELECT || SPI_tuptable == NULL ||
5043 SPI_processed != 1) {
5045 SPI_freetuptable(tuptable);
5047 SPI_freeplan(spi_plan);
5055 PG_FREE_IF_COPY(pgraster, 0);
5058 elog(ERROR,
"RASTER_mapAlgebraExpr: Error executing prepared plan");
5063 tupdesc = SPI_tuptable->tupdesc;
5064 tuptable = SPI_tuptable;
5066 tuple = tuptable->vals[0];
5067 datum = SPI_getbinval(tuple, tupdesc, 1, &isnull);
5068 if ( SPI_result == SPI_ERROR_NOATTRIBUTE ) {
5069 POSTGIS_RT_DEBUGF(3,
"Expression for pixel %d,%d (value %g) errored, skip setting",
x+1,
y+1,
r);
5070 newval = newinitialvalue;
5072 else if ( isnull ) {
5073 POSTGIS_RT_DEBUGF(3,
"Expression for pixel %d,%d (value %g) evaluated to NULL, skip setting",
x+1,
y+1,
r);
5074 newval = newinitialvalue;
5076 newval = DatumGetFloat8(datum);
5079 SPI_freetuptable(tuptable);
5083 newval = newinitialvalue;
5096 if (initexpr != NULL) {
5097 SPI_freeplan(spi_plan);
5111 POSTGIS_RT_DEBUG(3,
"RASTER_mapAlgebraExpr: raster modified, serializing it.");
5115 PG_FREE_IF_COPY(pgraster, 0);
5122 SET_VARSIZE(pgrtn, pgrtn->
size);
5130 PG_RETURN_POINTER(pgrtn);
5145 int x,
y,
nband, width, height;
5147 double newnodatavalue = 0.0;
5148 double newinitialvalue = 0.0;
5149 double newval = 0.0;
5154 #if POSTGIS_PGSQL_VERSION < 120
5155 FunctionCallInfoData cbdata;
5157 LOCAL_FCINFO(cbdata, FUNC_MAX_ARGS);
5160 char * strFromText = NULL;
5166 if (PG_ARGISNULL(0)) {
5167 elog(WARNING,
"Raster is NULL. Returning NULL");
5173 pgraster = (
rt_pgraster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
5176 PG_FREE_IF_COPY(pgraster, 0);
5177 elog(ERROR,
"RASTER_mapAlgebraFct: Could not deserialize raster");
5185 if (PG_ARGISNULL(1))
5188 nband = PG_GETARG_INT32(1);
5204 if ( NULL == newrast ) {
5207 PG_FREE_IF_COPY(pgraster, 0);
5209 elog(ERROR,
"RASTER_mapAlgebraFct: Could not create a new raster");
5234 elog(NOTICE,
"Raster is empty. Returning an empty raster");
5236 PG_FREE_IF_COPY(pgraster, 0);
5240 if (NULL == pgrtn) {
5241 elog(ERROR,
"RASTER_mapAlgebraFct: Could not serialize raster");
5245 SET_VARSIZE(pgrtn, pgrtn->
size);
5246 PG_RETURN_POINTER(pgrtn);
5256 elog(NOTICE,
"Raster does not have the required band. Returning a raster "
5259 PG_FREE_IF_COPY(pgraster, 0);
5263 if (NULL == pgrtn) {
5264 elog(ERROR,
"RASTER_mapAlgebraFct: Could not serialize raster");
5268 SET_VARSIZE(pgrtn, pgrtn->
size);
5269 PG_RETURN_POINTER(pgrtn);
5274 if ( NULL ==
band ) {
5275 elog(NOTICE,
"Could not get the required band. Returning a raster "
5278 PG_FREE_IF_COPY(pgraster, 0);
5282 if (NULL == pgrtn) {
5283 elog(ERROR,
"RASTER_mapAlgebraFct: Could not serialize raster");
5287 SET_VARSIZE(pgrtn, pgrtn->
size);
5288 PG_RETURN_POINTER(pgrtn);
5294 POSTGIS_RT_DEBUG(3,
"RASTER_mapAlgebraFct: Getting NODATA value for band...");
5311 newinitialvalue = newnodatavalue;
5318 if (PG_ARGISNULL(2)) {
5326 if (newpixeltype ==
PT_END)
5330 if (newpixeltype ==
PT_END) {
5333 PG_FREE_IF_COPY(pgraster, 0);
5336 elog(ERROR,
"RASTER_mapAlgebraFct: Invalid pixeltype");
5344 if (PG_ARGISNULL(3)) {
5347 PG_FREE_IF_COPY(pgraster, 0);
5350 elog(ERROR,
"RASTER_mapAlgebraFct: Required function is missing. Returning NULL");
5354 oid = PG_GETARG_OID(3);
5355 if (oid == InvalidOid) {
5358 PG_FREE_IF_COPY(pgraster, 0);
5361 elog(ERROR,
"RASTER_mapAlgebraFct: Got invalid function object id. Returning NULL");
5365 fmgr_info(oid, &cbinfo);
5368 if (cbinfo.fn_retset) {
5371 PG_FREE_IF_COPY(pgraster, 0);
5374 elog(ERROR,
"RASTER_mapAlgebraFct: Function provided must return double precision not resultset");
5378 else if (cbinfo.fn_nargs < 2 || cbinfo.fn_nargs > 3) {
5381 PG_FREE_IF_COPY(pgraster, 0);
5384 elog(ERROR,
"RASTER_mapAlgebraFct: Function does not have two or three input parameters");
5388 if (cbinfo.fn_nargs == 2)
5393 if (func_volatile(oid) ==
'v') {
5394 elog(NOTICE,
"Function provided is VOLATILE. Unless required and for best performance, function should be IMMUTABLE or STABLE");
5398 #if POSTGIS_PGSQL_VERSION < 120
5399 InitFunctionCallInfoData(cbdata, &cbinfo, 2, InvalidOid, NULL, NULL);
5401 cbdata.argnull[0] =
FALSE;
5402 cbdata.argnull[1] =
FALSE;
5403 cbdata.argnull[2] =
FALSE;
5405 InitFunctionCallInfoData(*cbdata, &cbinfo, 2, InvalidOid, NULL, NULL);
5407 cbdata->args[0].isnull =
FALSE;
5408 cbdata->args[1].isnull =
FALSE;
5409 cbdata->args[2].isnull =
FALSE;
5413 if (PG_ARGISNULL(4)) {
5414 if (cbinfo.fn_strict) {
5417 PG_FREE_IF_COPY(pgraster, 0);
5420 elog(ERROR,
"RASTER_mapAlgebraFct: Strict callback functions cannot have null parameters");
5424 #if POSTGIS_PGSQL_VERSION < 120
5425 cbdata.arg[k] = (Datum)NULL;
5426 cbdata.argnull[k] =
TRUE;
5428 cbdata->args[k].value = (Datum)NULL;
5429 cbdata->args[k].isnull =
TRUE;
5433 #if POSTGIS_PGSQL_VERSION < 120
5434 cbdata.arg[k] = PG_GETARG_DATUM(4);
5436 cbdata->args[k].value = PG_GETARG_DATUM(4);
5447 POSTGIS_RT_DEBUG(3,
"RASTER_mapAlgebraFct: Band is a nodata band, returning "
5448 "a raster filled with nodata");
5451 newinitialvalue,
TRUE, newnodatavalue, 0);
5454 PG_FREE_IF_COPY(pgraster, 0);
5459 if (NULL == pgrtn) {
5460 elog(ERROR,
"RASTER_mapAlgebraFct: Could not serialize raster");
5464 SET_VARSIZE(pgrtn, pgrtn->
size);
5465 PG_RETURN_POINTER(pgrtn);
5474 newinitialvalue,
TRUE, newnodatavalue, 0);
5478 if ( NULL == newband ) {
5479 elog(NOTICE,
"Could not modify band for new raster. Returning new "
5480 "raster with the original band");
5483 PG_FREE_IF_COPY(pgraster, 0);
5488 if (NULL == pgrtn) {
5489 elog(ERROR,
"RASTER_mapAlgebraFct: Could not serialize raster");
5493 SET_VARSIZE(pgrtn, pgrtn->
size);
5494 PG_RETURN_POINTER(pgrtn);
5500 for (
x = 0;
x < width;
x++) {
5501 for(
y = 0;
y < height;
y++) {
5509 if (
FLT_EQ(
r, newnodatavalue)) {
5510 if (cbinfo.fn_strict) {
5511 POSTGIS_RT_DEBUG(3,
"RASTER_mapAlgebraFct: Strict callbacks cannot accept NULL arguments, skipping NODATA cell.");
5514 #if POSTGIS_PGSQL_VERSION < 120
5515 cbdata.argnull[0] =
TRUE;
5516 cbdata.arg[0] = (Datum)NULL;
5518 cbdata->args[0].isnull =
TRUE;
5519 cbdata->args[0].value = (Datum)NULL;
5523 #if POSTGIS_PGSQL_VERSION < 120
5524 cbdata.argnull[0] =
FALSE;
5525 cbdata.arg[0] = Float8GetDatum(
r);
5527 cbdata->args[0].isnull =
FALSE;
5528 cbdata->args[0].value = Float8GetDatum(
r);
5533 if (cbinfo.fn_nargs == 3) {
5537 d[0] = Int32GetDatum(
x+1);
5538 d[1] = Int32GetDatum(
y+1);
5540 a = construct_array(d, 2, INT4OID,
sizeof(
int32),
true,
'i');
5542 #if POSTGIS_PGSQL_VERSION < 120
5543 cbdata.argnull[1] =
FALSE;
5544 cbdata.arg[1] = PointerGetDatum(a);
5546 cbdata->args[1].isnull =
FALSE;
5547 cbdata->args[1].value = PointerGetDatum(a);
5554 #if POSTGIS_PGSQL_VERSION < 120
5555 tmpnewval = FunctionCallInvoke(&cbdata);
5557 if (cbdata.isnull) {
5558 newval = newnodatavalue;
5561 tmpnewval = FunctionCallInvoke(cbdata);
5565 newval = newnodatavalue;
5569 newval = DatumGetFloat8(tmpnewval);
5583 POSTGIS_RT_DEBUG(3,
"RASTER_mapAlgebraFct: raster modified, serializing it.");
5587 PG_FREE_IF_COPY(pgraster, 0);
5598 SET_VARSIZE(pgrtn, pgrtn->
size);
5599 PG_RETURN_POINTER(pgrtn);
5614 int x,
y,
nband, width, height, ngbwidth, ngbheight, winwidth, winheight, u, v, nIndex, nNullItems;
5616 double newnodatavalue = 0.0;
5617 double newinitialvalue = 0.0;
5618 double newval = 0.0;
5623 #if POSTGIS_PGSQL_VERSION < 120
5624 FunctionCallInfoData cbdata;
5626 LOCAL_FCINFO(cbdata, FUNC_MAX_ARGS);
5629 ArrayType * neighborDatum;
5630 char * strFromText = NULL;
5631 text * txtNodataMode = NULL;
5632 text * txtCallbackParam = NULL;
5634 float fltReplace = 0;
5635 bool valuereplace =
false, pixelreplace, nNodataOnly =
true, nNullSkip =
false;
5636 Datum * neighborData = NULL;
5637 bool * neighborNulls = NULL;
5638 int neighborDims[2];
5647 if (PG_ARGISNULL(0)) {
5648 elog(WARNING,
"Raster is NULL. Returning NULL");
5654 pgraster = (
rt_pgraster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
5658 PG_FREE_IF_COPY(pgraster, 0);
5659 elog(ERROR,
"RASTER_mapAlgebraFctNgb: Could not deserialize raster");
5667 if (PG_ARGISNULL(1))
5670 nband = PG_GETARG_INT32(1);
5675 POSTGIS_RT_DEBUG(3,
"RASTER_mapAlgebraFctNgb: Creating new empty raster...");
5686 if ( NULL == newrast ) {
5688 PG_FREE_IF_COPY(pgraster, 0);
5689 elog(ERROR,
"RASTER_mapAlgebraFctNgb: Could not create a new raster");
5714 elog(NOTICE,
"Raster is empty. Returning an empty raster");
5716 PG_FREE_IF_COPY(pgraster, 0);
5720 if (NULL == pgrtn) {
5721 elog(ERROR,
"RASTER_mapAlgebraFctNgb: Could not serialize raster");
5725 SET_VARSIZE(pgrtn, pgrtn->
size);
5726 PG_RETURN_POINTER(pgrtn);
5736 elog(NOTICE,
"Raster does not have the required band. Returning a raster "
5739 PG_FREE_IF_COPY(pgraster, 0);
5743 if (NULL == pgrtn) {
5744 elog(ERROR,
"RASTER_mapAlgebraFctNgb: Could not serialize raster");
5748 SET_VARSIZE(pgrtn, pgrtn->
size);
5749 PG_RETURN_POINTER(pgrtn);
5754 if ( NULL ==
band ) {
5755 elog(NOTICE,
"Could not get the required band. Returning a raster "
5758 PG_FREE_IF_COPY(pgraster, 0);
5762 if (NULL == pgrtn) {
5763 elog(ERROR,
"RASTER_mapAlgebraFctNgb: Could not serialize raster");
5767 SET_VARSIZE(pgrtn, pgrtn->
size);
5768 PG_RETURN_POINTER(pgrtn);
5774 POSTGIS_RT_DEBUG(3,
"RASTER_mapAlgebraFctNgb: Getting NODATA value for band...");
5791 newinitialvalue = newnodatavalue;
5798 if (PG_ARGISNULL(2)) {
5804 POSTGIS_RT_DEBUGF(3,
"RASTER_mapAlgebraFctNgb: Pixeltype parameter: %s", strFromText);
5807 if (newpixeltype ==
PT_END)
5811 if (newpixeltype ==
PT_END) {
5814 PG_FREE_IF_COPY(pgraster, 0);
5817 elog(ERROR,
"RASTER_mapAlgebraFctNgb: Invalid pixeltype");
5825 if (PG_ARGISNULL(5)) {
5828 PG_FREE_IF_COPY(pgraster, 0);
5831 elog(ERROR,
"RASTER_mapAlgebraFctNgb: Required function is missing");
5835 oid = PG_GETARG_OID(5);
5836 if (oid == InvalidOid) {
5839 PG_FREE_IF_COPY(pgraster, 0);
5842 elog(ERROR,
"RASTER_mapAlgebraFctNgb: Got invalid function object id");
5846 fmgr_info(oid, &cbinfo);
5849 if (cbinfo.fn_retset) {
5852 PG_FREE_IF_COPY(pgraster, 0);
5855 elog(ERROR,
"RASTER_mapAlgebraFctNgb: Function provided must return double precision not resultset");
5859 else if (cbinfo.fn_nargs != 3) {
5862 PG_FREE_IF_COPY(pgraster, 0);
5865 elog(ERROR,
"RASTER_mapAlgebraFctNgb: Function does not have three input parameters");
5869 if (func_volatile(oid) ==
'v') {
5870 elog(NOTICE,
"Function provided is VOLATILE. Unless required and for best performance, function should be IMMUTABLE or STABLE");
5874 #if POSTGIS_PGSQL_VERSION < 120
5875 InitFunctionCallInfoData(cbdata, &cbinfo, 3, InvalidOid, NULL, NULL);
5876 memset(cbdata.argnull,
FALSE,
sizeof(
bool) * 3);
5878 InitFunctionCallInfoData(*cbdata, &cbinfo, 3, InvalidOid, NULL, NULL);
5879 cbdata->args[0].isnull =
FALSE;
5880 cbdata->args[1].isnull =
FALSE;
5881 cbdata->args[2].isnull =
FALSE;
5885 if (PG_ARGISNULL(7)) {
5886 if (cbinfo.fn_strict) {
5889 PG_FREE_IF_COPY(pgraster, 0);
5892 elog(ERROR,
"RASTER_mapAlgebraFctNgb: Strict callback functions cannot have NULL parameters");
5896 #if POSTGIS_PGSQL_VERSION < 120
5897 cbdata.arg[2] = (Datum)NULL;
5898 cbdata.argnull[2] =
TRUE;
5900 cbdata->args[2].value = (Datum)NULL;
5901 cbdata->args[2].isnull =
TRUE;
5905 #if POSTGIS_PGSQL_VERSION < 120
5906 cbdata.arg[2] = PG_GETARG_DATUM(7);
5908 cbdata->args[2].value = PG_GETARG_DATUM(7);
5919 POSTGIS_RT_DEBUG(3,
"RASTER_mapAlgebraFctNgb: Band is a nodata band, returning "
5920 "a raster filled with nodata");
5923 newinitialvalue,
TRUE, newnodatavalue, 0);
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);
5946 newinitialvalue,
TRUE, newnodatavalue, 0);
5950 if ( NULL == newband ) {
5951 elog(NOTICE,
"Could not modify band for new raster. Returning new "
5952 "raster with the original band");
5955 PG_FREE_IF_COPY(pgraster, 0);
5960 if (NULL == pgrtn) {
5961 elog(ERROR,
"RASTER_mapAlgebraFctNgb: Could not serialize raster");
5965 SET_VARSIZE(pgrtn, pgrtn->
size);
5966 PG_RETURN_POINTER(pgrtn);
5970 if (PG_ARGISNULL(3) || PG_GETARG_INT32(3) <= 0) {
5971 elog(NOTICE,
"Neighborhood width is NULL or <= 0. Returning new "
5972 "raster with the original band");
5975 PG_FREE_IF_COPY(pgraster, 0);
5980 if (NULL == pgrtn) {
5981 elog(ERROR,
"RASTER_mapAlgebraFctNgb: Could not serialize raster");
5985 SET_VARSIZE(pgrtn, pgrtn->
size);
5986 PG_RETURN_POINTER(pgrtn);
5989 ngbwidth = PG_GETARG_INT32(3);
5990 winwidth = ngbwidth * 2 + 1;
5993 if (PG_ARGISNULL(4) || PG_GETARG_INT32(4) <= 0) {
5994 elog(NOTICE,
"Neighborhood height is NULL or <= 0. Returning new "
5995 "raster with the original band");
5998 PG_FREE_IF_COPY(pgraster, 0);
6003 if (NULL == pgrtn) {
6004 elog(ERROR,
"RASTER_mapAlgebraFctNgb: Could not serialize raster");
6008 SET_VARSIZE(pgrtn, pgrtn->
size);
6009 PG_RETURN_POINTER(pgrtn);
6012 ngbheight = PG_GETARG_INT32(4);
6013 winheight = ngbheight * 2 + 1;
6016 if (PG_ARGISNULL(6)) {
6017 elog(NOTICE,
"Neighborhood NODATA behavior defaulting to 'ignore'");
6018 txtNodataMode = cstring_to_text(
"ignore");
6021 txtNodataMode = PG_GETARG_TEXT_P(6);
6024 txtCallbackParam = (text*)palloc(VARSIZE(txtNodataMode));
6025 SET_VARSIZE(txtCallbackParam, VARSIZE(txtNodataMode));
6026 memcpy((
void *)VARDATA(txtCallbackParam), (
void *)VARDATA(txtNodataMode), VARSIZE(txtNodataMode) - VARHDRSZ);
6029 #if POSTGIS_PGSQL_VERSION < 120
6030 cbdata.arg[1] = CStringGetDatum(txtCallbackParam);
6032 cbdata->args[1].value = CStringGetDatum(txtCallbackParam);
6038 if (strcmp(strFromText,
"VALUE") == 0)
6039 valuereplace =
true;
6040 else if (strcmp(strFromText,
"IGNORE") != 0 && strcmp(strFromText,
"NULL") != 0) {
6042 if (sscanf(strFromText,
"%d", &intReplace) <= 0 && sscanf(strFromText,
"%f", &fltReplace) <= 0) {
6044 elog(NOTICE,
"Neighborhood NODATA mode is not recognized. Must be one of 'value', 'ignore', "
6045 "'NULL', or a numeric value. Returning new raster with the original band");
6048 pfree(txtCallbackParam);
6052 PG_FREE_IF_COPY(pgraster, 0);
6057 if (NULL == pgrtn) {
6058 elog(ERROR,
"RASTER_mapAlgebraFctNgb: Could not serialize raster");
6062 SET_VARSIZE(pgrtn, pgrtn->
size);
6063 PG_RETURN_POINTER(pgrtn);
6066 else if (strcmp(strFromText,
"NULL") == 0) {
6075 neighborData = (Datum *)palloc(winwidth * winheight *
sizeof(Datum));
6076 neighborNulls = (
bool *)palloc(winwidth * winheight *
sizeof(
bool));
6079 neighborDims[0] = winwidth;
6080 neighborDims[1] = winheight;
6087 get_typlenbyvalalign(FLOAT8OID, &typlen, &typbyval, &typalign);
6089 for (
x = 0 + ngbwidth;
x < width - ngbwidth;
x++) {
6090 for(
y = 0 + ngbheight;
y < height - ngbheight;
y++) {
6095 pixelreplace =
false;
6099 pixelreplace =
true;
6102 for (u =
x - ngbwidth; u <=
x + ngbwidth; u++) {
6103 for (v =
y - ngbheight; v <=
y + ngbheight; v++) {
6108 neighborData[nIndex] = Float8GetDatum((
double)
r);
6109 neighborNulls[nIndex] =
false;
6110 nNodataOnly =
false;
6114 if (valuereplace && pixelreplace) {
6116 neighborData[nIndex] = Float8GetDatum((
double)rpix);
6117 neighborNulls[nIndex] =
false;
6122 neighborData[nIndex] = PointerGetDatum(NULL);
6123 neighborNulls[nIndex] =
true;
6130 neighborData[nIndex] = PointerGetDatum(NULL);
6131 neighborNulls[nIndex] =
true;
6143 if (!(nNodataOnly ||
6144 (nNullSkip && nNullItems > 0) ||
6145 (valuereplace && nNullItems > 0))) {
6147 x,
y, winwidth, winheight);
6149 neighborDatum = construct_md_array((
void *)neighborData, neighborNulls, 2, neighborDims, neighborLbs,
6150 FLOAT8OID, typlen, typbyval, typalign);
6152 #if POSTGIS_PGSQL_VERSION < 120
6154 cbdata.arg[0] = PointerGetDatum(neighborDatum);
6157 tmpnewval = FunctionCallInvoke(&cbdata);
6160 if (cbdata.isnull) {
6161 newval = newnodatavalue;
6165 cbdata->args[0].value = PointerGetDatum(neighborDatum);
6168 tmpnewval = FunctionCallInvoke(cbdata);
6173 newval = newnodatavalue;
6177 newval = DatumGetFloat8(tmpnewval);
6193 pfree(neighborNulls);
6194 pfree(neighborData);
6196 pfree(txtCallbackParam);
6199 PG_FREE_IF_COPY(pgraster, 0);
6203 POSTGIS_RT_DEBUG(3,
"RASTER_mapAlgebraFctNgb: raster modified, serializing it.");
6214 SET_VARSIZE(pgrtn, pgrtn->
size);
6215 PG_RETURN_POINTER(pgrtn);
6218 #define ARGKWCOUNT 8
6226 const uint32_t set_count = 2;
6228 int pgrastpos[2] = {-1, -1};
6231 int _isempty[2] = {0};
6232 uint32_t bandindex[2] = {0};
6235 int _hasnodata[2] = {0};
6236 double _nodataval[2] = {0};
6237 double _offset[4] = {0.};
6238 double _rastoffset[2][4] = {{0.}};
6239 int _haspixel[2] = {0};
6240 double _pixel[2] = {0};
6241 int _pos[2][2] = {{0}};
6242 uint16_t _dim[2][2] = {{0}};
6244 char *pixtypename = NULL;
6246 char *extenttypename = NULL;
6251 uint16_t dim[2] = {0};
6254 double nodataval = 0;
6255 double gt[6] = {0.};
6257 Oid calltype = InvalidOid;
6259 const uint32_t spi_count = 3;
6260 uint16_t spi_exprpos[3] = {4, 7, 8};
6261 uint32_t spi_argcount[3] = {0};
6264 SPIPlanPtr spi_plan[3] = {NULL};
6265 uint16_t spi_empty = 0;
6266 Oid *argtype = NULL;
6267 uint8_t argpos[3][8] = {{0}};
6268 char *argkw[] = {
"[rast1.x]",
"[rast1.y]",
"[rast1.val]",
"[rast1]",
"[rast2.x]",
"[rast2.y]",
"[rast2.val]",
"[rast2]"};
6272 SPITupleTable *tuptable = NULL;
6275 bool isnull =
FALSE;
6276 int hasargval[3] = {0};
6277 double argval[3] = {0.};
6278 int hasnodatanodataval = 0;
6279 double nodatanodataval = 0;
6282 Oid ufc_noid = InvalidOid;
6284 #if POSTGIS_PGSQL_VERSION < 120
6285 FunctionCallInfoData ufc_info;
6287 LOCAL_FCINFO(ufc_info, FUNC_MAX_ARGS);
6289 int ufc_nullcount = 0;
6305 for (i = 0, j = 0; i < set_count; i++) {
6306 if (!PG_ARGISNULL(j)) {
6307 pgrast[i] = (
rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(j));
6314 for (k = 0; k <= i; k++) {
6315 if (k < i &&
rast[k] != NULL)
6317 if (pgrastpos[k] != -1)
6318 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6320 elog(ERROR,
"RASTER_mapAlgebra2: Could not deserialize the %s raster", i < 1 ?
"first" :
"second");
6328 if (!PG_ARGISNULL(j)) {
6329 bandindex[i] = PG_GETARG_INT32(j);
6342 if (
rast[0] == NULL &&
rast[1] == NULL) {
6343 elog(NOTICE,
"The two rasters provided are NULL. Returning NULL");
6344 for (k = 0; k < set_count; k++) {
6345 if (pgrastpos[k] != -1)
6346 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6352 if (_isempty[0] && _isempty[1]) {
6353 elog(NOTICE,
"The two rasters provided are empty. Returning empty raster");
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 create empty raster");
6373 SET_VARSIZE(pgrtn, pgrtn->
size);
6374 PG_RETURN_POINTER(pgrtn);
6379 (
rast[0] == NULL || _isempty[0]) ||
6380 (
rast[1] == NULL || _isempty[1])
6383 if (
rast[0] == NULL || _isempty[0]) {
6396 if (_rast[i] != NULL)
6408 if (_rast[i] == NULL) {
6410 for (k = 0; k < set_count; k++) {
6411 if (pgrastpos[k] != -1)
6412 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6414 elog(ERROR,
"RASTER_mapAlgebra2: Could not create NODATA raster");
6448 for (k = 0; k < set_count; k++) {
6449 if (_rast[k] != NULL)
6451 if (pgrastpos[k] != -1)
6452 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6454 elog(ERROR,
"RASTER_mapAlgebra2: Could not test for alignment on the two rasters");
6458 elog(NOTICE,
"The two rasters provided do not have the same alignment. Returning NULL");
6459 for (k = 0; k < set_count; k++) {
6460 if (_rast[k] != NULL)
6462 if (pgrastpos[k] != -1)
6463 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6469 if (!PG_ARGISNULL(5)) {
6473 if (pixtype ==
PT_END ) {
6474 for (k = 0; k < set_count; k++) {
6475 if (_rast[k] != NULL)
6477 if (pgrastpos[k] != -1)
6478 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6480 elog(ERROR,
"RASTER_mapAlgebra2: Invalid pixel type: %s", pixtypename);
6486 if (!PG_ARGISNULL(6)) {
6499 for (k = 0; k < set_count; k++) {
6500 if (_rast[k] != NULL)
6502 if (pgrastpos[k] != -1)
6503 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6505 elog(ERROR,
"RASTER_mapAlgebra2: Could not get output raster of correct extent");
6510 _rastoffset[0][0] = _offset[0];
6511 _rastoffset[0][1] = _offset[1];
6512 _rastoffset[1][0] = _offset[2];
6513 _rastoffset[1][1] = _offset[3];
6521 switch (extenttype) {
6531 (extenttype ==
ET_FIRST && i == 0) ||
6535 elog(NOTICE,
"The %s raster is NULL. Returning NULL", (i != 1 ?
"FIRST" :
"SECOND"));
6536 for (k = 0; k < set_count; k++) {
6537 if (_rast[k] != NULL)
6539 if (pgrastpos[k] != -1)
6540 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6548 elog(NOTICE,
"The %s raster does not have the band at index %d. Returning no band raster of correct extent",
6549 (i != 1 ?
"FIRST" :
"SECOND"), bandindex[i]
6552 for (k = 0; k < set_count; k++) {
6553 if (_rast[k] != NULL)
6555 if (pgrastpos[k] != -1)
6556 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6561 if (!pgrtn) PG_RETURN_NULL();
6563 SET_VARSIZE(pgrtn, pgrtn->
size);
6564 PG_RETURN_POINTER(pgrtn);
6572 _isempty[0] || _isempty[1] ||
6575 elog(NOTICE,
"The two rasters provided have no intersection. Returning no band raster");
6578 if (dim[0] || dim[1]) {
6583 for (k = 0; k < set_count; k++) {
6584 if (_rast[k] != NULL)
6586 if (pgrastpos[k] != -1)
6587 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6589 elog(ERROR,
"RASTER_mapAlgebra2: Could not create no band raster");
6597 for (k = 0; k < set_count; k++) {
6598 if (_rast[k] != NULL)
6600 if (pgrastpos[k] != -1)
6601 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6606 if (!pgrtn) PG_RETURN_NULL();
6608 SET_VARSIZE(pgrtn, pgrtn->
size);
6609 PG_RETURN_POINTER(pgrtn);
6614 for (k = 0; k < set_count; k++) {
6615 if (_rast[k] != NULL)
6617 if (pgrastpos[k] != -1)
6618 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6620 elog(ERROR,
"RASTER_mapAlgebra2: ET_LAST and ET_CUSTOM are not implemented");
6630 elog(NOTICE,
"The two rasters provided do not have the respectively specified band indices. Returning no band raster of correct extent");
6632 for (k = 0; k < set_count; k++) {
6633 if (_rast[k] != NULL)
6635 if (pgrastpos[k] != -1)
6636 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6641 if (!pgrtn) PG_RETURN_NULL();
6643 SET_VARSIZE(pgrtn, pgrtn->
size);
6644 PG_RETURN_POINTER(pgrtn);
6648 for (i = 0; i < set_count; i++) {
6657 if (_band[i] == NULL) {
6658 for (k = 0; k < set_count; k++) {
6659 if (_rast[k] != NULL)
6661 if (pgrastpos[k] != -1)
6662 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6665 elog(ERROR,
"RASTER_mapAlgebra2: Could not get band %d of the %s raster",
6667 (i < 1 ?
"FIRST" :
"SECOND")
6679 if ((extenttype ==
ET_SECOND && !_isempty[1]) || _isempty[0])
6686 if (extenttype ==
ET_SECOND && !_isempty[1] && _hasnodata[1]) {
6687 nodataval = _nodataval[1];
6689 else if (!_isempty[0] && _hasnodata[0]) {
6690 nodataval = _nodataval[0];
6692 else if (!_isempty[1] && _hasnodata[1]) {
6693 nodataval = _nodataval[1];
6696 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));
6708 for (k = 0; k < set_count; k++) {
6709 if (_rast[k] != NULL)
6711 if (pgrastpos[k] != -1)
6712 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6715 elog(ERROR,
"RASTER_mapAlgebra2: Could not add new band to output raster");
6722 for (k = 0; k < set_count; k++) {
6723 if (_rast[k] != NULL)
6725 if (pgrastpos[k] != -1)
6726 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6729 elog(ERROR,
"RASTER_mapAlgebra2: Could not get newly added band of output raster");
6734 (
int) _rastoffset[0][0],
6735 (
int) _rastoffset[0][1],
6736 (
int) _rastoffset[1][0],
6737 (
int) _rastoffset[1][1]
6740 POSTGIS_RT_DEBUGF(4,
"metadata = (%f, %f, %d, %d, %f, %f, %f, %f, %d)",
6757 calltype = get_fn_expr_argtype(fcinfo->flinfo, 4);
6764 if (SPI_connect() != SPI_OK_CONNECT) {
6765 for (k = 0; k < set_count; k++) {
6766 if (_rast[k] != NULL)
6768 if (pgrastpos[k] != -1)
6769 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6772 elog(ERROR,
"RASTER_mapAlgebra2: Could not connect to the SPI manager");
6777 memset(hasargval, 0,
sizeof(
int) * spi_count);
6787 for (i = 0; i < spi_count; i++) {
6788 if (!PG_ARGISNULL(spi_exprpos[i])) {
6790 char place[5] =
"$1";
6805 sprintf(place,
"$%d", k);
6811 len = strlen(
"SELECT (") + strlen(expr) + strlen(
")::double precision");
6812 sql = (
char *) palloc(len + 1);
6815 for (k = 0; k < spi_count; k++) SPI_freeplan(spi_plan[k]);
6818 for (k = 0; k < set_count; k++) {
6819 if (_rast[k] != NULL)
6821 if (pgrastpos[k] != -1)
6822 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6826 elog(ERROR,
"RASTER_mapAlgebra2: Could not allocate memory for expression parameter %d", spi_exprpos[i]);
6830 memcpy(
sql,
"SELECT (", strlen(
"SELECT ("));
6831 memcpy(
sql + strlen(
"SELECT ("), expr, strlen(expr));
6832 memcpy(
sql + strlen(
"SELECT (") + strlen(expr),
")::double precision", strlen(
")::double precision"));
6838 if (spi_argcount[i]) {
6839 argtype = (Oid *) palloc(spi_argcount[i] *
sizeof(Oid));
6840 if (argtype == NULL) {
6843 for (k = 0; k < spi_count; k++) SPI_freeplan(spi_plan[k]);
6846 for (k = 0; k < set_count; k++) {
6847 if (_rast[k] != NULL)
6849 if (pgrastpos[k] != -1)
6850 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6854 elog(ERROR,
"RASTER_mapAlgebra2: Could not allocate memory for prepared plan argtypes of expression parameter %d", spi_exprpos[i]);
6860 if (argpos[i][j] < 1)
continue;
6864 (strstr(argkw[j],
"[rast1.x]") != NULL) ||
6865 (strstr(argkw[j],
"[rast1.y]") != NULL) ||
6866 (strstr(argkw[j],
"[rast2.x]") != NULL) ||
6867 (strstr(argkw[j],
"[rast2.y]") != NULL)
6869 argtype[k] = INT4OID;
6873 argtype[k] = FLOAT8OID;
6879 spi_plan[i] = SPI_prepare(
sql, spi_argcount[i], argtype);
6882 if (spi_plan[i] == NULL) {
6885 for (k = 0; k < spi_count; k++) SPI_freeplan(spi_plan[k]);
6888 for (k = 0; k < set_count; k++) {
6889 if (_rast[k] != NULL)
6891 if (pgrastpos[k] != -1)
6892 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6896 elog(ERROR,
"RASTER_mapAlgebra2: Could not create prepared plan of expression parameter %d", spi_exprpos[i]);
6902 err = SPI_execute(
sql,
TRUE, 0);
6903 if (err != SPI_OK_SELECT || SPI_tuptable == NULL || SPI_processed != 1) {
6906 for (k = 0; k < spi_count; k++) SPI_freeplan(spi_plan[k]);
6909 for (k = 0; k < set_count; k++) {
6910 if (_rast[k] != NULL)
6912 if (pgrastpos[k] != -1)
6913 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6917 elog(ERROR,
"RASTER_mapAlgebra2: Could not evaluate expression parameter %d", spi_exprpos[i]);
6922 tupdesc = SPI_tuptable->tupdesc;
6923 tuptable = SPI_tuptable;
6924 tuple = tuptable->vals[0];
6926 datum = SPI_getbinval(tuple, tupdesc, 1, &isnull);
6927 if (SPI_result == SPI_ERROR_NOATTRIBUTE) {
6930 if (SPI_tuptable) SPI_freetuptable(tuptable);
6931 for (k = 0; k < spi_count; k++) SPI_freeplan(spi_plan[k]);
6934 for (k = 0; k < set_count; k++) {
6935 if (_rast[k] != NULL)
6937 if (pgrastpos[k] != -1)
6938 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6942 elog(ERROR,
"RASTER_mapAlgebra2: Could not get result of expression parameter %d", spi_exprpos[i]);
6948 argval[i] = DatumGetFloat8(datum);
6951 if (SPI_tuptable) SPI_freetuptable(tuptable);
6961 if (!PG_ARGISNULL(9)) {
6962 hasnodatanodataval = 1;
6963 nodatanodataval = PG_GETARG_FLOAT8(9);
6966 hasnodatanodataval = 0;
6969 case REGPROCEDUREOID: {
6971 if (!PG_ARGISNULL(4)) {
6974 ufc_noid = PG_GETARG_OID(4);
6977 fmgr_info(ufc_noid, &ufl_info);
6981 if (ufl_info.fn_retset) {
6985 else if (ufl_info.fn_nargs < 3 || ufl_info.fn_nargs > 4) {
6995 for (k = 0; k < set_count; k++) {
6996 if (_rast[k] != NULL)
6998 if (pgrastpos[k] != -1)
6999 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
7004 elog(ERROR,
"RASTER_mapAlgebra2: Function provided must have three or four input parameters");
7006 elog(ERROR,
"RASTER_mapAlgebra2: Function provided must return double precision not resultset");
7010 if (func_volatile(ufc_noid) ==
'v') {
7011 elog(NOTICE,
"Function provided is VOLATILE. Unless required and for best performance, function should be IMMUTABLE or STABLE");
7015 #if POSTGIS_PGSQL_VERSION < 120
7016 InitFunctionCallInfoData(ufc_info, &ufl_info, ufl_info.fn_nargs, InvalidOid, NULL, NULL);
7017 memset(ufc_info.argnull,
FALSE,
sizeof(
bool) * ufl_info.fn_nargs);
7019 InitFunctionCallInfoData(
7020 *ufc_info, &ufl_info, ufl_info.fn_nargs, InvalidOid, NULL, NULL);
7021 ufc_info->args[0].isnull =
FALSE;
7022 ufc_info->args[1].isnull =
FALSE;
7023 ufc_info->args[2].isnull =
FALSE;
7024 if (ufl_info.fn_nargs == 4)
7025 ufc_info->args[3].isnull =
FALSE;
7028 if (ufl_info.fn_nargs != 4)
7032 #if POSTGIS_PGSQL_VERSION < 120
7033 if (!PG_ARGISNULL(7)) {
7034 ufc_info.arg[k] = PG_GETARG_DATUM(7);
7037 ufc_info.arg[k] = (Datum) NULL;
7038 ufc_info.argnull[k] =
TRUE;
7042 if (!PG_ARGISNULL(7))
7044 ufc_info->args[k].value = PG_GETARG_DATUM(7);
7048 ufc_info->args[k].value = (Datum)NULL;
7049 ufc_info->args[k].isnull =
TRUE;
7057 for (k = 0; k < set_count; k++) {
7058 if (_rast[k] != NULL)
7060 if (pgrastpos[k] != -1)
7061 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
7064 elog(ERROR,
"RASTER_mapAlgebra2: Invalid data type for expression or userfunction");
7072 (calltype == TEXTOID) && (
7073 (spi_empty != spi_count) || hasnodatanodataval
7076 (calltype == REGPROCEDUREOID) && (ufc_noid != InvalidOid)
7078 for (
x = 0;
x < dim[0];
x++) {
7079 for (
y = 0;
y < dim[1];
y++) {
7082 for (i = 0; i < set_count; i++) {
7087 _x = (int)
x - (
int)_rastoffset[i][0];
7088 _y = (int)
y - (
int)_rastoffset[i][1];
7091 _pos[i][0] = _x + 1;
7092 _pos[i][1] = _y + 1;
7095 if (_band[i] == NULL) {
7096 if (!_hasnodata[i]) {
7098 _pixel[i] = _nodataval[i];
7103 (_x >= 0 && _x < _dim[i][0]) &&
7104 (_y >= 0 && _y < _dim[i][1])
7109 if (calltype == TEXTOID) {
7110 for (k = 0; k < spi_count; k++) SPI_freeplan(spi_plan[k]);
7114 for (k = 0; k < set_count; k++) {
7115 if (_rast[k] != NULL)
7117 if (pgrastpos[k] != -1)
7118 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
7122 elog(ERROR,
"RASTER_mapAlgebra2: Could not get pixel of %s raster", (i < 1 ?
"FIRST" :
"SECOND"));
7126 if (!_hasnodata[i] || !isnodata)
7145 if (!_haspixel[0] && !_haspixel[1])
7149 else if (_haspixel[0] && !_haspixel[1])
7153 else if (!_haspixel[0] && _haspixel[1])
7162 if (hasnodatanodataval) {
7164 pixel = nodatanodataval;
7168 else if (hasargval[i]) {
7173 else if (spi_plan[i] != NULL) {
7177 memset(values, (Datum) NULL,
sizeof(Datum) *
ARGKWCOUNT);
7182 if (spi_argcount[i]) {
7186 if (idx < 1)
continue;
7189 if (strstr(argkw[j],
"[rast1.x]") != NULL) {
7190 values[idx] = _pos[0][0];
7192 else if (strstr(argkw[j],
"[rast1.y]") != NULL) {
7193 values[idx] = _pos[0][1];
7196 (strstr(argkw[j],
"[rast1.val]") != NULL) ||
7197 (strstr(argkw[j],
"[rast1]") != NULL)
7199 if (_isempty[0] || !_haspixel[0])
7202 values[idx] = Float8GetDatum(_pixel[0]);
7204 else if (strstr(argkw[j],
"[rast2.x]") != NULL) {
7205 values[idx] = _pos[1][0];
7207 else if (strstr(argkw[j],
"[rast2.y]") != NULL) {
7208 values[idx] = _pos[1][1];
7211 (strstr(argkw[j],
"[rast2.val]") != NULL) ||
7212 (strstr(argkw[j],
"[rast2]") != NULL)
7214 if (_isempty[1] || !_haspixel[1])
7217 values[idx] = Float8GetDatum(_pixel[1]);
7223 err = SPI_execute_plan(spi_plan[i], values, nulls,
TRUE, 1);
7224 if (err != SPI_OK_SELECT || SPI_tuptable == NULL || SPI_processed != 1) {
7226 for (k = 0; k < spi_count; k++) SPI_freeplan(spi_plan[k]);
7229 for (k = 0; k < set_count; k++) {
7230 if (_rast[k] != NULL)
7232 if (pgrastpos[k] != -1)
7233 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
7237 elog(ERROR,
"RASTER_mapAlgebra2: Unexpected error when running prepared statement %d", i);
7242 tupdesc = SPI_tuptable->tupdesc;
7243 tuptable = SPI_tuptable;
7244 tuple = tuptable->vals[0];
7246 datum = SPI_getbinval(tuple, tupdesc, 1, &isnull);
7247 if (SPI_result == SPI_ERROR_NOATTRIBUTE) {
7249 if (SPI_tuptable) SPI_freetuptable(tuptable);
7250 for (k = 0; k < spi_count; k++) SPI_freeplan(spi_plan[k]);
7253 for (k = 0; k < set_count; k++) {
7254 if (_rast[k] != NULL)
7256 if (pgrastpos[k] != -1)
7257 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
7261 elog(ERROR,
"RASTER_mapAlgebra2: Could not get result of prepared statement %d", i);
7267 pixel = DatumGetFloat8(datum);
7270 if (SPI_tuptable) SPI_freetuptable(tuptable);
7273 case REGPROCEDUREOID: {
7278 for (i = 0; i < set_count; i++) {
7279 #if POSTGIS_PGSQL_VERSION < 120
7280 ufc_info.arg[i] = Float8GetDatum(_pixel[i]);
7282 ufc_info->args[i].value = Float8GetDatum(_pixel[i]);
7286 #if POSTGIS_PGSQL_VERSION < 120
7287 ufc_info.argnull[i] =
FALSE;
7289 ufc_info->args[i].isnull =
FALSE;
7294 #if POSTGIS_PGSQL_VERSION < 120
7295 ufc_info.argnull[i] =
TRUE;
7297 ufc_info->args[i].isnull =
TRUE;
7305 if (ufl_info.fn_strict && ufc_nullcount)
7309 if (ufl_info.fn_nargs == 4) {
7312 for (i = 0; i < set_count; i++) {
7314 d[0] = Int32GetDatum(_pos[i][0]);
7315 d[1] = Int32GetDatum(_pos[i][1]);
7318 d[2] = Int32GetDatum(_pos[i][0]);
7319 d[3] = Int32GetDatum(_pos[i][1]);
7323 a = construct_array(d, 4, INT4OID,
sizeof(
int32),
true,
'i');
7324 #if POSTGIS_PGSQL_VERSION < 120
7325 ufc_info.arg[2] = PointerGetDatum(a);
7326 ufc_info.argnull[2] =
FALSE;
7328 ufc_info->args[2].value = PointerGetDatum(a);
7329 ufc_info->args[2].isnull =
FALSE;
7333 #if POSTGIS_PGSQL_VERSION < 120
7334 datum = FunctionCallInvoke(&ufc_info);
7337 if (!ufc_info.isnull) {
7339 pixel = DatumGetFloat8(datum);
7342 datum = FunctionCallInvoke(ufc_info);
7345 if (!ufc_info->isnull)
7348 pixel = DatumGetFloat8(datum);
7358 if (calltype == TEXTOID) {
7359 for (k = 0; k < spi_count; k++) SPI_freeplan(spi_plan[k]);
7363 for (k = 0; k < set_count; k++) {
7364 if (_rast[k] != NULL)
7366 if (pgrastpos[k] != -1)
7367 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
7371 elog(ERROR,
"RASTER_mapAlgebra2: Could not set pixel value of output raster");
7383 if (calltype == TEXTOID) {
7384 for (i = 0; i < spi_count; i++) {
7385 if (spi_plan[i] != NULL) SPI_freeplan(spi_plan[i]);
7390 for (k = 0; k < set_count; k++) {
7391 if (_rast[k] != NULL)
7393 if (pgrastpos[k] != -1)
7394 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
7399 if (!pgrtn) PG_RETURN_NULL();
7403 SET_VARSIZE(pgrtn, pgrtn->
size);
7404 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)
uint8_t * lwgeom_to_wkb(const LWGEOM *geom, uint8_t variant, size_t *size_out)
Convert LWGEOM to a char* in WKB format.
#define SRID_UNKNOWN
Unknown SRID value.
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_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.
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 * 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