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"
81 #if defined(__clang__)
82 # pragma clang diagnostic push
83 # pragma clang diagnostic ignored "-Wgnu-variable-sized-type-not-at-end"
93 char fcinfo_data[SizeForFunctionCallInfo(FUNC_MAX_ARGS)];
98 #if defined(__clang__)
99 # pragma clang diagnostic pop
131 elog(ERROR,
"rtpg_nmapalgebra_arg_init: Could not allocate memory for arguments");
166 if (arg->
raster != NULL) {
183 if( arg->
mask != NULL )
206 if (arg == NULL || array == NULL) {
207 elog(ERROR,
"rtpg_nmapalgebra_rastbandarg_process: NULL values not permitted for parameters");
211 etype = ARR_ELEMTYPE(array);
212 get_typlenbyvalalign(etype, &typlen, &typbyval, &typalign);
217 typlen, typbyval, typalign,
222 elog(ERROR,
"rtpg_nmapalgebra_rastbandarg_process: Invalid argument for rastbandarg");
240 arg->
nband == NULL ||
244 elog(ERROR,
"rtpg_nmapalgebra_rastbandarg_process: Could not allocate memory for processing rastbandarg");
253 for (i = 0; i < n; i++) {
268 tup = (HeapTupleHeader) DatumGetPointer(e[i]);
270 elog(ERROR,
"rtpg_nmapalgebra_rastbandarg_process: Invalid argument for rastbandarg at index %d", i);
276 tupv = GetAttributeByName(tup,
"rast", &isnull);
278 elog(NOTICE,
"First argument (nband) of rastbandarg at index %d is NULL. Assuming NULL raster", i);
292 for (j = 0; j < i; j++) {
304 if (arg->
raster[i] == NULL) {
305 elog(ERROR,
"rtpg_nmapalgebra_rastbandarg_process: Could not deserialize raster at index %d", i);
321 tupv = GetAttributeByName(tup,
"nband", &isnull);
324 elog(NOTICE,
"First argument (nband) of rastbandarg at index %d is NULL. Assuming nband = %d", i,
nband);
327 nband = DatumGetInt32(tupv);
330 elog(ERROR,
"rtpg_nmapalgebra_rastbandarg_process: Band number provided for rastbandarg at index %d must be greater than zero (1-based)", i);
354 arg->
nband == NULL ||
357 elog(ERROR,
"rtpg_nmapalgebra_rastbandarg_process: Could not reallocate memory for processed rastbandarg");
372 double *
value,
int *nodata
380 ArrayType *mdValues = NULL;
381 Datum *_values = NULL;
382 bool *_nodata = NULL;
384 ArrayType *mdPos = NULL;
393 int lbound[3] = {1, 1, 1};
394 Datum datum = (Datum) NULL;
408 if (_values == NULL || _nodata == NULL) {
409 elog(ERROR,
"rtpg_nmapalgebra_callback: Could not allocate memory for values array");
416 for (z = 0; z < arg->
rasters; z++) {
418 for (
y = 0;
y < arg->
rows;
y++) {
424 _nodata[i] = (bool) arg->
nodata[z][
y][
x];
426 _values[i] = Float8GetDatum(arg->
values[z][
y][
x]);
428 _values[i] = (Datum) NULL;
436 get_typlenbyvalalign(FLOAT8OID, &typlen, &typbyval, &typalign);
439 mdValues = construct_md_array(
443 typlen, typbyval, typalign
448 _pos = palloc(
sizeof(Datum) * (arg->
rasters + 1) * 2);
449 _null = palloc(
sizeof(
bool) * (arg->
rasters + 1) * 2);
450 if (_pos == NULL || _null == NULL) {
452 elog(ERROR,
"rtpg_nmapalgebra_callback: Could not allocate memory for position array");
455 memset(_null, 0,
sizeof(
bool) * (arg->
rasters + 1) * 2);
464 for (z = 0; z < arg->
rasters; z++) {
465 _pos[i] = (Datum)arg->
src_pixel[z][0] + 1;
468 _pos[i] = (Datum)arg->
src_pixel[z][1] + 1;
473 get_typlenbyvalalign(INT4OID, &typlen, &typbyval, &typalign);
481 mdPos = construct_md_array(
485 typlen, typbyval, typalign
490 callback->
ufc_info->args[0].value = PointerGetDatum(mdValues);
491 callback->
ufc_info->args[1].value = PointerGetDatum(mdPos);
494 datum = FunctionCallInvoke(callback->
ufc_info);
503 *
value = DatumGetFloat8(datum);
506 *
value = (double) DatumGetFloat4(datum);
509 *
value = (double) DatumGetInt32(datum);
512 *
value = (double) DatumGetInt16(datum);
530 ArrayType *maskArray;
561 elog(ERROR,
"RASTER_nMapAlgebra: Could not initialize argument structure");
568 elog(ERROR,
"RASTER_nMapAlgebra: Could not process rastbandarg");
572 POSTGIS_RT_DEBUGF(4,
"allnull, allempty, noband = %d, %d, %d", allnull, allempty, noband);
576 elog(NOTICE,
"All input rasters are NULL. Returning NULL");
582 if (!PG_ARGISNULL(2)) {
583 char *pixtypename = text_to_cstring(PG_GETARG_TEXT_P(2));
589 elog(ERROR,
"RASTER_nMapAlgebra: Invalid pixel type: %s", pixtypename);
595 if (!PG_ARGISNULL(3)){
596 arg->
distance[0] = PG_GETARG_INT32(3);
601 if (!PG_ARGISNULL(4)){
602 arg->
distance[1] = PG_GETARG_INT32(4);
608 elog(ERROR,
"RASTER_nMapAlgebra: Distance for X and Y axis must be greater than or equal to zero");
613 if (!PG_ARGISNULL(5)) {
621 if (PG_ARGISNULL(6)) {
622 elog(NOTICE,
"Custom extent is NULL. Returning NULL");
633 elog(ERROR,
"RASTER_nMapAlgebra: Could not deserialize custom extent");
637 elog(NOTICE,
"Custom extent is an empty raster. Returning empty raster");
642 elog(ERROR,
"RASTER_nMapAlgebra: Could not create empty raster");
648 if (!pgraster) PG_RETURN_NULL();
650 SET_VARSIZE(pgraster, pgraster->
size);
651 PG_RETURN_POINTER(pgraster);
657 if( PG_ARGISNULL(7) ){
662 maskArray = PG_GETARG_ARRAYTYPE_P(7);
663 etype = ARR_ELEMTYPE(maskArray);
664 get_typlenbyvalalign(etype,&typlen,&typbyval,&typalign);
672 elog(ERROR,
"RASTER_nMapAlgebra: Mask data type must be FLOAT8 or FLOAT4");
676 ndims = ARR_NDIM(maskArray);
679 elog(ERROR,
"RASTER_nMapAlgebra: Mask Must be a 2D array");
684 maskDims = ARR_DIMS(maskArray);
686 if (maskDims[0] % 2 == 0 || maskDims[1] % 2 == 0) {
687 elog(ERROR,
"RASTER_nMapAlgebra: Mask dimensions must be odd");
695 typlen, typbyval,typalign,
696 &maskElements,&maskNulls,&num
699 if (num < 1 || num != (maskDims[0] * maskDims[1])) {
704 elog(ERROR,
"RASTER_nMapAlgebra: Could not deconstruct new values array");
710 arg->
mask->
values = palloc(
sizeof(
double*)* maskDims[0]);
711 arg->
mask->
nodata = palloc(
sizeof(
int*)*maskDims[0]);
712 for (i = 0; i < maskDims[0]; i++) {
713 arg->
mask->
values[i] = (
double*) palloc(
sizeof(
double) * maskDims[1]);
714 arg->
mask->
nodata[i] = (
int*) palloc(
sizeof(
int) * maskDims[1]);
719 for (
y = 0;
y < maskDims[0];
y++) {
720 for (
x = 0;
x < maskDims[1];
x++) {
728 arg->
mask->
values[
y][
x] = (double) DatumGetFloat4(maskElements[i]);
732 arg->
mask->
values[
y][
x] = (double) DatumGetFloat8(maskElements[i]);
743 if (maskDims[0] == 1 && maskDims[1] == 1) {
754 if (PG_ARGISNULL(8) || !PG_GETARG_BOOL(8)) {
755 if (arg->
mask != NULL)
766 elog(NOTICE,
"All input rasters are empty. Returning empty raster");
771 elog(NOTICE,
"All input rasters do not have bands at indicated indexes. Returning empty raster");
779 elog(ERROR,
"RASTER_nMapAlgebra: Could not create empty raster");
785 if (!pgraster) PG_RETURN_NULL();
787 SET_VARSIZE(pgraster, pgraster->
size);
788 PG_RETURN_POINTER(pgraster);
792 if (!PG_ARGISNULL(1) || get_fn_expr_argtype(fcinfo->flinfo, 1) == REGPROCEDUREOID) {
811 get_func_result_type(
838 elog(ERROR,
"RASTER_nMapAlgebra: Function provided must return a double precision, float, int or smallint");
841 elog(ERROR,
"RASTER_nMapAlgebra: Function provided must return scalar (double precision, float, int, smallint)");
844 elog(ERROR,
"RASTER_nMapAlgebra: Function provided must have three input parameters");
847 elog(ERROR,
"RASTER_nMapAlgebra: Function provided must return double precision, not resultset");
854 elog(NOTICE,
"Function provided is VOLATILE. Unless required and for best performance, function should be IMMUTABLE or STABLE");
868 if (!PG_ARGISNULL(9))
874 arg->
callback.
ufc_info->args[2].value = PointerGetDatum(construct_empty_array(TEXTOID));
885 elog(ERROR,
"RASTER_nMapAlgebra: callbackfunc must be provided");
928 if (itrset == NULL) {
930 elog(ERROR,
"RASTER_nMapAlgebra: Could not allocate memory for iterator arguments");
959 elog(ERROR,
"RASTER_nMapAlgebra: Could not run raster iterator function");
973 SET_VARSIZE(pgraster, pgraster->
size);
974 PG_RETURN_POINTER(pgraster);
1018 elog(ERROR,
"rtpg_nmapalgebraexpr_arg_init: Could not allocate memory for arguments");
1024 elog(ERROR,
"rtpg_nmapalgebraexpr_arg_init: Could not allocate memory for arg->bandarg");
1037 elog(ERROR,
"rtpg_nmapalgebraexpr_arg_init: Could not allocate memory for spi_argpos");
1068 double *
value,
int *nodata
1071 SPIPlanPtr plan = NULL;
1091 else if (arg->
nodata[0][0][0] && !arg->
nodata[1][0][0]) {
1101 else if (!arg->
nodata[0][0][0] && arg->
nodata[1][0][0]) {
1128 if (arg->
nodata[0][0][0]) {
1164 SPITupleTable *tuptable = NULL;
1167 bool isnull =
FALSE;
1172 memset(values, (Datum) NULL,
sizeof(Datum) * callback->
kw.
count);
1173 memset(nulls,
FALSE,
sizeof(
char) * callback->
kw.
count);
1178 for (i = 0; i < callback->
kw.
count; i++) {
1180 if (idx < 1)
continue;
1183 if (arg->
rasters == 1 && i > 7) {
1184 elog(ERROR,
"rtpg_nmapalgebraexpr_callback: rast2 argument specified in single-raster invocation");
1191 values[idx] = Int32GetDatum(arg->
src_pixel[0][0] + 1);
1195 values[idx] = Int32GetDatum(arg->
src_pixel[0][1] + 1);
1201 if (!arg->
nodata[0][0][0])
1202 values[idx] = Float8GetDatum(arg->
values[0][0][0]);
1209 values[idx] = Int32GetDatum(arg->
src_pixel[0][0] + 1);
1213 values[idx] = Int32GetDatum(arg->
src_pixel[0][1] + 1);
1219 if (!arg->
nodata[0][0][0])
1220 values[idx] = Float8GetDatum(arg->
values[0][0][0]);
1227 values[idx] = Int32GetDatum(arg->
src_pixel[1][0] + 1);
1231 values[idx] = Int32GetDatum(arg->
src_pixel[1][1] + 1);
1237 if (!arg->
nodata[1][0][0])
1238 values[idx] = Float8GetDatum(arg->
values[1][0][0]);
1248 err = SPI_execute_plan(plan, values, nulls,
TRUE, 1);
1249 if (err != SPI_OK_SELECT || SPI_tuptable == NULL || SPI_processed != 1) {
1250 elog(ERROR,
"rtpg_nmapalgebraexpr_callback: Unexpected error when running prepared statement %d",
id);
1255 tupdesc = SPI_tuptable->tupdesc;
1256 tuptable = SPI_tuptable;
1257 tuple = tuptable->vals[0];
1259 datum = SPI_getbinval(tuple, tupdesc, 1, &isnull);
1260 if (SPI_result == SPI_ERROR_NOATTRIBUTE) {
1261 if (SPI_tuptable) SPI_freetuptable(tuptable);
1262 elog(ERROR,
"rtpg_nmapalgebraexpr_callback: Could not get result of prepared statement %d",
id);
1267 *
value = DatumGetFloat8(datum);
1287 if (SPI_tuptable) SPI_freetuptable(tuptable);
1297 MemoryContext mainMemCtx = CurrentMemoryContext;
1300 uint16_t exprpos[3] = {1, 4, 5};
1314 SPITupleTable *tuptable = NULL;
1317 bool isnull =
FALSE;
1323 const int argkwcount = 12;
1339 if (PG_ARGISNULL(0))
1345 elog(ERROR,
"RASTER_nMapAlgebraExpr: Could not initialize argument structure");
1352 elog(ERROR,
"RASTER_nMapAlgebra: Could not process rastbandarg");
1356 POSTGIS_RT_DEBUGF(4,
"allnull, allempty, noband = %d, %d, %d", allnull, allempty, noband);
1360 elog(NOTICE,
"All input rasters are NULL. Returning NULL");
1372 if (!PG_ARGISNULL(2)) {
1373 char *pixtypename = text_to_cstring(PG_GETARG_TEXT_P(2));
1379 elog(ERROR,
"RASTER_nMapAlgebraExpr: Invalid pixel type: %s", pixtypename);
1386 if (!PG_ARGISNULL(3)) {
1392 if (numraster < 2) {
1393 elog(NOTICE,
"CUSTOM extent type not supported. Defaulting to FIRST");
1397 elog(NOTICE,
"CUSTOM extent type not supported. Defaulting to INTERSECTION");
1401 else if (numraster < 2)
1407 if (!PG_ARGISNULL(6)) {
1415 elog(NOTICE,
"All input rasters are empty. Returning empty raster");
1420 elog(NOTICE,
"All input rasters do not have bands at indicated indexes. Returning empty raster");
1428 elog(ERROR,
"RASTER_nMapAlgebraExpr: Could not create empty raster");
1434 if (!pgraster) PG_RETURN_NULL();
1436 SET_VARSIZE(pgraster, pgraster->
size);
1437 PG_RETURN_POINTER(pgraster);
1441 if (SPI_connect() != SPI_OK_CONNECT) {
1443 elog(ERROR,
"RASTER_nMapAlgebraExpr: Could not connect to the SPI manager");
1459 char place[12] =
"$1";
1461 if (PG_ARGISNULL(exprpos[i]))
1464 expr = text_to_cstring(PG_GETARG_TEXT_P(exprpos[i]));
1467 for (j = 0, k = 1; j < argkwcount; j++) {
1479 sprintf(place,
"$%d", k);
1485 len = strlen(
"SELECT (") + strlen(expr) + strlen(
")::double precision");
1486 sql = (
char *) palloc(len + 1);
1490 elog(ERROR,
"RASTER_nMapAlgebraExpr: Could not allocate memory for expression parameter %d", exprpos[i]);
1494 memcpy(
sql,
"SELECT (", strlen(
"SELECT ("));
1495 memcpy(
sql + strlen(
"SELECT ("), expr, strlen(expr));
1496 memcpy(
sql + strlen(
"SELECT (") + strlen(expr),
")::double precision", strlen(
")::double precision"));
1505 if (argtype == NULL) {
1509 elog(ERROR,
"RASTER_nMapAlgebraExpr: Could not allocate memory for prepared plan argtypes of expression parameter %d", exprpos[i]);
1514 for (j = 0, k = 0; j < argkwcount; j++) {
1519 (strstr(argkw[j],
"[rast.x]") != NULL) ||
1520 (strstr(argkw[j],
"[rast.y]") != NULL) ||
1521 (strstr(argkw[j],
"[rast1.x]") != NULL) ||
1522 (strstr(argkw[j],
"[rast1.y]") != NULL) ||
1523 (strstr(argkw[j],
"[rast2.x]") != NULL) ||
1524 (strstr(argkw[j],
"[rast2.y]") != NULL)
1526 argtype[k] = INT4OID;
1529 argtype[k] = FLOAT8OID;
1541 elog(ERROR,
"RASTER_nMapAlgebraExpr: Could not create prepared plan of expression parameter %d", exprpos[i]);
1547 POSTGIS_RT_DEBUGF(3,
"expression parameter %d has no args, simply executing", exprpos[i]);
1548 err = SPI_execute(
sql,
TRUE, 0);
1551 if (err != SPI_OK_SELECT || SPI_tuptable == NULL || SPI_processed != 1) {
1554 elog(ERROR,
"RASTER_nMapAlgebraExpr: Could not evaluate expression parameter %d", exprpos[i]);
1559 tupdesc = SPI_tuptable->tupdesc;
1560 tuptable = SPI_tuptable;
1561 tuple = tuptable->vals[0];
1563 datum = SPI_getbinval(tuple, tupdesc, 1, &isnull);
1564 if (SPI_result == SPI_ERROR_NOATTRIBUTE) {
1565 if (SPI_tuptable) SPI_freetuptable(tuptable);
1568 elog(ERROR,
"RASTER_nMapAlgebraExpr: Could not get result of expression parameter %d", exprpos[i]);
1577 if (SPI_tuptable) SPI_freetuptable(tuptable);
1597 for (i = 0; i < numraster; i++) {
1621 if (itrset == NULL) {
1624 elog(ERROR,
"RASTER_nMapAlgebra: Could not allocate memory for iterator arguments");
1629 for (i = 0; i < numraster; i++) {
1653 elog(ERROR,
"RASTER_nMapAlgebraExpr: Could not run raster iterator function");
1656 else if (
raster == NULL) {
1662 MemoryContextSwitchTo(mainMemCtx);
1673 SET_VARSIZE(pgraster, pgraster->
size);
1674 PG_RETURN_POINTER(pgraster);
1694 assert(cutype && strlen(cutype) > 0);
1696 if (strcmp(cutype,
"LAST") == 0)
1698 else if (strcmp(cutype,
"FIRST") == 0)
1700 else if (strcmp(cutype,
"MIN") == 0)
1702 else if (strcmp(cutype,
"MAX") == 0)
1704 else if (strcmp(cutype,
"COUNT") == 0)
1706 else if (strcmp(cutype,
"SUM") == 0)
1708 else if (strcmp(cutype,
"MEAN") == 0)
1710 else if (strcmp(cutype,
"RANGE") == 0)
1737 for (i = 0; i < arg->
numband; i++) {
1761 double *
value,
int *nodata
1773 elog(ERROR,
"rtpg_union_callback: Invalid arguments passed to callback");
1789 else if (!arg->
nodata[0][0][0] && arg->
nodata[1][0][0]) {
1795 else if (arg->
nodata[0][0][0] && !arg->
nodata[1][0][0]) {
1823 else if (!arg->
nodata[0][0][0] && arg->
nodata[1][0][0])
1826 else if (arg->
nodata[0][0][0] && !arg->
nodata[1][0][0])
1852 double *
value,
int *nodata
1862 elog(ERROR,
"rtpg_union_mean_callback: Invalid arguments passed to callback");
1885 double *
value,
int *nodata
1895 elog(ERROR,
"rtpg_union_range_callback: Invalid arguments passed to callback");
1928 HeapTupleHeader tup;
1934 char *utypename = NULL;
1937 etype = ARR_ELEMTYPE(array);
1938 get_typlenbyvalalign(etype, &typlen, &typbyval, &typalign);
1943 typlen, typbyval, typalign,
1948 elog(ERROR,
"rtpg_union_unionarg_process: Invalid argument for unionarg");
1956 elog(ERROR,
"rtpg_union_unionarg_process: Could not allocate memory for band information");
1961 for (i = 0; i < n; i++) {
1970 tup = (HeapTupleHeader) DatumGetPointer(e[i]);
1972 elog(ERROR,
"rtpg_union_unionarg_process: Invalid argument for unionarg");
1977 tupv = GetAttributeByName(tup,
"nband", &isnull);
1980 elog(NOTICE,
"First argument (nband) of unionarg is NULL. Assuming nband = %d",
nband);
1983 nband = DatumGetInt32(tupv);
1986 elog(ERROR,
"rtpg_union_unionarg_process: Band number must be greater than zero (1-based)");
1991 tupv = GetAttributeByName(tup,
"uniontype", &isnull);
1993 elog(NOTICE,
"Second argument (uniontype) of unionarg is NULL. Assuming uniontype = LAST");
1997 utypename = text_to_cstring((text *) DatumGetPointer(tupv));
2018 elog(ERROR,
"rtpg_union_unionarg_process: Could not reallocate memory for band information");
2035 if (numbands <= arg->numband)
2045 elog(ERROR,
"rtpg_union_noarg: Could not reallocate memory for band information");
2051 for (; i < arg->
numband; i++) {
2059 elog(ERROR,
"rtpg_union_noarg: Could not allocate memory for working rasters");
2068 elog(ERROR,
"rtpg_union_noarg: Could not create working raster");
2081 MemoryContext aggcontext;
2082 MemoryContext oldcontext;
2092 int isempty[2] = {0};
2093 int hasband[2] = {0};
2095 double _offset[4] = {0.};
2103 char *utypename = NULL;
2107 double nodataval = 0;
2113 uint16_t _dim[2] = {0};
2120 if (!AggCheckCallContext(fcinfo, &aggcontext)) {
2121 elog(ERROR,
"RASTER_union_transfn: Cannot be called in a non-aggregate context");
2126 oldcontext = MemoryContextSwitchTo(aggcontext);
2128 if (PG_ARGISNULL(0)) {
2133 MemoryContextSwitchTo(oldcontext);
2134 elog(ERROR,
"RASTER_union_transfn: Could not allocate memory for state variable");
2150 if (!PG_ARGISNULL(1)) {
2152 pgraster = (
rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
2159 PG_FREE_IF_COPY(pgraster, 1);
2161 MemoryContextSwitchTo(oldcontext);
2162 elog(ERROR,
"RASTER_union_transfn: Could not deserialize raster");
2175 if (!PG_ARGISNULL(2)) {
2176 Oid calltype = get_fn_expr_argtype(fcinfo->flinfo, 2);
2187 utypename = text_to_cstring(PG_GETARG_TEXT_P(2));
2211 if (numband > idx) {
2229 PG_FREE_IF_COPY(pgraster, 1);
2232 MemoryContextSwitchTo(oldcontext);
2233 elog(ERROR,
"RASTER_union_transfn: Could not allocate memory for band information");
2238 for (i = idx; i < iwr->
numband; i++) {
2262 nband = PG_GETARG_INT32(2);
2268 PG_FREE_IF_COPY(pgraster, 1);
2271 MemoryContextSwitchTo(oldcontext);
2272 elog(ERROR,
"RASTER_union_transfn: Band number must be greater than zero (1-based)");
2283 PG_FREE_IF_COPY(pgraster, 1);
2286 MemoryContextSwitchTo(oldcontext);
2287 elog(ERROR,
"RASTER_union_transfn: Could not allocate memory for band information");
2308 PG_FREE_IF_COPY(pgraster, 1);
2311 MemoryContextSwitchTo(oldcontext);
2312 elog(ERROR,
"RASTER_union_transfn: Could not process unionarg");
2321 if (nargs > 3 && !PG_ARGISNULL(3)) {
2322 utypename = text_to_cstring(PG_GETARG_TEXT_P(3));
2336 for (i = 0; i < iwr->
numband; i++) {
2351 PG_FREE_IF_COPY(pgraster, 1);
2354 MemoryContextSwitchTo(oldcontext);
2355 elog(ERROR,
"RASTER_union_transfn: Could not allocate memory for working raster(s)");
2370 PG_FREE_IF_COPY(pgraster, 1);
2373 MemoryContextSwitchTo(oldcontext);
2374 elog(ERROR,
"RASTER_union_transfn: Could not create working raster");
2391 PG_FREE_IF_COPY(pgraster, 1);
2394 MemoryContextSwitchTo(oldcontext);
2395 elog(ERROR,
"RASTER_union_transfn: Could not check and balance number of bands");
2402 if (itrset == NULL) {
2407 PG_FREE_IF_COPY(pgraster, 1);
2410 MemoryContextSwitchTo(oldcontext);
2411 elog(ERROR,
"RASTER_union_transfn: Could not allocate memory for iterator arguments");
2416 for (i = 0; i < iwr->
numband; i++) {
2436 if (!isempty[0] && hasband[0])
2438 else if (!isempty[1] && hasband[1])
2445 if (_band != NULL) {
2483 itrset[0].
nband = 0;
2499 if (!isempty[0] && !isempty[1] && hasband[0] && hasband[1]) {
2513 PG_FREE_IF_COPY(pgraster, 1);
2516 MemoryContextSwitchTo(oldcontext);
2517 elog(ERROR,
"RASTER_union_transfn: Could not create internal raster");
2521 _offset[0], _offset[1], _offset[2], _offset[3]);
2528 double igt[6] = {0};
2544 hasnodata, nodataval,
2553 PG_FREE_IF_COPY(pgraster, 1);
2556 MemoryContextSwitchTo(oldcontext);
2557 elog(ERROR,
"RASTER_union_transfn: Could not add new band to internal raster");
2565 for (
y = 0;
y < _dim[1];
y++) {
2566 POSTGIS_RT_DEBUGF(4,
"Getting pixel line of working raster at (x, y, length) = (0, %d, %d)",
y, _dim[0]);
2580 PG_FREE_IF_COPY(pgraster, 1);
2583 MemoryContextSwitchTo(oldcontext);
2584 elog(ERROR,
"RASTER_union_transfn: Could not get pixel line from band of working raster");
2588 POSTGIS_RT_DEBUGF(4,
"Setting pixel line at (x, y, length) = (%d, %d, %d)", (
int) _offset[0], (
int) _offset[1] +
y, nvals);
2591 (
int) _offset[0], (
int) _offset[1] +
y,
2601 PG_FREE_IF_COPY(pgraster, 1);
2604 MemoryContextSwitchTo(oldcontext);
2605 elog(ERROR,
"RASTER_union_transfn: Could not set pixel line to band of internal raster");
2621 hasnodata, nodataval,
2638 PG_FREE_IF_COPY(pgraster, 1);
2641 MemoryContextSwitchTo(oldcontext);
2642 elog(ERROR,
"RASTER_union_transfn: Could not run raster iterator function");
2650 for (
y = 0;
y < _dim[1];
y++) {
2651 POSTGIS_RT_DEBUGF(4,
"Getting pixel line of iterator raster at (x, y, length) = (0, %d, %d)",
y, _dim[0]);
2667 PG_FREE_IF_COPY(pgraster, 1);
2670 MemoryContextSwitchTo(oldcontext);
2671 elog(ERROR,
"RASTER_union_transfn: Could not get pixel line from band of working raster");
2675 POSTGIS_RT_DEBUGF(4,
"Setting pixel line at (x, y, length) = (%d, %d, %d)", (
int) _offset[2], (
int) _offset[3] +
y, nvals);
2678 (
int) _offset[2], (
int) _offset[3] +
y,
2690 PG_FREE_IF_COPY(pgraster, 1);
2693 MemoryContextSwitchTo(oldcontext);
2694 elog(ERROR,
"RASTER_union_transfn: Could not set pixel line to band of internal raster");
2714 hasnodata, nodataval,
2728 PG_FREE_IF_COPY(pgraster, 1);
2731 MemoryContextSwitchTo(oldcontext);
2732 elog(ERROR,
"RASTER_union_transfn: Could not run raster iterator function");
2751 PG_FREE_IF_COPY(pgraster, 1);
2755 MemoryContextSwitchTo(oldcontext);
2759 PG_RETURN_POINTER(iwr);
2779 double nodataval = 0;
2784 if (!AggCheckCallContext(fcinfo, NULL)) {
2785 elog(ERROR,
"RASTER_union_finalfn: Cannot be called in a non-aggregate context");
2790 if (PG_ARGISNULL(0))
2797 if (itrset == NULL) {
2799 elog(ERROR,
"RASTER_union_finalfn: Could not allocate memory for iterator arguments");
2803 for (i = 0; i < iwr->
numband; i++) {
2818 itrset[0].
nband = 0;
2820 itrset[1].
nband = 0;
2828 hasnodata, nodataval,
2841 hasnodata, nodataval,
2855 elog(ERROR,
"RASTER_union_finalfn: Could not run raster iterator function");
2861 if (_raster == NULL)
2867 uint32_t bandNums[1] = {0};
2869 status = (_rtn == NULL) ? -1 : 0;
2894 elog(ERROR,
"RASTER_union_finalfn: Could not add band to final raster");
2907 if (!_rtn) PG_RETURN_NULL();
2917 SET_VARSIZE(pgraster, pgraster->
size);
2918 PG_RETURN_POINTER(pgraster);
2947 elog(ERROR,
"rtpg_clip_arg_init: Could not allocate memory for function arguments");
2961 if (arg->
band != NULL)
2966 if (arg->
mask != NULL)
3006 double offset[4] = {0.};
3016 int mask_height = 0;
3023 char **options = NULL;
3025 int options_len = 0;
3030 if (PG_ARGISNULL(0) || PG_ARGISNULL(2))
3036 elog(ERROR,
"RASTER_clip: Could not initialize argument structure");
3041 pgraster = (
rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
3045 if (arg->
raster == NULL) {
3047 PG_FREE_IF_COPY(pgraster, 0);
3048 elog(ERROR,
"RASTER_clip: Could not deserialize raster");
3054 elog(NOTICE,
"Input raster is empty or has no bands. Returning empty raster");
3057 PG_FREE_IF_COPY(pgraster, 0);
3061 elog(ERROR,
"RASTER_clip: Could not create empty raster");
3070 SET_VARSIZE(pgrtn, pgrtn->
size);
3071 PG_RETURN_POINTER(pgrtn);
3079 gser = PG_GETARG_GSERIALIZED_P(2);
3091 elog(NOTICE,
"Geometry provided does not have the same SRID as the raster. Returning NULL");
3094 PG_FREE_IF_COPY(pgraster, 0);
3096 PG_FREE_IF_COPY(gser, 2);
3102 if (!PG_ARGISNULL(4) && !PG_GETARG_BOOL(4))
3109 PG_FREE_IF_COPY(pgraster, 0);
3111 PG_FREE_IF_COPY(gser, 2);
3113 elog(ERROR,
"RASTER_clip: Could not get convex hull of raster");
3120 PG_FREE_IF_COPY(gser, 2);
3125 elog(NOTICE,
"The input raster and input geometry do not intersect. Returning empty raster");
3128 PG_FREE_IF_COPY(pgraster, 0);
3133 elog(ERROR,
"RASTER_clip: Could not create empty raster");
3142 SET_VARSIZE(pgrtn, pgrtn->
size);
3143 PG_RETURN_POINTER(pgrtn);
3147 if (!PG_ARGISNULL(1)) {
3148 array = PG_GETARG_ARRAYTYPE_P(1);
3149 etype = ARR_ELEMTYPE(array);
3150 get_typlenbyvalalign(etype, &typlen, &typbyval, &typalign);
3158 PG_FREE_IF_COPY(pgraster, 0);
3160 elog(ERROR,
"RASTER_clip: Invalid data type for band indexes");
3167 typlen, typbyval, typalign,
3172 if (arg->
band == NULL) {
3174 PG_FREE_IF_COPY(pgraster, 0);
3176 elog(ERROR,
"RASTER_clip: Could not allocate memory for band arguments");
3180 for (i = 0, j = 0; i < arg->
numbands; i++) {
3181 if (nulls[i])
continue;
3185 arg->
band[j].
nband = DatumGetInt16(e[i]) - 1;
3188 arg->
band[j].
nband = DatumGetInt32(e[i]) - 1;
3195 if (j < arg->numbands) {
3197 if (arg->
band == NULL) {
3199 PG_FREE_IF_COPY(pgraster, 0);
3201 elog(ERROR,
"RASTER_clip: Could not reallocate memory for band arguments");
3209 for (i = 0; i < arg->
numbands; i++) {
3211 elog(NOTICE,
"Band at index %d not found in raster", arg->
band[i].
nband + 1);
3213 PG_FREE_IF_COPY(pgraster, 0);
3228 if (arg->
band == NULL) {
3231 PG_FREE_IF_COPY(pgraster, 0);
3234 elog(ERROR,
"RASTER_clip: Could not allocate memory for band arguments");
3238 for (i = 0; i < arg->
numbands; i++) {
3247 if (!PG_ARGISNULL(3)) {
3248 array = PG_GETARG_ARRAYTYPE_P(3);
3249 etype = ARR_ELEMTYPE(array);
3250 get_typlenbyvalalign(etype, &typlen, &typbyval, &typalign);
3258 PG_FREE_IF_COPY(pgraster, 0);
3260 elog(ERROR,
"RASTER_clip: Invalid data type for NODATA values");
3267 typlen, typbyval, typalign,
3272 for (i = 0, j = 0; i < arg->
numbands; i++, j++) {
3298 if (!PG_ARGISNULL(5) && PG_GETARG_BOOL(5) ==
TRUE){
3299 options_len = options_len + 1;
3300 options = (
char **) palloc(
sizeof(
char *) * options_len);
3301 options[options_len - 1] = palloc(
sizeof(
char*) * (strlen(
"ALL_TOUCHED=TRUE") + 1));
3302 strcpy(options[options_len - 1],
"ALL_TOUCHED=TRUE");
3307 options = (
char **) repalloc(options,
sizeof(
char *) * options_len);
3308 options[options_len - 1] = NULL;
3335 if (options_len) pfree(options);
3337 if (arg->
mask == NULL) {
3339 PG_FREE_IF_COPY(pgraster, 0);
3340 elog(ERROR,
"RASTER_clip: Could not rasterize intersection geometry");
3352 PG_FREE_IF_COPY(pgraster, 0);
3353 elog(ERROR,
"RASTER_clip: Could not compute extent of rasters");
3362 for (i = 0; i < arg->
numbands; i++) {
3383 PG_FREE_IF_COPY(pgraster, 0);
3384 elog(ERROR,
"RASTER_clip: Could not add new band in output raster");
3398 for (
y = 0;
y < height;
y++) {
3399 for (
x = 0;
x < width;
x++) {
3400 mask_x =
x - (int)offset[2];
3401 mask_y =
y - (int)offset[3];
3405 mask_x < mask_width &&
3407 mask_y < mask_height
3414 PG_FREE_IF_COPY(pgraster, 0);
3415 elog(ERROR,
"RASTER_clip: Could not get pixel value");
3423 input_x =
x - (int)offset[0];
3424 input_y =
y - (int)offset[1];
3428 PG_FREE_IF_COPY(pgraster, 0);
3429 elog(ERROR,
"RASTER_clip: Could not get pixel value");
3439 PG_FREE_IF_COPY(pgraster, 0);
3440 elog(ERROR,
"RASTER_clip: Could not set pixel value");
3448 PG_FREE_IF_COPY(pgraster, 0);
3458 SET_VARSIZE(pgrtn, pgrtn->
size);
3459 PG_RETURN_POINTER(pgrtn);
3472 uint32_t numBands = 0;
3492 HeapTupleHeader tup;
3497 text *exprtext = NULL;
3502 char *pixeltype = NULL;
3503 text *pixeltypetext = NULL;
3505 double nodataval = 0;
3506 bool hasnodata =
FALSE;
3508 char **comma_set = NULL;
3509 uint32_t comma_n = 0;
3510 char **colon_set = NULL;
3511 uint32_t colon_n = 0;
3512 char **dash_set = NULL;
3513 uint32_t dash_n = 0;
3518 if (PG_ARGISNULL(0))
3520 pgraster = (
rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
3525 PG_FREE_IF_COPY(pgraster, 0);
3526 elog(ERROR,
"RASTER_reclass: Could not deserialize raster");
3530 POSTGIS_RT_DEBUGF(3,
"RASTER_reclass: %d possible bands to be reclassified", numBands);
3534 array = PG_GETARG_ARRAYTYPE_P(1);
3535 etype = ARR_ELEMTYPE(array);
3536 get_typlenbyvalalign(etype, &typlen, &typbyval, &typalign);
3538 deconstruct_array(array, etype, typlen, typbyval, typalign, &e,
3542 elog(NOTICE,
"Invalid argument for reclassargset. Returning original raster");
3546 PG_FREE_IF_COPY(pgraster, 0);
3550 SET_VARSIZE(pgrtn, pgrtn->
size);
3551 PG_RETURN_POINTER(pgrtn);
3559 for (i = 0; i < n; i++) {
3560 if (nulls[i])
continue;
3563 tup = (HeapTupleHeader) DatumGetPointer(e[i]);
3565 elog(NOTICE,
"Invalid argument for reclassargset. Returning original raster");
3569 PG_FREE_IF_COPY(pgraster, 0);
3573 SET_VARSIZE(pgrtn, pgrtn->
size);
3574 PG_RETURN_POINTER(pgrtn);
3578 tupv = GetAttributeByName(tup,
"nband", &isnull);
3580 elog(NOTICE,
"Invalid argument for reclassargset. Missing value of nband for reclassarg of index %d . Returning original raster", i);
3584 PG_FREE_IF_COPY(pgraster, 0);
3588 SET_VARSIZE(pgrtn, pgrtn->
size);
3589 PG_RETURN_POINTER(pgrtn);
3591 nband = DatumGetInt32(tupv);
3595 if (nband < 1 || nband > numBands) {
3596 elog(NOTICE,
"Invalid argument for reclassargset. Invalid band index (must use 1-based) for reclassarg of index %d . Returning original raster", i);
3600 PG_FREE_IF_COPY(pgraster, 0);
3604 SET_VARSIZE(pgrtn, pgrtn->
size);
3605 PG_RETURN_POINTER(pgrtn);
3609 tupv = GetAttributeByName(tup,
"reclassexpr", &isnull);
3611 elog(NOTICE,
"Invalid argument for reclassargset. Missing value of reclassexpr for reclassarg of index %d . Returning original raster", i);
3615 PG_FREE_IF_COPY(pgraster, 0);
3619 SET_VARSIZE(pgrtn, pgrtn->
size);
3620 PG_RETURN_POINTER(pgrtn);
3622 exprtext = (text *) DatumGetPointer(tupv);
3623 if (NULL == exprtext) {
3624 elog(NOTICE,
"Invalid argument for reclassargset. Missing value of reclassexpr for reclassarg of index %d . Returning original raster", i);
3628 PG_FREE_IF_COPY(pgraster, 0);
3632 SET_VARSIZE(pgrtn, pgrtn->
size);
3633 PG_RETURN_POINTER(pgrtn);
3635 expr = text_to_cstring(exprtext);
3644 elog(NOTICE,
"Invalid argument for reclassargset. Invalid expression of reclassexpr for reclassarg of index %d . Returning original raster", i);
3648 PG_FREE_IF_COPY(pgraster, 0);
3652 SET_VARSIZE(pgrtn, pgrtn->
size);
3653 PG_RETURN_POINTER(pgrtn);
3660 for (a = 0, j = 0; a < comma_n; a++) {
3666 elog(NOTICE,
"Invalid argument for reclassargset. Invalid expression of reclassexpr for reclassarg of index %d . Returning original raster", i);
3667 for (k = 0; k < j; k++) pfree(exprset[k]);
3672 PG_FREE_IF_COPY(pgraster, 0);
3676 SET_VARSIZE(pgrtn, pgrtn->
size);
3677 PG_RETURN_POINTER(pgrtn);
3683 for (b = 0; b < colon_n; b++) {
3688 if (dash_n < 1 || dash_n > 3) {
3689 elog(NOTICE,
"Invalid argument for reclassargset. Invalid expression of reclassexpr for reclassarg of index %d . Returning original raster", i);
3690 for (k = 0; k < j; k++) pfree(exprset[k]);
3695 PG_FREE_IF_COPY(pgraster, 0);
3699 SET_VARSIZE(pgrtn, pgrtn->
size);
3700 PG_RETURN_POINTER(pgrtn);
3703 for (c = 0; c < dash_n; c++) {
3707 strlen(dash_set[c]) == 1 && (
3708 strchr(dash_set[c],
'(') != NULL ||
3709 strchr(dash_set[c],
'[') != NULL ||
3710 strchr(dash_set[c],
')') != NULL ||
3711 strchr(dash_set[c],
']') != NULL
3715 junk = palloc(
sizeof(
char) * (strlen(dash_set[c + 1]) + 2));
3717 for (k = 0; k <= j; k++) pfree(exprset[k]);
3720 PG_FREE_IF_COPY(pgraster, 0);
3722 elog(ERROR,
"RASTER_reclass: Could not allocate memory");
3726 sprintf(junk,
"%s%s", dash_set[c], dash_set[c + 1]);
3728 dash_set[c] = repalloc(dash_set[c],
sizeof(
char) * (strlen(junk) + 1));
3729 strcpy(dash_set[c], junk);
3733 for (dash_it = 1; dash_it < dash_n; dash_it++) {
3734 dash_set[dash_it - 1] = repalloc(dash_set[dash_it - 1], (strlen(dash_set[dash_it]) + 1) *
sizeof(
char));
3735 strcpy(dash_set[dash_it - 1], dash_set[dash_it]);
3739 pfree(dash_set[dash_n]);
3740 dash_set = repalloc(dash_set,
sizeof(
char *) * dash_n);
3744 if (c < 1 && dash_n > 2) {
3745 elog(NOTICE,
"Invalid argument for reclassargset. Invalid expression of reclassexpr for reclassarg of index %d . Returning original raster", i);
3746 for (k = 0; k < j; k++) pfree(exprset[k]);
3751 PG_FREE_IF_COPY(pgraster, 0);
3755 SET_VARSIZE(pgrtn, pgrtn->
size);
3756 PG_RETURN_POINTER(pgrtn);
3767 strchr(dash_set[c],
')') != NULL ||
3768 strchr(dash_set[c],
']') != NULL
3773 else if (strchr(dash_set[c],
'(') != NULL){
3783 strrchr(dash_set[c],
'(') != NULL ||
3784 strrchr(dash_set[c],
'[') != NULL
3789 else if (strrchr(dash_set[c],
']') != NULL) {
3797 POSTGIS_RT_DEBUGF(4,
"RASTER_reclass: exc_val %d inc_val %d", exc_val, inc_val);
3805 val = strtod(dash_set[c], &junk);
3806 if (errno != 0 || dash_set[c] == junk) {
3807 elog(NOTICE,
"Invalid argument for reclassargset. Invalid expression of reclassexpr for reclassarg of index %d . Returning original raster", i);
3808 for (k = 0; k < j; k++) pfree(exprset[k]);
3813 PG_FREE_IF_COPY(pgraster, 0);
3817 SET_VARSIZE(pgrtn, pgrtn->
size);
3818 PG_RETURN_POINTER(pgrtn);
3824 junk = strstr(colon_set[b], dash_set[c]);
3829 "(colon_set[%d], dash_set[%d], junk) = (%s, %s, %s)",
3830 b, c, colon_set[b], dash_set[c], junk
3833 if (junk != colon_set[b]) {
3835 if (*(junk - 1) ==
'-') {
3838 ((junk - 1) == colon_set[b]) ||
3839 (*(junk - 2) ==
'-') ||
3840 (*(junk - 2) ==
'[') ||
3841 (*(junk - 2) ==
'(')
3861 exprset[j]->
src.
min = val;
3867 exprset[j]->
src.
max = val;
3877 exprset[j]->
dst.
min = val;
3880 exprset[j]->
dst.
max = val;
3888 , exprset[j]->src.min
3898 tupv = GetAttributeByName(tup,
"pixeltype", &isnull);
3900 elog(NOTICE,
"Invalid argument for reclassargset. Missing value of pixeltype for reclassarg of index %d . Returning original raster", i);
3904 PG_FREE_IF_COPY(pgraster, 0);
3908 SET_VARSIZE(pgrtn, pgrtn->
size);
3909 PG_RETURN_POINTER(pgrtn);
3911 pixeltypetext = (text *) DatumGetPointer(tupv);
3912 if (NULL == pixeltypetext) {
3913 elog(NOTICE,
"Invalid argument for reclassargset. Missing value of pixeltype for reclassarg of index %d . Returning original raster", i);
3917 PG_FREE_IF_COPY(pgraster, 0);
3921 SET_VARSIZE(pgrtn, pgrtn->
size);
3922 PG_RETURN_POINTER(pgrtn);
3924 pixeltype = text_to_cstring(pixeltypetext);
3929 tupv = GetAttributeByName(tup,
"nodataval", &isnull);
3935 nodataval = DatumGetFloat8(tupv);
3944 elog(NOTICE,
"Could not find raster band of index %d. Returning original raster",
nband);
3945 for (k = 0; k < j; k++) pfree(exprset[k]);
3950 PG_FREE_IF_COPY(pgraster, 0);
3954 SET_VARSIZE(pgrtn, pgrtn->
size);
3955 PG_RETURN_POINTER(pgrtn);
3959 for (k = 0; k < j; k++) pfree(exprset[k]);
3963 PG_FREE_IF_COPY(pgraster, 0);
3965 elog(ERROR,
"RASTER_reclass: Could not reclassify raster band of index %d",
nband);
3971 for (k = 0; k < j; k++) pfree(exprset[k]);
3976 PG_FREE_IF_COPY(pgraster, 0);
3978 elog(ERROR,
"RASTER_reclass: Could not replace raster band of index %d with reclassified band",
nband);
3986 for (k = 0; k < j; k++) pfree(exprset[k]);
3992 PG_FREE_IF_COPY(pgraster, 0);
3998 SET_VARSIZE(pgrtn, pgrtn->
size);
3999 PG_RETURN_POINTER(pgrtn);
4019 uint32_t bandnum = 0, numbands = 0;
4020 ArrayType *arraysrc, *arraydst;
4021 bool hasnodata =
false;
4022 double nodataval = 0.0;
4023 text *pixeltype = NULL;
4024 uint32_t szsrc, szdst;
4025 ArrayIterator itersrc, iterdst;
4027 bool nullsrc, nulldst;
4032 if (PG_ARGISNULL(0) || PG_ARGISNULL(1) || PG_ARGISNULL(2) || PG_ARGISNULL(3) || PG_ARGISNULL(4))
4036 pgraster = (
rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
4037 arraysrc = PG_GETARG_ARRAYTYPE_P(1);
4038 arraydst = PG_GETARG_ARRAYTYPE_P(2);
4039 bandnum = PG_GETARG_INT32(3);
4040 pixeltype = PG_GETARG_TEXT_P(4);
4041 if (!PG_ARGISNULL(5)) {
4043 nodataval = PG_GETARG_FLOAT8(5);
4047 szsrc = ArrayGetNItems(ARR_NDIM(arraysrc), ARR_DIMS(arraysrc));
4048 szdst = ArrayGetNItems(ARR_NDIM(arraydst), ARR_DIMS(arraydst));
4050 elog(ERROR,
"array lengths must be the same");
4055 PG_FREE_IF_COPY(pgraster, 0);
4056 elog(ERROR,
"%s: Could not deserialize raster", __func__);
4062 elog(ERROR,
"Raster has no bands");
4063 if (bandnum < 1 || bandnum > numbands)
4064 elog(ERROR,
"Invalid band index %d, input raster has %d bands. Band indexes are one-based.",
4070 elog(ERROR,
"Could not find raster band of index %d", bandnum);
4077 elog(ERROR,
"Unknown output pixel type '%s'", text_to_cstring(pixeltype));
4081 elog(ERROR,
"Unsupported pixtype");
4085 reclassmap->
count = 0;
4086 reclassmap->
srctype = pixtypsrc;
4087 reclassmap->
dsttype = pixtypdst;
4094 itersrc = array_create_iterator(arraysrc, 0, NULL);
4095 iterdst = array_create_iterator(arraydst, 0, NULL);
4096 while(array_iterate(itersrc, &dsrc, &nullsrc) &&
4097 array_iterate(iterdst, &ddst, &nulldst))
4099 double valsrc = nullsrc ? nodataval : (double)DatumGetFloat8(dsrc);
4100 double valdst = nulldst ? nodataval : (double)DatumGetFloat8(ddst);
4101 Assert(szdst > reclassmap->
count);
4104 reclassmap->
count++;
4106 array_free_iterator(itersrc);
4107 array_free_iterator(iterdst);
4112 pfree(reclassmap->
pairs);
4115 elog(ERROR,
"Band reclassification failed");
4121 PG_FREE_IF_COPY(pgraster, 0);
4122 elog(ERROR,
"Could not replace raster band of index %d with reclassified band", bandnum);
4127 PG_FREE_IF_COPY(pgraster, 0);
4131 SET_VARSIZE(pgrtn, pgrtn->
size);
4132 PG_RETURN_POINTER(pgrtn);
4165 elog(ERROR,
"rtpg_colormap_arg: Could not allocate memory for function arguments");
4176 elog(ERROR,
"rtpg_colormap_arg: Could not allocate memory for function arguments");
4209 for (i = 0; i < arg->
nentry; i++) {
4210 if (arg->
entry[i] != NULL)
4211 pfree(arg->
entry[i]);
4217 for (i = 0; i < arg->
nelement; i++)
4237 if (PG_ARGISNULL(0))
4243 elog(ERROR,
"RASTER_colorMap: Could not initialize argument structure");
4248 pgraster = (
rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
4254 PG_FREE_IF_COPY(pgraster, 0);
4255 elog(ERROR,
"RASTER_colorMap: Could not deserialize raster");
4260 if (!PG_ARGISNULL(1))
4261 arg->
nband = PG_GETARG_INT32(1);
4266 elog(NOTICE,
"Raster does not have band at index %d. Returning empty raster", arg->
nband);
4271 PG_FREE_IF_COPY(pgraster, 0);
4272 elog(ERROR,
"RASTER_colorMap: Could not create empty raster");
4277 PG_FREE_IF_COPY(pgraster, 0);
4281 if (pgraster == NULL)
4284 SET_VARSIZE(pgraster, ((
rt_pgraster*) pgraster)->size);
4285 PG_RETURN_POINTER(pgraster);
4290 if (arg->
band == NULL) {
4293 PG_FREE_IF_COPY(pgraster, 0);
4294 elog(ERROR,
"RASTER_colorMap: Could not get band at index %d",
nband);
4299 if (!PG_ARGISNULL(3)) {
4300 char *method = NULL;
4301 char *tmp = text_to_cstring(PG_GETARG_TEXT_P(3));
4308 if (strcmp(method,
"INTERPOLATE") == 0)
4310 else if (strcmp(method,
"EXACT") == 0)
4312 else if (strcmp(method,
"NEAREST") == 0)
4315 elog(NOTICE,
"Unknown value provided for method. Defaulting to INTERPOLATE");
4325 if (PG_ARGISNULL(2)) {
4327 PG_FREE_IF_COPY(pgraster, 0);
4328 elog(ERROR,
"RASTER_colorMap: Value must be provided for colormap");
4333 char *colormap = text_to_cstring(PG_GETARG_TEXT_P(2));
4342 if (!strlen(colormap)) {
4344 PG_FREE_IF_COPY(pgraster, 0);
4345 elog(ERROR,
"RASTER_colorMap: Value must be provided for colormap");
4353 PG_FREE_IF_COPY(pgraster, 0);
4354 elog(ERROR,
"RASTER_colorMap: Could not process the value provided for colormap");
4362 PG_FREE_IF_COPY(pgraster, 0);
4363 elog(ERROR,
"RASTER_colorMap: Could not allocate memory for colormap entries");
4369 for (i = 0; i < arg->
nentry; i++) {
4383 if (!strlen(_entry)) {
4393 PG_FREE_IF_COPY(pgraster, 0);
4394 elog(ERROR,
"RASTER_colorMap: Could not process colormap entry %d", i + 1);
4398 elog(NOTICE,
"More than five elements in colormap entry %d. Using at most five elements", i + 1);
4407 for (j = 0; j < arg->
nelement; j++) {
4416 char *percent = NULL;
4420 strcmp(_element,
"NV") == 0 ||
4421 strcmp(_element,
"NULL") == 0 ||
4422 strcmp(_element,
"NODATA") == 0
4427 elog(NOTICE,
"More than one NODATA entry found. Using only the first one");
4436 else if ((percent = strchr(_element,
'%')) != NULL) {
4448 PG_FREE_IF_COPY(pgraster, 0);
4449 elog(ERROR,
"RASTER_colorMap: Could not get band's summary stats to process percentages");
4455 tmp = palloc(
sizeof(
char) * (percent - _element + 1));
4459 PG_FREE_IF_COPY(pgraster, 0);
4460 elog(ERROR,
"RASTER_colorMap: Could not allocate memory for value of percentage");
4464 memcpy(tmp, _element, percent - _element);
4465 tmp[percent - _element] =
'\0';
4470 value = strtod(tmp, NULL);
4472 if (errno != 0 || _element == junk) {
4475 PG_FREE_IF_COPY(pgraster, 0);
4476 elog(ERROR,
"RASTER_colorMap: Could not process percent string to value");
4482 elog(NOTICE,
"Percentage values cannot be less than zero. Defaulting to zero");
4485 else if (
value > 100.) {
4486 elog(NOTICE,
"Percentage values cannot be greater than 100. Defaulting to 100");
4498 if (errno != 0 || _element == junk) {
4501 PG_FREE_IF_COPY(pgraster, 0);
4502 elog(ERROR,
"RASTER_colorMap: Could not process string to value");
4513 value = (int) strtod(_element, &junk);
4514 if (errno != 0 || _element == junk) {
4517 PG_FREE_IF_COPY(pgraster, 0);
4518 elog(ERROR,
"RASTER_colorMap: Could not process string to value");
4523 elog(NOTICE,
"RGBA value cannot be greater than 255. Defaulting to 255");
4526 else if (
value < 0) {
4527 elog(NOTICE,
"RGBA value cannot be less than zero. Defaulting to zero");
4536 POSTGIS_RT_DEBUGF(4,
"colormap->entry[%d] (isnodata, value, R, G, B, A) = (%d, %f, %d, %d, %d, %d)",
4557 PG_FREE_IF_COPY(pgraster, 0);
4558 elog(ERROR,
"RASTER_colorMap: Could not create new raster with applied colormap");
4563 PG_FREE_IF_COPY(pgraster, 0);
4569 if (pgraster == NULL)
4572 SET_VARSIZE(pgraster, ((
rt_pgraster*) pgraster)->size);
4573 PG_RETURN_POINTER(pgraster);
4585 int x,
y,
nband, width, height;
4587 double newnodatavalue = 0.0;
4588 double newinitialvalue = 0.0;
4589 double newval = 0.0;
4590 char *newexpr = NULL;
4591 char *initexpr = NULL;
4592 char *expression = NULL;
4593 int hasnodataval = 0;
4594 double nodataval = 0.;
4596 int skipcomputation = 0;
4598 const int argkwcount = 3;
4599 enum KEYWORDS { kVAL=0, kX=1, kY=2 };
4600 char *argkw[] = {
"[rast]",
"[rast.x]",
"[rast.y]"};
4601 Oid argkwtypes[] = { FLOAT8OID, INT4OID, INT4OID };
4603 Oid argtype[] = { FLOAT8OID, INT4OID, INT4OID };
4604 uint8_t argpos[3] = {0};
4609 SPIPlanPtr spi_plan = NULL;
4610 SPITupleTable * tuptable = NULL;
4612 char * strFromText = NULL;
4613 Datum *values = NULL;
4614 Datum datum = (Datum)NULL;
4616 bool isnull =
FALSE;
4623 if (PG_ARGISNULL(0)) {
4624 elog(NOTICE,
"Raster is NULL. Returning NULL");
4630 pgraster = (
rt_pgraster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
4633 PG_FREE_IF_COPY(pgraster, 0);
4634 elog(ERROR,
"RASTER_mapAlgebraExpr: Could not deserialize raster");
4640 if (PG_ARGISNULL(1))
4643 nband = PG_GETARG_INT32(1);
4649 POSTGIS_RT_DEBUG(3,
"RASTER_mapAlgebraExpr: Creating new empty raster...");
4660 if ( NULL == newrast ) {
4661 PG_FREE_IF_COPY(pgraster, 0);
4662 elog(ERROR,
"RASTER_mapAlgebraExpr: Could not create a new raster");
4687 elog(NOTICE,
"Raster is empty. Returning an empty raster");
4689 PG_FREE_IF_COPY(pgraster, 0);
4693 if (NULL == pgrtn) {
4695 elog(ERROR,
"RASTER_mapAlgebraExpr: Could not serialize raster");
4699 SET_VARSIZE(pgrtn, pgrtn->
size);
4700 PG_RETURN_POINTER(pgrtn);
4711 elog(NOTICE,
"Raster does not have the required band. Returning a raster "
4714 PG_FREE_IF_COPY(pgraster, 0);
4718 if (NULL == pgrtn) {
4719 elog(ERROR,
"RASTER_mapAlgebraExpr: Could not serialize raster");
4723 SET_VARSIZE(pgrtn, pgrtn->
size);
4724 PG_RETURN_POINTER(pgrtn);
4729 if ( NULL ==
band ) {
4730 elog(NOTICE,
"Could not get the required band. Returning a raster "
4733 PG_FREE_IF_COPY(pgraster, 0);
4737 if (NULL == pgrtn) {
4738 elog(ERROR,
"RASTER_mapAlgebraExpr: Could not serialize raster");
4742 SET_VARSIZE(pgrtn, pgrtn->
size);
4743 PG_RETURN_POINTER(pgrtn);
4749 POSTGIS_RT_DEBUG(3,
"RASTER_mapAlgebraExpr: Getting NODATA value for band...");
4767 newinitialvalue = newnodatavalue;
4774 if (PG_ARGISNULL(2)) {
4779 strFromText = text_to_cstring(PG_GETARG_TEXT_P(2));
4782 if (newpixeltype ==
PT_END)
4786 if (newpixeltype ==
PT_END) {
4787 PG_FREE_IF_COPY(pgraster, 0);
4788 elog(ERROR,
"RASTER_mapAlgebraExpr: Invalid pixeltype");
4797 if (!PG_ARGISNULL(3)) {
4798 expression = text_to_cstring(PG_GETARG_TEXT_P(3));
4799 len = strlen(
"SELECT (") + strlen(expression) + strlen(
")::double precision");
4800 initexpr = (
char *)palloc(len + 1);
4802 memcpy(initexpr,
"SELECT (", strlen(
"SELECT ("));
4803 memcpy(initexpr + strlen(
"SELECT ("), expression, strlen(expression));
4804 memcpy(initexpr + strlen(
"SELECT (") + strlen(expression),
")::double precision", strlen(
")::double precision"));
4805 initexpr[len] =
'\0';
4823 if (!PG_ARGISNULL(4)) {
4825 nodataval = PG_GETARG_FLOAT8(4);
4826 newinitialvalue = nodataval;
4843 POSTGIS_RT_DEBUG(3,
"RASTER_mapAlgebraExpr: Band is a nodata band, returning "
4844 "a raster filled with nodata");
4847 newinitialvalue,
TRUE, newnodatavalue, 0);
4853 PG_FREE_IF_COPY(pgraster, 0);
4858 if (NULL == pgrtn) {
4859 elog(ERROR,
"RASTER_mapAlgebraExpr: Could not serialize raster");
4863 SET_VARSIZE(pgrtn, pgrtn->
size);
4864 PG_RETURN_POINTER(pgrtn);
4872 if (initexpr != NULL && ( !strcmp(initexpr,
"SELECT [rast]") || !strcmp(initexpr,
"SELECT [rast.val]") ) && !hasnodataval) {
4875 "Returning raster with band %d from original raster",
nband);
4887 PG_FREE_IF_COPY(pgraster, 0);
4892 if (NULL == pgrtn) {
4893 elog(ERROR,
"RASTER_mapAlgebraExpr: Could not serialize raster");
4897 SET_VARSIZE(pgrtn, pgrtn->
size);
4898 PG_RETURN_POINTER(pgrtn);
4905 if (initexpr != NULL && strstr(initexpr,
"[rast") == NULL) {
4906 ret = SPI_connect();
4907 if (ret != SPI_OK_CONNECT) {
4908 PG_FREE_IF_COPY(pgraster, 0);
4909 elog(ERROR,
"RASTER_mapAlgebraExpr: Could not connect to the SPI manager");
4914 ret = SPI_execute(initexpr,
FALSE, 0);
4916 if (ret != SPI_OK_SELECT || SPI_tuptable == NULL || SPI_processed != 1) {
4920 SPI_freetuptable(tuptable);
4921 PG_FREE_IF_COPY(pgraster, 0);
4924 elog(ERROR,
"RASTER_mapAlgebraExpr: Invalid construction for expression");
4928 tupdesc = SPI_tuptable->tupdesc;
4929 tuptable = SPI_tuptable;
4931 tuple = tuptable->vals[0];
4932 newexpr = SPI_getvalue(tuple, tupdesc, 1);
4934 POSTGIS_RT_DEBUG(3,
"Constant expression evaluated to NULL, keeping initvalue");
4935 newval = newinitialvalue;
4937 newval = atof(newexpr);
4940 SPI_freetuptable(tuptable);
4947 skipcomputation = 1;
4953 if (!hasnodataval) {
4954 newinitialvalue = newval;
4955 skipcomputation = 2;
4959 else if (
FLT_NEQ(newval, newinitialvalue)) {
4960 skipcomputation = 2;
4969 newinitialvalue,
TRUE, newnodatavalue, 0);
4976 if (expression == NULL || skipcomputation == 2) {
4982 PG_FREE_IF_COPY(pgraster, 0);
4987 if (NULL == pgrtn) {
4988 elog(ERROR,
"RASTER_mapAlgebraExpr: Could not serialize raster");
4992 SET_VARSIZE(pgrtn, pgrtn->
size);
4993 PG_RETURN_POINTER(pgrtn);
4996 RASTER_DEBUG(3,
"RASTER_mapAlgebraExpr: Creating new raster band...");
5000 if ( NULL == newband ) {
5001 elog(NOTICE,
"Could not modify band for new raster. Returning new "
5002 "raster with the original band");
5007 PG_FREE_IF_COPY(pgraster, 0);
5012 if (NULL == pgrtn) {
5013 elog(ERROR,
"RASTER_mapAlgebraExpr: Could not serialize raster");
5017 SET_VARSIZE(pgrtn, pgrtn->
size);
5018 PG_RETURN_POINTER(pgrtn);
5024 if (initexpr != NULL) {
5027 pfree(initexpr); initexpr=newexpr;
5029 sprintf(place,
"$1");
5030 for (i = 0, j = 1; i < argkwcount; i++) {
5033 pfree(initexpr); initexpr=newexpr;
5035 argtype[argcount] = argkwtypes[i];
5039 sprintf(place,
"$%d", j);
5049 values = (Datum *) palloc(
sizeof(Datum) * argcount);
5050 if (values == NULL) {
5055 PG_FREE_IF_COPY(pgraster, 0);
5058 elog(ERROR,
"RASTER_mapAlgebraExpr: Could not allocate memory for value parameters of prepared statement");
5063 nulls = (
char *)palloc(argcount);
5064 if (nulls == NULL) {
5069 PG_FREE_IF_COPY(pgraster, 0);
5072 elog(ERROR,
"RASTER_mapAlgebraExpr: Could not allocate memory for null parameters of prepared statement");
5077 ret = SPI_connect();
5078 if (ret != SPI_OK_CONNECT) {
5083 PG_FREE_IF_COPY(pgraster, 0);
5086 elog(ERROR,
"RASTER_mapAlgebraExpr: Could not connect to the SPI manager");
5091 spi_plan = SPI_prepare(initexpr, argcount, argtype);
5093 if (spi_plan == NULL) {
5096 PG_FREE_IF_COPY(pgraster, 0);
5103 elog(ERROR,
"RASTER_mapAlgebraExpr: Could not prepare expression");
5108 for (
x = 0;
x < width;
x++) {
5109 for(
y = 0;
y < height;
y++) {
5117 if (skipcomputation == 0) {
5118 if (initexpr != NULL) {
5120 memset(nulls,
'n', argcount);
5122 for (i = 0; i < argkwcount; i++) {
5124 if (idx < 1)
continue;
5129 values[idx] = Int32GetDatum(
x+1);
5134 values[idx] = Int32GetDatum(
y+1);
5137 else if (i == kVAL ) {
5138 values[idx] = Float8GetDatum(
r);
5144 ret = SPI_execute_plan(spi_plan, values, nulls,
FALSE, 0);
5145 if (ret != SPI_OK_SELECT || SPI_tuptable == NULL ||
5146 SPI_processed != 1) {
5148 SPI_freetuptable(tuptable);
5150 SPI_freeplan(spi_plan);
5158 PG_FREE_IF_COPY(pgraster, 0);
5161 elog(ERROR,
"RASTER_mapAlgebraExpr: Error executing prepared plan");
5166 tupdesc = SPI_tuptable->tupdesc;
5167 tuptable = SPI_tuptable;
5169 tuple = tuptable->vals[0];
5170 datum = SPI_getbinval(tuple, tupdesc, 1, &isnull);
5171 if ( SPI_result == SPI_ERROR_NOATTRIBUTE ) {
5172 POSTGIS_RT_DEBUGF(3,
"Expression for pixel %d,%d (value %g) errored, skip setting",
x+1,
y+1,
r);
5173 newval = newinitialvalue;
5175 else if ( isnull ) {
5176 POSTGIS_RT_DEBUGF(3,
"Expression for pixel %d,%d (value %g) evaluated to NULL, skip setting",
x+1,
y+1,
r);
5177 newval = newinitialvalue;
5179 newval = DatumGetFloat8(datum);
5182 SPI_freetuptable(tuptable);
5186 newval = newinitialvalue;
5199 if (initexpr != NULL) {
5200 SPI_freeplan(spi_plan);
5214 POSTGIS_RT_DEBUG(3,
"RASTER_mapAlgebraExpr: raster modified, serializing it.");
5218 PG_FREE_IF_COPY(pgraster, 0);
5225 SET_VARSIZE(pgrtn, pgrtn->
size);
5233 PG_RETURN_POINTER(pgrtn);
5248 int x,
y,
nband, width, height;
5250 double newnodatavalue = 0.0;
5251 double newinitialvalue = 0.0;
5252 double newval = 0.0;
5257 LOCAL_FCINFO(cbdata, FUNC_MAX_ARGS);
5260 char * strFromText = NULL;
5266 if (PG_ARGISNULL(0)) {
5267 elog(WARNING,
"Raster is NULL. Returning NULL");
5273 pgraster = (
rt_pgraster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
5276 PG_FREE_IF_COPY(pgraster, 0);
5277 elog(ERROR,
"RASTER_mapAlgebraFct: Could not deserialize raster");
5285 if (PG_ARGISNULL(1))
5288 nband = PG_GETARG_INT32(1);
5304 if ( NULL == newrast ) {
5307 PG_FREE_IF_COPY(pgraster, 0);
5309 elog(ERROR,
"RASTER_mapAlgebraFct: Could not create a new raster");
5334 elog(NOTICE,
"Raster is empty. Returning an empty raster");
5336 PG_FREE_IF_COPY(pgraster, 0);
5340 if (NULL == pgrtn) {
5341 elog(ERROR,
"RASTER_mapAlgebraFct: Could not serialize raster");
5345 SET_VARSIZE(pgrtn, pgrtn->
size);
5346 PG_RETURN_POINTER(pgrtn);
5356 elog(NOTICE,
"Raster does not have the required band. Returning a raster "
5359 PG_FREE_IF_COPY(pgraster, 0);
5363 if (NULL == pgrtn) {
5364 elog(ERROR,
"RASTER_mapAlgebraFct: Could not serialize raster");
5368 SET_VARSIZE(pgrtn, pgrtn->
size);
5369 PG_RETURN_POINTER(pgrtn);
5374 if ( NULL ==
band ) {
5375 elog(NOTICE,
"Could not get the required band. Returning a raster "
5378 PG_FREE_IF_COPY(pgraster, 0);
5382 if (NULL == pgrtn) {
5383 elog(ERROR,
"RASTER_mapAlgebraFct: Could not serialize raster");
5387 SET_VARSIZE(pgrtn, pgrtn->
size);
5388 PG_RETURN_POINTER(pgrtn);
5394 POSTGIS_RT_DEBUG(3,
"RASTER_mapAlgebraFct: Getting NODATA value for band...");
5411 newinitialvalue = newnodatavalue;
5418 if (PG_ARGISNULL(2)) {
5423 strFromText = text_to_cstring(PG_GETARG_TEXT_P(2));
5426 if (newpixeltype ==
PT_END)
5430 if (newpixeltype ==
PT_END) {
5433 PG_FREE_IF_COPY(pgraster, 0);
5436 elog(ERROR,
"RASTER_mapAlgebraFct: Invalid pixeltype");
5444 if (PG_ARGISNULL(3)) {
5447 PG_FREE_IF_COPY(pgraster, 0);
5450 elog(ERROR,
"RASTER_mapAlgebraFct: Required function is missing. Returning NULL");
5454 oid = PG_GETARG_OID(3);
5455 if (oid == InvalidOid) {
5458 PG_FREE_IF_COPY(pgraster, 0);
5461 elog(ERROR,
"RASTER_mapAlgebraFct: Got invalid function object id. Returning NULL");
5465 fmgr_info(oid, &cbinfo);
5468 if (cbinfo.fn_retset) {
5471 PG_FREE_IF_COPY(pgraster, 0);
5474 elog(ERROR,
"RASTER_mapAlgebraFct: Function provided must return double precision not resultset");
5478 else if (cbinfo.fn_nargs < 2 || cbinfo.fn_nargs > 3) {
5481 PG_FREE_IF_COPY(pgraster, 0);
5484 elog(ERROR,
"RASTER_mapAlgebraFct: Function does not have two or three input parameters");
5488 if (cbinfo.fn_nargs == 2)
5493 if (func_volatile(oid) ==
'v') {
5494 elog(NOTICE,
"Function provided is VOLATILE. Unless required and for best performance, function should be IMMUTABLE or STABLE");
5498 InitFunctionCallInfoData(*cbdata, &cbinfo, 2, InvalidOid, NULL, NULL);
5500 cbdata->args[0].isnull =
FALSE;
5501 cbdata->args[1].isnull =
FALSE;
5502 cbdata->args[2].isnull =
FALSE;
5505 if (PG_ARGISNULL(4)) {
5506 if (cbinfo.fn_strict) {
5509 PG_FREE_IF_COPY(pgraster, 0);
5512 elog(ERROR,
"RASTER_mapAlgebraFct: Strict callback functions cannot have null parameters");
5516 cbdata->args[k].value = (Datum)NULL;
5517 cbdata->args[k].isnull =
TRUE;
5520 cbdata->args[k].value = PG_GETARG_DATUM(4);
5530 POSTGIS_RT_DEBUG(3,
"RASTER_mapAlgebraFct: Band is a nodata band, returning "
5531 "a raster filled with nodata");
5534 newinitialvalue,
TRUE, newnodatavalue, 0);
5537 PG_FREE_IF_COPY(pgraster, 0);
5542 if (NULL == pgrtn) {
5543 elog(ERROR,
"RASTER_mapAlgebraFct: Could not serialize raster");
5547 SET_VARSIZE(pgrtn, pgrtn->
size);
5548 PG_RETURN_POINTER(pgrtn);
5557 newinitialvalue,
TRUE, newnodatavalue, 0);
5561 if ( NULL == newband ) {
5562 elog(NOTICE,
"Could not modify band for new raster. Returning new "
5563 "raster with the original band");
5566 PG_FREE_IF_COPY(pgraster, 0);
5571 if (NULL == pgrtn) {
5572 elog(ERROR,
"RASTER_mapAlgebraFct: Could not serialize raster");
5576 SET_VARSIZE(pgrtn, pgrtn->
size);
5577 PG_RETURN_POINTER(pgrtn);
5583 for (
x = 0;
x < width;
x++) {
5584 for(
y = 0;
y < height;
y++) {
5592 if (
FLT_EQ(
r, newnodatavalue)) {
5593 if (cbinfo.fn_strict) {
5594 POSTGIS_RT_DEBUG(3,
"RASTER_mapAlgebraFct: Strict callbacks cannot accept NULL arguments, skipping NODATA cell.");
5597 cbdata->args[0].isnull =
TRUE;
5598 cbdata->args[0].value = (Datum)NULL;
5601 cbdata->args[0].isnull =
FALSE;
5602 cbdata->args[0].value = Float8GetDatum(
r);
5606 if (cbinfo.fn_nargs == 3) {
5610 d[0] = Int32GetDatum(
x+1);
5611 d[1] = Int32GetDatum(
y+1);
5613 a = construct_array(d, 2, INT4OID,
sizeof(
int32),
true,
'i');
5615 cbdata->args[1].isnull =
FALSE;
5616 cbdata->args[1].value = PointerGetDatum(a);
5622 tmpnewval = FunctionCallInvoke(cbdata);
5626 newval = newnodatavalue;
5629 newval = DatumGetFloat8(tmpnewval);
5643 POSTGIS_RT_DEBUG(3,
"RASTER_mapAlgebraFct: raster modified, serializing it.");
5647 PG_FREE_IF_COPY(pgraster, 0);
5658 SET_VARSIZE(pgrtn, pgrtn->
size);
5659 PG_RETURN_POINTER(pgrtn);
5674 int x,
y,
nband, width, height, ngbwidth, ngbheight, winwidth, winheight, u, v, nIndex, nNullItems;
5676 double newnodatavalue = 0.0;
5677 double newinitialvalue = 0.0;
5678 double newval = 0.0;
5683 LOCAL_FCINFO(cbdata, FUNC_MAX_ARGS);
5685 ArrayType * neighborDatum;
5686 char * strFromText = NULL;
5687 text * txtNodataMode = NULL;
5688 text * txtCallbackParam = NULL;
5690 float fltReplace = 0;
5691 bool valuereplace =
false, pixelreplace, nNodataOnly =
true, nNullSkip =
false;
5692 Datum * neighborData = NULL;
5693 bool * neighborNulls = NULL;
5694 int neighborDims[2];
5703 if (PG_ARGISNULL(0)) {
5704 elog(WARNING,
"Raster is NULL. Returning NULL");
5710 pgraster = (
rt_pgraster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
5714 PG_FREE_IF_COPY(pgraster, 0);
5715 elog(ERROR,
"RASTER_mapAlgebraFctNgb: Could not deserialize raster");
5723 if (PG_ARGISNULL(1))
5726 nband = PG_GETARG_INT32(1);
5731 POSTGIS_RT_DEBUG(3,
"RASTER_mapAlgebraFctNgb: Creating new empty raster...");
5742 if ( NULL == newrast ) {
5744 PG_FREE_IF_COPY(pgraster, 0);
5745 elog(ERROR,
"RASTER_mapAlgebraFctNgb: Could not create a new raster");
5770 elog(NOTICE,
"Raster is empty. Returning an empty raster");
5772 PG_FREE_IF_COPY(pgraster, 0);
5776 if (NULL == pgrtn) {
5777 elog(ERROR,
"RASTER_mapAlgebraFctNgb: Could not serialize raster");
5781 SET_VARSIZE(pgrtn, pgrtn->
size);
5782 PG_RETURN_POINTER(pgrtn);
5792 elog(NOTICE,
"Raster does not have the required band. Returning a raster "
5795 PG_FREE_IF_COPY(pgraster, 0);
5799 if (NULL == pgrtn) {
5800 elog(ERROR,
"RASTER_mapAlgebraFctNgb: Could not serialize raster");
5804 SET_VARSIZE(pgrtn, pgrtn->
size);
5805 PG_RETURN_POINTER(pgrtn);
5810 if ( NULL ==
band ) {
5811 elog(NOTICE,
"Could not get the required band. Returning a raster "
5814 PG_FREE_IF_COPY(pgraster, 0);
5818 if (NULL == pgrtn) {
5819 elog(ERROR,
"RASTER_mapAlgebraFctNgb: Could not serialize raster");
5823 SET_VARSIZE(pgrtn, pgrtn->
size);
5824 PG_RETURN_POINTER(pgrtn);
5830 POSTGIS_RT_DEBUG(3,
"RASTER_mapAlgebraFctNgb: Getting NODATA value for band...");
5847 newinitialvalue = newnodatavalue;
5854 if (PG_ARGISNULL(2)) {
5859 strFromText = text_to_cstring(PG_GETARG_TEXT_P(2));
5860 POSTGIS_RT_DEBUGF(3,
"RASTER_mapAlgebraFctNgb: Pixeltype parameter: %s", strFromText);
5863 if (newpixeltype ==
PT_END)
5867 if (newpixeltype ==
PT_END) {
5870 PG_FREE_IF_COPY(pgraster, 0);
5873 elog(ERROR,
"RASTER_mapAlgebraFctNgb: Invalid pixeltype");
5881 if (PG_ARGISNULL(5)) {
5884 PG_FREE_IF_COPY(pgraster, 0);
5887 elog(ERROR,
"RASTER_mapAlgebraFctNgb: Required function is missing");
5891 oid = PG_GETARG_OID(5);
5892 if (oid == InvalidOid) {
5895 PG_FREE_IF_COPY(pgraster, 0);
5898 elog(ERROR,
"RASTER_mapAlgebraFctNgb: Got invalid function object id");
5902 fmgr_info(oid, &cbinfo);
5905 if (cbinfo.fn_retset) {
5908 PG_FREE_IF_COPY(pgraster, 0);
5911 elog(ERROR,
"RASTER_mapAlgebraFctNgb: Function provided must return double precision not resultset");
5915 else if (cbinfo.fn_nargs != 3) {
5918 PG_FREE_IF_COPY(pgraster, 0);
5921 elog(ERROR,
"RASTER_mapAlgebraFctNgb: Function does not have three input parameters");
5925 if (func_volatile(oid) ==
'v') {
5926 elog(NOTICE,
"Function provided is VOLATILE. Unless required and for best performance, function should be IMMUTABLE or STABLE");
5930 InitFunctionCallInfoData(*cbdata, &cbinfo, 3, InvalidOid, NULL, NULL);
5931 cbdata->args[0].isnull =
FALSE;
5932 cbdata->args[1].isnull =
FALSE;
5933 cbdata->args[2].isnull =
FALSE;
5936 if (PG_ARGISNULL(7)) {
5937 if (cbinfo.fn_strict) {
5940 PG_FREE_IF_COPY(pgraster, 0);
5943 elog(ERROR,
"RASTER_mapAlgebraFctNgb: Strict callback functions cannot have NULL parameters");
5947 cbdata->args[2].value = (Datum)NULL;
5948 cbdata->args[2].isnull =
TRUE;
5951 cbdata->args[2].value = PG_GETARG_DATUM(7);
5961 POSTGIS_RT_DEBUG(3,
"RASTER_mapAlgebraFctNgb: Band is a nodata band, returning "
5962 "a raster filled with nodata");
5965 newinitialvalue,
TRUE, newnodatavalue, 0);
5968 PG_FREE_IF_COPY(pgraster, 0);
5973 if (NULL == pgrtn) {
5974 elog(ERROR,
"RASTER_mapAlgebraFctNgb: Could not serialize raster");
5978 SET_VARSIZE(pgrtn, pgrtn->
size);
5979 PG_RETURN_POINTER(pgrtn);
5988 newinitialvalue,
TRUE, newnodatavalue, 0);
5992 if ( NULL == newband ) {
5993 elog(NOTICE,
"Could not modify band for new raster. Returning new "
5994 "raster with the original band");
5997 PG_FREE_IF_COPY(pgraster, 0);
6002 if (NULL == pgrtn) {
6003 elog(ERROR,
"RASTER_mapAlgebraFctNgb: Could not serialize raster");
6007 SET_VARSIZE(pgrtn, pgrtn->
size);
6008 PG_RETURN_POINTER(pgrtn);
6012 if (PG_ARGISNULL(3) || PG_GETARG_INT32(3) <= 0) {
6013 elog(NOTICE,
"Neighborhood width is NULL or <= 0. Returning new "
6014 "raster with the original band");
6017 PG_FREE_IF_COPY(pgraster, 0);
6022 if (NULL == pgrtn) {
6023 elog(ERROR,
"RASTER_mapAlgebraFctNgb: Could not serialize raster");
6027 SET_VARSIZE(pgrtn, pgrtn->
size);
6028 PG_RETURN_POINTER(pgrtn);
6031 ngbwidth = PG_GETARG_INT32(3);
6032 winwidth = ngbwidth * 2 + 1;
6035 if (PG_ARGISNULL(4) || PG_GETARG_INT32(4) <= 0) {
6036 elog(NOTICE,
"Neighborhood height is NULL or <= 0. Returning new "
6037 "raster with the original band");
6040 PG_FREE_IF_COPY(pgraster, 0);
6045 if (NULL == pgrtn) {
6046 elog(ERROR,
"RASTER_mapAlgebraFctNgb: Could not serialize raster");
6050 SET_VARSIZE(pgrtn, pgrtn->
size);
6051 PG_RETURN_POINTER(pgrtn);
6054 ngbheight = PG_GETARG_INT32(4);
6055 winheight = ngbheight * 2 + 1;
6058 if (PG_ARGISNULL(6)) {
6059 elog(NOTICE,
"Neighborhood NODATA behavior defaulting to 'ignore'");
6060 txtNodataMode = cstring_to_text(
"ignore");
6063 txtNodataMode = PG_GETARG_TEXT_P(6);
6066 txtCallbackParam = (text*)palloc(VARSIZE(txtNodataMode));
6067 SET_VARSIZE(txtCallbackParam, VARSIZE(txtNodataMode));
6068 memcpy((
void *)VARDATA(txtCallbackParam), (
void *)VARDATA(txtNodataMode), VARSIZE(txtNodataMode) - VARHDRSZ);
6071 cbdata->args[1].value = PointerGetDatum(txtCallbackParam);
6073 strFromText = text_to_cstring(txtNodataMode);
6076 if (strcmp(strFromText,
"VALUE") == 0)
6077 valuereplace =
true;
6078 else if (strcmp(strFromText,
"IGNORE") != 0 && strcmp(strFromText,
"NULL") != 0) {
6080 if (sscanf(strFromText,
"%d", &intReplace) <= 0 && sscanf(strFromText,
"%f", &fltReplace) <= 0) {
6082 elog(NOTICE,
"Neighborhood NODATA mode is not recognized. Must be one of 'value', 'ignore', "
6083 "'NULL', or a numeric value. Returning new raster with the original band");
6086 pfree(txtCallbackParam);
6090 PG_FREE_IF_COPY(pgraster, 0);
6095 if (NULL == pgrtn) {
6096 elog(ERROR,
"RASTER_mapAlgebraFctNgb: Could not serialize raster");
6100 SET_VARSIZE(pgrtn, pgrtn->
size);
6101 PG_RETURN_POINTER(pgrtn);
6104 else if (strcmp(strFromText,
"NULL") == 0) {
6113 neighborData = (Datum *)palloc(
sizeof(Datum) * winwidth * winheight);
6114 neighborNulls = (
bool *)palloc(
sizeof(
bool) * winwidth * winheight);
6117 neighborDims[0] = winwidth;
6118 neighborDims[1] = winheight;
6125 get_typlenbyvalalign(FLOAT8OID, &typlen, &typbyval, &typalign);
6127 for (
x = 0 + ngbwidth;
x < width - ngbwidth;
x++) {
6128 for(
y = 0 + ngbheight;
y < height - ngbheight;
y++) {
6133 pixelreplace =
false;
6137 pixelreplace =
true;
6140 for (u =
x - ngbwidth; u <=
x + ngbwidth; u++) {
6141 for (v =
y - ngbheight; v <=
y + ngbheight; v++) {
6146 neighborData[nIndex] = Float8GetDatum((
double)
r);
6147 neighborNulls[nIndex] =
false;
6148 nNodataOnly =
false;
6152 if (valuereplace && pixelreplace) {
6154 neighborData[nIndex] = Float8GetDatum((
double)rpix);
6155 neighborNulls[nIndex] =
false;
6160 neighborData[nIndex] = PointerGetDatum(NULL);
6161 neighborNulls[nIndex] =
true;
6168 neighborData[nIndex] = PointerGetDatum(NULL);
6169 neighborNulls[nIndex] =
true;
6181 if (!(nNodataOnly ||
6182 (nNullSkip && nNullItems > 0) ||
6183 (valuereplace && nNullItems > 0))) {
6185 x,
y, winwidth, winheight);
6187 neighborDatum = construct_md_array((
void *)neighborData, neighborNulls, 2, neighborDims, neighborLbs,
6188 FLOAT8OID, typlen, typbyval, typalign);
6191 cbdata->args[0].value = PointerGetDatum(neighborDatum);
6194 tmpnewval = FunctionCallInvoke(cbdata);
6199 newval = newnodatavalue;
6202 newval = DatumGetFloat8(tmpnewval);
6218 pfree(neighborNulls);
6219 pfree(neighborData);
6221 pfree(txtCallbackParam);
6224 PG_FREE_IF_COPY(pgraster, 0);
6228 POSTGIS_RT_DEBUG(3,
"RASTER_mapAlgebraFctNgb: raster modified, serializing it.");
6239 SET_VARSIZE(pgrtn, pgrtn->
size);
6240 PG_RETURN_POINTER(pgrtn);
6243 #define ARGKWCOUNT 8
6251 const uint32_t set_count = 2;
6253 int pgrastpos[2] = {-1, -1};
6256 int _isempty[2] = {0};
6257 uint32_t bandindex[2] = {0};
6260 int _hasnodata[2] = {0};
6261 double _nodataval[2] = {0};
6262 double _offset[4] = {0.};
6263 double _rastoffset[2][4] = {{0.}};
6264 int _haspixel[2] = {0};
6265 double _pixel[2] = {0};
6266 int _pos[2][2] = {{0}};
6267 uint16_t _dim[2][2] = {{0}};
6269 char *pixtypename = NULL;
6271 char *extenttypename = NULL;
6276 uint16_t dim[2] = {0};
6279 double nodataval = 0;
6280 double gt[6] = {0.};
6282 Oid calltype = InvalidOid;
6284 const uint32_t spi_count = 3;
6285 uint16_t spi_exprpos[3] = {4, 7, 8};
6286 uint32_t spi_argcount[3] = {0};
6289 SPIPlanPtr spi_plan[3] = {NULL};
6290 uint16_t spi_empty = 0;
6291 Oid *argtype = NULL;
6292 uint8_t argpos[3][8] = {{0}};
6293 char *argkw[] = {
"[rast1.x]",
"[rast1.y]",
"[rast1.val]",
"[rast1]",
"[rast2.x]",
"[rast2.y]",
"[rast2.val]",
"[rast2]"};
6297 SPITupleTable *tuptable = NULL;
6300 bool isnull =
FALSE;
6301 int hasargval[3] = {0};
6302 double argval[3] = {0.};
6303 int hasnodatanodataval = 0;
6304 double nodatanodataval = 0;
6307 Oid ufc_noid = InvalidOid;
6309 LOCAL_FCINFO(ufc_info, FUNC_MAX_ARGS);
6311 int ufc_nullcount = 0;
6327 for (i = 0, j = 0; i < set_count; i++) {
6328 if (!PG_ARGISNULL(j)) {
6329 pgrast[i] = (
rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(j));
6336 for (k = 0; k <= i; k++) {
6337 if (k < i &&
rast[k] != NULL)
6339 if (pgrastpos[k] != -1)
6340 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6342 elog(ERROR,
"RASTER_mapAlgebra2: Could not deserialize the %s raster", i < 1 ?
"first" :
"second");
6350 if (!PG_ARGISNULL(j)) {
6351 bandindex[i] = PG_GETARG_INT32(j);
6364 if (
rast[0] == NULL &&
rast[1] == NULL) {
6365 elog(NOTICE,
"The two rasters provided are NULL. Returning NULL");
6366 for (k = 0; k < set_count; k++) {
6367 if (pgrastpos[k] != -1)
6368 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6374 if (_isempty[0] && _isempty[1]) {
6375 elog(NOTICE,
"The two rasters provided are empty. Returning empty raster");
6379 for (k = 0; k < set_count; k++) {
6380 if (
rast[k] != NULL)
6382 if (pgrastpos[k] != -1)
6383 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6385 elog(ERROR,
"RASTER_mapAlgebra2: Could not create empty raster");
6395 SET_VARSIZE(pgrtn, pgrtn->
size);
6396 PG_RETURN_POINTER(pgrtn);
6401 (
rast[0] == NULL || _isempty[0]) ||
6402 (
rast[1] == NULL || _isempty[1])
6405 if (
rast[0] == NULL || _isempty[0]) {
6418 if (_rast[i] != NULL)
6430 if (_rast[i] == NULL) {
6432 for (k = 0; k < set_count; k++) {
6433 if (pgrastpos[k] != -1)
6434 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6436 elog(ERROR,
"RASTER_mapAlgebra2: Could not create NODATA raster");
6470 for (k = 0; k < set_count; k++) {
6471 if (_rast[k] != NULL)
6473 if (pgrastpos[k] != -1)
6474 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6476 elog(ERROR,
"RASTER_mapAlgebra2: Could not test for alignment on the two rasters");
6480 elog(NOTICE,
"The two rasters provided do not have the same alignment. Returning NULL");
6481 for (k = 0; k < set_count; k++) {
6482 if (_rast[k] != NULL)
6484 if (pgrastpos[k] != -1)
6485 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6491 if (!PG_ARGISNULL(5)) {
6492 pixtypename = text_to_cstring(PG_GETARG_TEXT_P(5));
6495 if (pixtype ==
PT_END ) {
6496 for (k = 0; k < set_count; k++) {
6497 if (_rast[k] != NULL)
6499 if (pgrastpos[k] != -1)
6500 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6502 elog(ERROR,
"RASTER_mapAlgebra2: Invalid pixel type: %s", pixtypename);
6508 if (!PG_ARGISNULL(6)) {
6521 for (k = 0; k < set_count; k++) {
6522 if (_rast[k] != NULL)
6524 if (pgrastpos[k] != -1)
6525 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6527 elog(ERROR,
"RASTER_mapAlgebra2: Could not get output raster of correct extent");
6532 _rastoffset[0][0] = _offset[0];
6533 _rastoffset[0][1] = _offset[1];
6534 _rastoffset[1][0] = _offset[2];
6535 _rastoffset[1][1] = _offset[3];
6543 switch (extenttype) {
6553 (extenttype ==
ET_FIRST && i == 0) ||
6557 elog(NOTICE,
"The %s raster is NULL. Returning NULL", (i != 1 ?
"FIRST" :
"SECOND"));
6558 for (k = 0; k < set_count; k++) {
6559 if (_rast[k] != NULL)
6561 if (pgrastpos[k] != -1)
6562 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6570 elog(NOTICE,
"The %s raster does not have the band at index %d. Returning no band raster of correct extent",
6571 (i != 1 ?
"FIRST" :
"SECOND"), bandindex[i]
6574 for (k = 0; k < set_count; k++) {
6575 if (_rast[k] != NULL)
6577 if (pgrastpos[k] != -1)
6578 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6583 if (!pgrtn) PG_RETURN_NULL();
6585 SET_VARSIZE(pgrtn, pgrtn->
size);
6586 PG_RETURN_POINTER(pgrtn);
6594 _isempty[0] || _isempty[1] ||
6597 elog(NOTICE,
"The two rasters provided have no intersection. Returning no band raster");
6600 if (dim[0] || dim[1]) {
6605 for (k = 0; k < set_count; k++) {
6606 if (_rast[k] != NULL)
6608 if (pgrastpos[k] != -1)
6609 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6611 elog(ERROR,
"RASTER_mapAlgebra2: Could not create no band raster");
6619 for (k = 0; k < set_count; k++) {
6620 if (_rast[k] != NULL)
6622 if (pgrastpos[k] != -1)
6623 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6628 if (!pgrtn) PG_RETURN_NULL();
6630 SET_VARSIZE(pgrtn, pgrtn->
size);
6631 PG_RETURN_POINTER(pgrtn);
6636 for (k = 0; k < set_count; k++) {
6637 if (_rast[k] != NULL)
6639 if (pgrastpos[k] != -1)
6640 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6642 elog(ERROR,
"RASTER_mapAlgebra2: ET_LAST and ET_CUSTOM are not implemented");
6652 elog(NOTICE,
"The two rasters provided do not have the respectively specified band indices. Returning no band raster of correct extent");
6654 for (k = 0; k < set_count; k++) {
6655 if (_rast[k] != NULL)
6657 if (pgrastpos[k] != -1)
6658 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6663 if (!pgrtn) PG_RETURN_NULL();
6665 SET_VARSIZE(pgrtn, pgrtn->
size);
6666 PG_RETURN_POINTER(pgrtn);
6670 for (i = 0; i < set_count; i++) {
6679 if (_band[i] == NULL) {
6680 for (k = 0; k < set_count; k++) {
6681 if (_rast[k] != NULL)
6683 if (pgrastpos[k] != -1)
6684 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6687 elog(ERROR,
"RASTER_mapAlgebra2: Could not get band %d of the %s raster",
6689 (i < 1 ?
"FIRST" :
"SECOND")
6701 if ((extenttype ==
ET_SECOND && !_isempty[1]) || _isempty[0])
6708 if (extenttype ==
ET_SECOND && !_isempty[1] && _hasnodata[1]) {
6709 nodataval = _nodataval[1];
6711 else if (!_isempty[0] && _hasnodata[0]) {
6712 nodataval = _nodataval[0];
6714 else if (!_isempty[1] && _hasnodata[1]) {
6715 nodataval = _nodataval[1];
6718 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));
6730 for (k = 0; k < set_count; k++) {
6731 if (_rast[k] != NULL)
6733 if (pgrastpos[k] != -1)
6734 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6737 elog(ERROR,
"RASTER_mapAlgebra2: Could not add new band to output raster");
6744 for (k = 0; k < set_count; k++) {
6745 if (_rast[k] != NULL)
6747 if (pgrastpos[k] != -1)
6748 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6751 elog(ERROR,
"RASTER_mapAlgebra2: Could not get newly added band of output raster");
6756 (
int) _rastoffset[0][0],
6757 (
int) _rastoffset[0][1],
6758 (
int) _rastoffset[1][0],
6759 (
int) _rastoffset[1][1]
6762 POSTGIS_RT_DEBUGF(4,
"metadata = (%f, %f, %d, %d, %f, %f, %f, %f, %d)",
6779 calltype = get_fn_expr_argtype(fcinfo->flinfo, 4);
6786 if (SPI_connect() != SPI_OK_CONNECT) {
6787 for (k = 0; k < set_count; k++) {
6788 if (_rast[k] != NULL)
6790 if (pgrastpos[k] != -1)
6791 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6794 elog(ERROR,
"RASTER_mapAlgebra2: Could not connect to the SPI manager");
6799 memset(hasargval, 0,
sizeof(
int) * spi_count);
6809 for (i = 0; i < spi_count; i++) {
6810 if (!PG_ARGISNULL(spi_exprpos[i])) {
6812 char place[5] =
"$1";
6813 expr = text_to_cstring(PG_GETARG_TEXT_P(spi_exprpos[i]));
6827 sprintf(place,
"$%d", k);
6833 len = strlen(
"SELECT (") + strlen(expr) + strlen(
")::double precision");
6834 sql = (
char *) palloc(len + 1);
6837 for (k = 0; k < spi_count; k++) SPI_freeplan(spi_plan[k]);
6840 for (k = 0; k < set_count; k++) {
6841 if (_rast[k] != NULL)
6843 if (pgrastpos[k] != -1)
6844 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6848 elog(ERROR,
"RASTER_mapAlgebra2: Could not allocate memory for expression parameter %d", spi_exprpos[i]);
6852 memcpy(
sql,
"SELECT (", strlen(
"SELECT ("));
6853 memcpy(
sql + strlen(
"SELECT ("), expr, strlen(expr));
6854 memcpy(
sql + strlen(
"SELECT (") + strlen(expr),
")::double precision", strlen(
")::double precision"));
6860 if (spi_argcount[i]) {
6861 argtype = (Oid *) palloc(spi_argcount[i] *
sizeof(Oid));
6862 if (argtype == NULL) {
6865 for (k = 0; k < spi_count; k++) SPI_freeplan(spi_plan[k]);
6868 for (k = 0; k < set_count; k++) {
6869 if (_rast[k] != NULL)
6871 if (pgrastpos[k] != -1)
6872 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6876 elog(ERROR,
"RASTER_mapAlgebra2: Could not allocate memory for prepared plan argtypes of expression parameter %d", spi_exprpos[i]);
6882 if (argpos[i][j] < 1)
continue;
6886 (strstr(argkw[j],
"[rast1.x]") != NULL) ||
6887 (strstr(argkw[j],
"[rast1.y]") != NULL) ||
6888 (strstr(argkw[j],
"[rast2.x]") != NULL) ||
6889 (strstr(argkw[j],
"[rast2.y]") != NULL)
6891 argtype[k] = INT4OID;
6895 argtype[k] = FLOAT8OID;
6901 spi_plan[i] = SPI_prepare(
sql, spi_argcount[i], argtype);
6904 if (spi_plan[i] == NULL) {
6907 for (k = 0; k < spi_count; k++) SPI_freeplan(spi_plan[k]);
6910 for (k = 0; k < set_count; k++) {
6911 if (_rast[k] != NULL)
6913 if (pgrastpos[k] != -1)
6914 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6918 elog(ERROR,
"RASTER_mapAlgebra2: Could not create prepared plan of expression parameter %d", spi_exprpos[i]);
6924 err = SPI_execute(
sql,
TRUE, 0);
6925 if (err != SPI_OK_SELECT || SPI_tuptable == NULL || SPI_processed != 1) {
6928 for (k = 0; k < spi_count; k++) SPI_freeplan(spi_plan[k]);
6931 for (k = 0; k < set_count; k++) {
6932 if (_rast[k] != NULL)
6934 if (pgrastpos[k] != -1)
6935 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
6939 elog(ERROR,
"RASTER_mapAlgebra2: Could not evaluate expression parameter %d", spi_exprpos[i]);
6944 tupdesc = SPI_tuptable->tupdesc;
6945 tuptable = SPI_tuptable;
6946 tuple = tuptable->vals[0];
6948 datum = SPI_getbinval(tuple, tupdesc, 1, &isnull);
6949 if (SPI_result == SPI_ERROR_NOATTRIBUTE) {
6952 if (SPI_tuptable) SPI_freetuptable(tuptable);
6953 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 result of expression parameter %d", spi_exprpos[i]);
6970 argval[i] = DatumGetFloat8(datum);
6973 if (SPI_tuptable) SPI_freetuptable(tuptable);
6983 if (!PG_ARGISNULL(9)) {
6984 hasnodatanodataval = 1;
6985 nodatanodataval = PG_GETARG_FLOAT8(9);
6988 hasnodatanodataval = 0;
6991 case REGPROCEDUREOID: {
6993 if (!PG_ARGISNULL(4)) {
6996 ufc_noid = PG_GETARG_OID(4);
6999 fmgr_info(ufc_noid, &ufl_info);
7003 if (ufl_info.fn_retset) {
7007 else if (ufl_info.fn_nargs < 3 || ufl_info.fn_nargs > 4) {
7017 for (k = 0; k < set_count; k++) {
7018 if (_rast[k] != NULL)
7020 if (pgrastpos[k] != -1)
7021 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
7026 elog(ERROR,
"RASTER_mapAlgebra2: Function provided must have three or four input parameters");
7028 elog(ERROR,
"RASTER_mapAlgebra2: Function provided must return double precision not resultset");
7032 if (func_volatile(ufc_noid) ==
'v') {
7033 elog(NOTICE,
"Function provided is VOLATILE. Unless required and for best performance, function should be IMMUTABLE or STABLE");
7037 InitFunctionCallInfoData(
7038 *ufc_info, &ufl_info, ufl_info.fn_nargs, InvalidOid, NULL, NULL);
7039 ufc_info->args[0].isnull =
FALSE;
7040 ufc_info->args[1].isnull =
FALSE;
7041 ufc_info->args[2].isnull =
FALSE;
7042 if (ufl_info.fn_nargs == 4)
7043 ufc_info->args[3].isnull =
FALSE;
7045 if (ufl_info.fn_nargs != 4)
7049 if (!PG_ARGISNULL(7))
7051 ufc_info->args[k].value = PG_GETARG_DATUM(7);
7055 ufc_info->args[k].value = (Datum)NULL;
7056 ufc_info->args[k].isnull =
TRUE;
7063 for (k = 0; k < set_count; k++) {
7064 if (_rast[k] != NULL)
7066 if (pgrastpos[k] != -1)
7067 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
7070 elog(ERROR,
"RASTER_mapAlgebra2: Invalid data type for expression or userfunction");
7078 (calltype == TEXTOID) && (
7079 (spi_empty != spi_count) || hasnodatanodataval
7082 (calltype == REGPROCEDUREOID) && (ufc_noid != InvalidOid)
7084 for (
x = 0;
x < dim[0];
x++) {
7085 for (
y = 0;
y < dim[1];
y++) {
7088 for (i = 0; i < set_count; i++) {
7093 _x = (int)
x - (
int)_rastoffset[i][0];
7094 _y = (int)
y - (
int)_rastoffset[i][1];
7097 _pos[i][0] = _x + 1;
7098 _pos[i][1] = _y + 1;
7101 if (_band[i] == NULL) {
7102 if (!_hasnodata[i]) {
7104 _pixel[i] = _nodataval[i];
7109 (_x >= 0 && _x < _dim[i][0]) &&
7110 (_y >= 0 && _y < _dim[i][1])
7115 if (calltype == TEXTOID) {
7116 for (k = 0; k < spi_count; k++) SPI_freeplan(spi_plan[k]);
7120 for (k = 0; k < set_count; k++) {
7121 if (_rast[k] != NULL)
7123 if (pgrastpos[k] != -1)
7124 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
7128 elog(ERROR,
"RASTER_mapAlgebra2: Could not get pixel of %s raster", (i < 1 ?
"FIRST" :
"SECOND"));
7132 if (!_hasnodata[i] || !isnodata)
7151 if (!_haspixel[0] && !_haspixel[1])
7155 else if (_haspixel[0] && !_haspixel[1])
7159 else if (!_haspixel[0] && _haspixel[1])
7168 if (hasnodatanodataval) {
7170 pixel = nodatanodataval;
7174 else if (hasargval[i]) {
7179 else if (spi_plan[i] != NULL) {
7183 memset(values, (Datum) NULL,
sizeof(Datum) *
ARGKWCOUNT);
7188 if (spi_argcount[i]) {
7192 if (idx < 1)
continue;
7195 if (strstr(argkw[j],
"[rast1.x]") != NULL) {
7196 values[idx] = _pos[0][0];
7198 else if (strstr(argkw[j],
"[rast1.y]") != NULL) {
7199 values[idx] = _pos[0][1];
7202 (strstr(argkw[j],
"[rast1.val]") != NULL) ||
7203 (strstr(argkw[j],
"[rast1]") != NULL)
7205 if (_isempty[0] || !_haspixel[0])
7208 values[idx] = Float8GetDatum(_pixel[0]);
7210 else if (strstr(argkw[j],
"[rast2.x]") != NULL) {
7211 values[idx] = _pos[1][0];
7213 else if (strstr(argkw[j],
"[rast2.y]") != NULL) {
7214 values[idx] = _pos[1][1];
7217 (strstr(argkw[j],
"[rast2.val]") != NULL) ||
7218 (strstr(argkw[j],
"[rast2]") != NULL)
7220 if (_isempty[1] || !_haspixel[1])
7223 values[idx] = Float8GetDatum(_pixel[1]);
7229 err = SPI_execute_plan(spi_plan[i], values, nulls,
TRUE, 1);
7230 if (err != SPI_OK_SELECT || SPI_tuptable == NULL || SPI_processed != 1) {
7232 for (k = 0; k < spi_count; k++) SPI_freeplan(spi_plan[k]);
7235 for (k = 0; k < set_count; k++) {
7236 if (_rast[k] != NULL)
7238 if (pgrastpos[k] != -1)
7239 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
7243 elog(ERROR,
"RASTER_mapAlgebra2: Unexpected error when running prepared statement %d", i);
7248 tupdesc = SPI_tuptable->tupdesc;
7249 tuptable = SPI_tuptable;
7250 tuple = tuptable->vals[0];
7252 datum = SPI_getbinval(tuple, tupdesc, 1, &isnull);
7253 if (SPI_result == SPI_ERROR_NOATTRIBUTE) {
7255 if (SPI_tuptable) SPI_freetuptable(tuptable);
7256 for (k = 0; k < spi_count; k++) SPI_freeplan(spi_plan[k]);
7259 for (k = 0; k < set_count; k++) {
7260 if (_rast[k] != NULL)
7262 if (pgrastpos[k] != -1)
7263 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
7267 elog(ERROR,
"RASTER_mapAlgebra2: Could not get result of prepared statement %d", i);
7273 pixel = DatumGetFloat8(datum);
7276 if (SPI_tuptable) SPI_freetuptable(tuptable);
7279 case REGPROCEDUREOID: {
7284 for (i = 0; i < set_count; i++) {
7285 ufc_info->args[i].value = Float8GetDatum(_pixel[i]);
7287 ufc_info->args[i].isnull =
FALSE;
7291 ufc_info->args[i].isnull =
TRUE;
7298 if (ufl_info.fn_strict && ufc_nullcount)
7302 if (ufl_info.fn_nargs == 4) {
7305 for (i = 0; i < set_count; i++) {
7307 d[0] = Int32GetDatum(_pos[i][0]);
7308 d[1] = Int32GetDatum(_pos[i][1]);
7311 d[2] = Int32GetDatum(_pos[i][0]);
7312 d[3] = Int32GetDatum(_pos[i][1]);
7316 a = construct_array(d, 4, INT4OID,
sizeof(
int32),
true,
'i');
7317 ufc_info->args[2].value = PointerGetDatum(a);
7318 ufc_info->args[2].isnull =
FALSE;
7321 datum = FunctionCallInvoke(ufc_info);
7324 if (!ufc_info->isnull)
7327 pixel = DatumGetFloat8(datum);
7336 if (calltype == TEXTOID) {
7337 for (k = 0; k < spi_count; k++) SPI_freeplan(spi_plan[k]);
7341 for (k = 0; k < set_count; k++) {
7342 if (_rast[k] != NULL)
7344 if (pgrastpos[k] != -1)
7345 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
7349 elog(ERROR,
"RASTER_mapAlgebra2: Could not set pixel value of output raster");
7361 if (calltype == TEXTOID) {
7362 for (i = 0; i < spi_count; i++) {
7363 if (spi_plan[i] != NULL) SPI_freeplan(spi_plan[i]);
7368 for (k = 0; k < set_count; k++) {
7369 if (_rast[k] != NULL)
7371 if (pgrastpos[k] != -1)
7372 PG_FREE_IF_COPY(pgrast[k], pgrastpos[k]);
7377 if (!pgrtn) PG_RETURN_NULL();
7381 SET_VARSIZE(pgrtn, pgrtn->
size);
7382 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_band rt_band_reclass_exact(rt_band srcband, rt_reclassmap map, uint32_t hasnodata, double nodataval)
Returns new band with values reclassified.
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_colormap_arg rtpg_colormap_arg_init(void)
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)
Datum RASTER_reclass_exact(PG_FUNCTION_ARGS)
static rtpg_nmapalgebra_arg rtpg_nmapalgebra_arg_init(void)
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
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)
static rtpg_clip_arg rtpg_clip_arg_init(void)
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
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::@13 method
struct rt_reclassexpr_t::rt_reclassrange src
struct rt_reclassexpr_t::rt_reclassrange dst
struct rt_classpair_t * pairs
rtpg_nmapalgebra_callback_arg callback
FunctionCallInfo ufc_info
FunctionCallInfoBaseData fcinfo
union rtpg_nmapalgebra_callback_arg::@24 ufc_info_data
rtpg_nmapalgebraexpr_callback_arg callback
rtpg_nmapalgebra_arg bandarg
struct rtpg_nmapalgebraexpr_callback_arg::@26 nodatanodata
struct rtpg_nmapalgebraexpr_callback_arg::@25 expr[3]
struct rtpg_nmapalgebraexpr_callback_arg::@27 kw
rtpg_union_band_arg bandarg
rtpg_union_type uniontype