Create a new empty raster with having the same georeference as the provided raster
If this new raster is empty (width = 0 OR height = 0) then there is nothing to compute and we return it right now
Check if the raster has the required band. Otherwise, return a raster without band
We set the initial value of the future band to nodata value. If nodata value is null, then the raster will be initialized to rt_band_get_min_value but all the values should be recomputed anyway
Optimization: If the raster is only filled with nodata values return right now a raster filled with the nodatavalueexpr TODO: Call rt_band_check_isnodata instead?
Create the raster receiving all the computed values. Initialize it to the new initial value
We compute a value only for the withdata value neighborhood since the nodata value has already been set by the first optimization
5504 int x,
y,
nband,
width,
height, ngbwidth, ngbheight, winwidth, winheight, u, v, nIndex, nNullItems;
5506 double newnodatavalue = 0.0;
5507 double newinitialvalue = 0.0;
5508 double newval = 0.0;
5513 FunctionCallInfoData cbdata;
5515 ArrayType * neighborDatum;
5516 char * strFromText = NULL;
5517 text * txtNodataMode = NULL;
5518 text * txtCallbackParam = NULL;
5520 float fltReplace = 0;
5521 bool valuereplace =
false, pixelreplace, nNodataOnly =
true, nNullSkip =
false;
5522 Datum * neighborData = NULL;
5523 bool * neighborNulls = NULL;
5524 int neighborDims[2];
5533 if (PG_ARGISNULL(0)) {
5534 elog(WARNING,
"Raster is NULL. Returning NULL");
5540 pgraster = (
rt_pgraster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
5544 PG_FREE_IF_COPY(pgraster, 0);
5545 elog(ERROR,
"RASTER_mapAlgebraFctNgb: Could not deserialize raster");
5553 if (PG_ARGISNULL(1))
5556 nband = PG_GETARG_INT32(1);
5561 POSTGIS_RT_DEBUG(3,
"RASTER_mapAlgebraFctNgb: Creating new empty raster...");
5572 if ( NULL == newrast ) {
5574 PG_FREE_IF_COPY(pgraster, 0);
5575 elog(ERROR,
"RASTER_mapAlgebraFctNgb: Could not create a new raster");
5600 elog(NOTICE,
"Raster is empty. Returning an empty raster");
5602 PG_FREE_IF_COPY(pgraster, 0);
5606 if (NULL == pgrtn) {
5607 elog(ERROR,
"RASTER_mapAlgebraFctNgb: Could not serialize raster");
5611 SET_VARSIZE(pgrtn, pgrtn->
size);
5612 PG_RETURN_POINTER(pgrtn);
5615 POSTGIS_RT_DEBUGF(3,
"RASTER_mapAlgebraFctNgb: Getting raster band %d...", nband);
5622 elog(NOTICE,
"Raster does not have the required band. Returning a raster " 5625 PG_FREE_IF_COPY(pgraster, 0);
5629 if (NULL == pgrtn) {
5630 elog(ERROR,
"RASTER_mapAlgebraFctNgb: Could not serialize raster");
5634 SET_VARSIZE(pgrtn, pgrtn->
size);
5635 PG_RETURN_POINTER(pgrtn);
5640 if ( NULL == band ) {
5641 elog(NOTICE,
"Could not get the required band. Returning a raster " 5644 PG_FREE_IF_COPY(pgraster, 0);
5648 if (NULL == pgrtn) {
5649 elog(ERROR,
"RASTER_mapAlgebraFctNgb: Could not serialize raster");
5653 SET_VARSIZE(pgrtn, pgrtn->
size);
5654 PG_RETURN_POINTER(pgrtn);
5660 POSTGIS_RT_DEBUG(3,
"RASTER_mapAlgebraFctNgb: Getting NODATA value for band...");
5677 newinitialvalue = newnodatavalue;
5684 if (PG_ARGISNULL(2)) {
5690 POSTGIS_RT_DEBUGF(3,
"RASTER_mapAlgebraFctNgb: Pixeltype parameter: %s", strFromText);
5693 if (newpixeltype ==
PT_END)
5697 if (newpixeltype ==
PT_END) {
5700 PG_FREE_IF_COPY(pgraster, 0);
5703 elog(ERROR,
"RASTER_mapAlgebraFctNgb: Invalid pixeltype");
5711 if (PG_ARGISNULL(5)) {
5714 PG_FREE_IF_COPY(pgraster, 0);
5717 elog(ERROR,
"RASTER_mapAlgebraFctNgb: Required function is missing");
5721 oid = PG_GETARG_OID(5);
5722 if (oid == InvalidOid) {
5725 PG_FREE_IF_COPY(pgraster, 0);
5728 elog(ERROR,
"RASTER_mapAlgebraFctNgb: Got invalid function object id");
5732 fmgr_info(oid, &cbinfo);
5735 if (cbinfo.fn_retset) {
5738 PG_FREE_IF_COPY(pgraster, 0);
5741 elog(ERROR,
"RASTER_mapAlgebraFctNgb: Function provided must return double precision not resultset");
5745 else if (cbinfo.fn_nargs != 3) {
5748 PG_FREE_IF_COPY(pgraster, 0);
5751 elog(ERROR,
"RASTER_mapAlgebraFctNgb: Function does not have three input parameters");
5755 if (func_volatile(oid) ==
'v') {
5756 elog(NOTICE,
"Function provided is VOLATILE. Unless required and for best performance, function should be IMMUTABLE or STABLE");
5760 #if POSTGIS_PGSQL_VERSION <= 90 5761 InitFunctionCallInfoData(cbdata, &cbinfo, 3, InvalidOid, NULL);
5763 InitFunctionCallInfoData(cbdata, &cbinfo, 3, InvalidOid, NULL, NULL);
5765 memset(cbdata.argnull,
FALSE,
sizeof(
bool) * 3);
5768 if (PG_ARGISNULL(7)) {
5769 if (cbinfo.fn_strict) {
5772 PG_FREE_IF_COPY(pgraster, 0);
5775 elog(ERROR,
"RASTER_mapAlgebraFctNgb: Strict callback functions cannot have NULL parameters");
5779 cbdata.arg[2] = (Datum)NULL;
5780 cbdata.argnull[2] =
TRUE;
5783 cbdata.arg[2] = PG_GETARG_DATUM(7);
5793 POSTGIS_RT_DEBUG(3,
"RASTER_mapAlgebraFctNgb: Band is a nodata band, returning " 5794 "a raster filled with nodata");
5797 newinitialvalue,
TRUE, newnodatavalue, 0);
5800 PG_FREE_IF_COPY(pgraster, 0);
5805 if (NULL == pgrtn) {
5806 elog(ERROR,
"RASTER_mapAlgebraFctNgb: Could not serialize raster");
5810 SET_VARSIZE(pgrtn, pgrtn->
size);
5811 PG_RETURN_POINTER(pgrtn);
5820 newinitialvalue,
TRUE, newnodatavalue, 0);
5824 if ( NULL == newband ) {
5825 elog(NOTICE,
"Could not modify band for new raster. Returning new " 5826 "raster with the original band");
5829 PG_FREE_IF_COPY(pgraster, 0);
5834 if (NULL == pgrtn) {
5835 elog(ERROR,
"RASTER_mapAlgebraFctNgb: Could not serialize raster");
5839 SET_VARSIZE(pgrtn, pgrtn->
size);
5840 PG_RETURN_POINTER(pgrtn);
5844 if (PG_ARGISNULL(3) || PG_GETARG_INT32(3) <= 0) {
5845 elog(NOTICE,
"Neighborhood width is NULL or <= 0. Returning new " 5846 "raster with the original band");
5849 PG_FREE_IF_COPY(pgraster, 0);
5854 if (NULL == pgrtn) {
5855 elog(ERROR,
"RASTER_mapAlgebraFctNgb: Could not serialize raster");
5859 SET_VARSIZE(pgrtn, pgrtn->
size);
5860 PG_RETURN_POINTER(pgrtn);
5863 ngbwidth = PG_GETARG_INT32(3);
5864 winwidth = ngbwidth * 2 + 1;
5867 if (PG_ARGISNULL(4) || PG_GETARG_INT32(4) <= 0) {
5868 elog(NOTICE,
"Neighborhood height is NULL or <= 0. Returning new " 5869 "raster with the original band");
5872 PG_FREE_IF_COPY(pgraster, 0);
5877 if (NULL == pgrtn) {
5878 elog(ERROR,
"RASTER_mapAlgebraFctNgb: Could not serialize raster");
5882 SET_VARSIZE(pgrtn, pgrtn->
size);
5883 PG_RETURN_POINTER(pgrtn);
5886 ngbheight = PG_GETARG_INT32(4);
5887 winheight = ngbheight * 2 + 1;
5890 if (PG_ARGISNULL(6)) {
5891 elog(NOTICE,
"Neighborhood NODATA behavior defaulting to 'ignore'");
5892 txtNodataMode = cstring_to_text(
"ignore");
5895 txtNodataMode = PG_GETARG_TEXT_P(6);
5898 txtCallbackParam = (text*)palloc(VARSIZE(txtNodataMode));
5899 SET_VARSIZE(txtCallbackParam, VARSIZE(txtNodataMode));
5900 memcpy((
void *)VARDATA(txtCallbackParam), (
void *)VARDATA(txtNodataMode), VARSIZE(txtNodataMode) - VARHDRSZ);
5903 cbdata.arg[1] = CStringGetDatum(txtCallbackParam);
5908 if (strcmp(strFromText,
"VALUE") == 0)
5909 valuereplace =
true;
5910 else if (strcmp(strFromText,
"IGNORE") != 0 && strcmp(strFromText,
"NULL") != 0) {
5912 if (sscanf(strFromText,
"%d", &intReplace) <= 0 && sscanf(strFromText,
"%f", &fltReplace) <= 0) {
5914 elog(NOTICE,
"Neighborhood NODATA mode is not recognized. Must be one of 'value', 'ignore', " 5915 "'NULL', or a numeric value. Returning new raster with the original band");
5918 pfree(txtCallbackParam);
5922 PG_FREE_IF_COPY(pgraster, 0);
5927 if (NULL == pgrtn) {
5928 elog(ERROR,
"RASTER_mapAlgebraFctNgb: Could not serialize raster");
5932 SET_VARSIZE(pgrtn, pgrtn->
size);
5933 PG_RETURN_POINTER(pgrtn);
5936 else if (strcmp(strFromText,
"NULL") == 0) {
5945 neighborData = (Datum *)palloc(winwidth * winheight *
sizeof(Datum));
5946 neighborNulls = (
bool *)palloc(winwidth * winheight *
sizeof(
bool));
5949 neighborDims[0] = winwidth;
5950 neighborDims[1] = winheight;
5957 get_typlenbyvalalign(FLOAT8OID, &typlen, &typbyval, &typalign);
5959 for (x = 0 + ngbwidth; x < width - ngbwidth; x++) {
5960 for(y = 0 + ngbheight; y < height - ngbheight; y++) {
5965 pixelreplace =
false;
5969 pixelreplace =
true;
5972 for (u = x - ngbwidth; u <= x + ngbwidth; u++) {
5973 for (v = y - ngbheight; v <= y + ngbheight; v++) {
5976 if (
FLT_NEQ(r, newnodatavalue)) {
5978 neighborData[nIndex] = Float8GetDatum((
double)r);
5979 neighborNulls[nIndex] =
false;
5980 nNodataOnly =
false;
5984 if (valuereplace && pixelreplace) {
5986 neighborData[nIndex] = Float8GetDatum((
double)rpix);
5987 neighborNulls[nIndex] =
false;
5992 neighborData[nIndex] = PointerGetDatum(NULL);
5993 neighborNulls[nIndex] =
true;
6000 neighborData[nIndex] = PointerGetDatum(NULL);
6001 neighborNulls[nIndex] =
true;
6013 if (!(nNodataOnly ||
6014 (nNullSkip && nNullItems > 0) ||
6015 (valuereplace && nNullItems > 0))) {
6017 x, y, winwidth, winheight);
6019 neighborDatum = construct_md_array((
void *)neighborData, neighborNulls, 2, neighborDims, neighborLbs,
6020 FLOAT8OID, typlen, typbyval, typalign);
6023 cbdata.arg[0] = PointerGetDatum(neighborDatum);
6026 tmpnewval = FunctionCallInvoke(&cbdata);
6029 if (cbdata.isnull) {
6030 newval = newnodatavalue;
6033 newval = DatumGetFloat8(tmpnewval);
6049 pfree(neighborNulls);
6050 pfree(neighborData);
6052 pfree(txtCallbackParam);
6055 PG_FREE_IF_COPY(pgraster, 0);
6059 POSTGIS_RT_DEBUG(3,
"RASTER_mapAlgebraFctNgb: raster modified, serializing it.");
6070 SET_VARSIZE(pgrtn, pgrtn->
size);
6071 PG_RETURN_POINTER(pgrtn);
char * text_to_cstring(const text *textptr)
void * rt_raster_serialize(rt_raster raster)
Return this raster in serialized form.
double rt_raster_get_x_offset(rt_raster raster)
Get raster x offset, in projection units.
double rt_raster_get_y_skew(rt_raster raster)
Get skew about the Y axis.
raster
Be careful!! Zeros function's input parameter can be a (height x width) array, not (width x height): ...
void rt_raster_set_skews(rt_raster raster, double skewX, double skewY)
Set skews about the X and Y axis.
int rt_raster_is_empty(rt_raster raster)
Return TRUE if the raster is empty.
#define POSTGIS_RT_DEBUGF(level, msg,...)
rt_errorstate rt_band_get_nodata(rt_band band, double *nodata)
Get NODATA value.
rt_errorstate rt_band_get_pixel(rt_band band, int x, int y, double *value, int *nodata)
Get pixel value.
void rt_raster_set_scale(rt_raster raster, double scaleX, double scaleY)
Set scale 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.
rt_pixtype rt_pixtype_index_from_name(const char *pixname)
void rt_raster_set_offsets(rt_raster raster, double x, double y)
Set insertion points in projection units.
rt_band rt_raster_get_band(rt_raster raster, int bandNum)
Return Nth band, or NULL if unavailable.
int rt_band_get_hasnodata_flag(rt_band band)
Get hasnodata flag value.
double rt_raster_get_x_scale(rt_raster raster)
Get scale X in projection units.
int rt_raster_has_band(rt_raster raster, int nband)
Return TRUE if the raster has a band of this number.
void rt_raster_set_srid(rt_raster raster, int32_t srid)
Set raster's SRID.
int32_t rt_raster_get_srid(rt_raster raster)
Get raster's SRID.
void rt_raster_destroy(rt_raster raster)
Release memory associated to a raster.
double rt_raster_get_y_scale(rt_raster raster)
Get scale Y in projection units.
rt_raster rt_raster_new(uint32_t width, uint32_t height)
Construct a raster with given dimensions.
double rt_raster_get_x_skew(rt_raster raster)
Get skew about the X axis.
uint16_t rt_raster_get_width(rt_raster raster)
double rt_band_get_min_value(rt_band band)
Returns the minimal possible value for the band according to the pixel type.
const char * rt_pixtype_name(rt_pixtype pixtype)
rt_pixtype rt_band_get_pixtype(rt_band band)
Return pixeltype of this band.
int rt_band_get_isnodata_flag(rt_band band)
Get isnodata flag value.
rt_errorstate rt_band_set_pixel(rt_band band, int x, int y, double val, int *converted)
Set single pixel's value.
uint16_t rt_raster_get_height(rt_raster raster)
#define POSTGIS_RT_DEBUG(level, msg)
rt_raster rt_raster_deserialize(void *serialized, int header_only)
Return a raster from a serialized form.
double rt_raster_get_y_offset(rt_raster raster)
Get raster y offset, in projection units.
char * rtpg_strtoupper(char *str)