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
5620 int x,
y,
nband, width, height, ngbwidth, ngbheight, winwidth, winheight, u, v, nIndex, nNullItems;
5622 double newnodatavalue = 0.0;
5623 double newinitialvalue = 0.0;
5624 double newval = 0.0;
5629 #if POSTGIS_PGSQL_VERSION < 120
5630 FunctionCallInfoData cbdata;
5632 LOCAL_FCINFO(cbdata, FUNC_MAX_ARGS);
5635 ArrayType * neighborDatum;
5636 char * strFromText = NULL;
5637 text * txtNodataMode = NULL;
5638 text * txtCallbackParam = NULL;
5640 float fltReplace = 0;
5641 bool valuereplace =
false, pixelreplace, nNodataOnly =
true, nNullSkip =
false;
5642 Datum * neighborData = NULL;
5643 bool * neighborNulls = NULL;
5644 int neighborDims[2];
5653 if (PG_ARGISNULL(0)) {
5654 elog(WARNING,
"Raster is NULL. Returning NULL");
5660 pgraster = (
rt_pgraster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
5664 PG_FREE_IF_COPY(pgraster, 0);
5665 elog(ERROR,
"RASTER_mapAlgebraFctNgb: Could not deserialize raster");
5673 if (PG_ARGISNULL(1))
5676 nband = PG_GETARG_INT32(1);
5681 POSTGIS_RT_DEBUG(3,
"RASTER_mapAlgebraFctNgb: Creating new empty raster...");
5692 if ( NULL == newrast ) {
5694 PG_FREE_IF_COPY(pgraster, 0);
5695 elog(ERROR,
"RASTER_mapAlgebraFctNgb: Could not create a new raster");
5720 elog(NOTICE,
"Raster is empty. Returning an empty raster");
5722 PG_FREE_IF_COPY(pgraster, 0);
5726 if (NULL == pgrtn) {
5727 elog(ERROR,
"RASTER_mapAlgebraFctNgb: Could not serialize raster");
5731 SET_VARSIZE(pgrtn, pgrtn->
size);
5732 PG_RETURN_POINTER(pgrtn);
5742 elog(NOTICE,
"Raster does not have the required band. Returning a raster "
5745 PG_FREE_IF_COPY(pgraster, 0);
5749 if (NULL == pgrtn) {
5750 elog(ERROR,
"RASTER_mapAlgebraFctNgb: Could not serialize raster");
5754 SET_VARSIZE(pgrtn, pgrtn->
size);
5755 PG_RETURN_POINTER(pgrtn);
5760 if ( NULL ==
band ) {
5761 elog(NOTICE,
"Could not get the required band. Returning a raster "
5764 PG_FREE_IF_COPY(pgraster, 0);
5768 if (NULL == pgrtn) {
5769 elog(ERROR,
"RASTER_mapAlgebraFctNgb: Could not serialize raster");
5773 SET_VARSIZE(pgrtn, pgrtn->
size);
5774 PG_RETURN_POINTER(pgrtn);
5780 POSTGIS_RT_DEBUG(3,
"RASTER_mapAlgebraFctNgb: Getting NODATA value for band...");
5797 newinitialvalue = newnodatavalue;
5804 if (PG_ARGISNULL(2)) {
5809 strFromText = text_to_cstring(PG_GETARG_TEXT_P(2));
5810 POSTGIS_RT_DEBUGF(3,
"RASTER_mapAlgebraFctNgb: Pixeltype parameter: %s", strFromText);
5813 if (newpixeltype ==
PT_END)
5817 if (newpixeltype ==
PT_END) {
5820 PG_FREE_IF_COPY(pgraster, 0);
5823 elog(ERROR,
"RASTER_mapAlgebraFctNgb: Invalid pixeltype");
5831 if (PG_ARGISNULL(5)) {
5834 PG_FREE_IF_COPY(pgraster, 0);
5837 elog(ERROR,
"RASTER_mapAlgebraFctNgb: Required function is missing");
5841 oid = PG_GETARG_OID(5);
5842 if (oid == InvalidOid) {
5845 PG_FREE_IF_COPY(pgraster, 0);
5848 elog(ERROR,
"RASTER_mapAlgebraFctNgb: Got invalid function object id");
5852 fmgr_info(oid, &cbinfo);
5855 if (cbinfo.fn_retset) {
5858 PG_FREE_IF_COPY(pgraster, 0);
5861 elog(ERROR,
"RASTER_mapAlgebraFctNgb: Function provided must return double precision not resultset");
5865 else if (cbinfo.fn_nargs != 3) {
5868 PG_FREE_IF_COPY(pgraster, 0);
5871 elog(ERROR,
"RASTER_mapAlgebraFctNgb: Function does not have three input parameters");
5875 if (func_volatile(oid) ==
'v') {
5876 elog(NOTICE,
"Function provided is VOLATILE. Unless required and for best performance, function should be IMMUTABLE or STABLE");
5880 #if POSTGIS_PGSQL_VERSION < 120
5881 InitFunctionCallInfoData(cbdata, &cbinfo, 3, InvalidOid, NULL, NULL);
5882 memset(cbdata.argnull,
FALSE,
sizeof(
bool) * 3);
5884 InitFunctionCallInfoData(*cbdata, &cbinfo, 3, InvalidOid, NULL, NULL);
5885 cbdata->args[0].isnull =
FALSE;
5886 cbdata->args[1].isnull =
FALSE;
5887 cbdata->args[2].isnull =
FALSE;
5891 if (PG_ARGISNULL(7)) {
5892 if (cbinfo.fn_strict) {
5895 PG_FREE_IF_COPY(pgraster, 0);
5898 elog(ERROR,
"RASTER_mapAlgebraFctNgb: Strict callback functions cannot have NULL parameters");
5902 #if POSTGIS_PGSQL_VERSION < 120
5903 cbdata.arg[2] = (Datum)NULL;
5904 cbdata.argnull[2] =
TRUE;
5906 cbdata->args[2].value = (Datum)NULL;
5907 cbdata->args[2].isnull =
TRUE;
5911 #if POSTGIS_PGSQL_VERSION < 120
5912 cbdata.arg[2] = PG_GETARG_DATUM(7);
5914 cbdata->args[2].value = PG_GETARG_DATUM(7);
5925 POSTGIS_RT_DEBUG(3,
"RASTER_mapAlgebraFctNgb: Band is a nodata band, returning "
5926 "a raster filled with nodata");
5929 newinitialvalue,
TRUE, newnodatavalue, 0);
5932 PG_FREE_IF_COPY(pgraster, 0);
5937 if (NULL == pgrtn) {
5938 elog(ERROR,
"RASTER_mapAlgebraFctNgb: Could not serialize raster");
5942 SET_VARSIZE(pgrtn, pgrtn->
size);
5943 PG_RETURN_POINTER(pgrtn);
5952 newinitialvalue,
TRUE, newnodatavalue, 0);
5956 if ( NULL == newband ) {
5957 elog(NOTICE,
"Could not modify band for new raster. Returning new "
5958 "raster with the original band");
5961 PG_FREE_IF_COPY(pgraster, 0);
5966 if (NULL == pgrtn) {
5967 elog(ERROR,
"RASTER_mapAlgebraFctNgb: Could not serialize raster");
5971 SET_VARSIZE(pgrtn, pgrtn->
size);
5972 PG_RETURN_POINTER(pgrtn);
5976 if (PG_ARGISNULL(3) || PG_GETARG_INT32(3) <= 0) {
5977 elog(NOTICE,
"Neighborhood width is NULL or <= 0. Returning new "
5978 "raster with the original band");
5981 PG_FREE_IF_COPY(pgraster, 0);
5986 if (NULL == pgrtn) {
5987 elog(ERROR,
"RASTER_mapAlgebraFctNgb: Could not serialize raster");
5991 SET_VARSIZE(pgrtn, pgrtn->
size);
5992 PG_RETURN_POINTER(pgrtn);
5995 ngbwidth = PG_GETARG_INT32(3);
5996 winwidth = ngbwidth * 2 + 1;
5999 if (PG_ARGISNULL(4) || PG_GETARG_INT32(4) <= 0) {
6000 elog(NOTICE,
"Neighborhood height is NULL or <= 0. Returning new "
6001 "raster with the original band");
6004 PG_FREE_IF_COPY(pgraster, 0);
6009 if (NULL == pgrtn) {
6010 elog(ERROR,
"RASTER_mapAlgebraFctNgb: Could not serialize raster");
6014 SET_VARSIZE(pgrtn, pgrtn->
size);
6015 PG_RETURN_POINTER(pgrtn);
6018 ngbheight = PG_GETARG_INT32(4);
6019 winheight = ngbheight * 2 + 1;
6022 if (PG_ARGISNULL(6)) {
6023 elog(NOTICE,
"Neighborhood NODATA behavior defaulting to 'ignore'");
6024 txtNodataMode = cstring_to_text(
"ignore");
6027 txtNodataMode = PG_GETARG_TEXT_P(6);
6030 txtCallbackParam = (text*)palloc(VARSIZE(txtNodataMode));
6031 SET_VARSIZE(txtCallbackParam, VARSIZE(txtNodataMode));
6032 memcpy((
void *)VARDATA(txtCallbackParam), (
void *)VARDATA(txtNodataMode), VARSIZE(txtNodataMode) - VARHDRSZ);
6035 #if POSTGIS_PGSQL_VERSION < 120
6036 cbdata.arg[1] = CStringGetDatum(txtCallbackParam);
6038 cbdata->args[1].value = CStringGetDatum(txtCallbackParam);
6041 strFromText = text_to_cstring(txtNodataMode);
6044 if (strcmp(strFromText,
"VALUE") == 0)
6045 valuereplace =
true;
6046 else if (strcmp(strFromText,
"IGNORE") != 0 && strcmp(strFromText,
"NULL") != 0) {
6048 if (sscanf(strFromText,
"%d", &intReplace) <= 0 && sscanf(strFromText,
"%f", &fltReplace) <= 0) {
6050 elog(NOTICE,
"Neighborhood NODATA mode is not recognized. Must be one of 'value', 'ignore', "
6051 "'NULL', or a numeric value. Returning new raster with the original band");
6054 pfree(txtCallbackParam);
6058 PG_FREE_IF_COPY(pgraster, 0);
6063 if (NULL == pgrtn) {
6064 elog(ERROR,
"RASTER_mapAlgebraFctNgb: Could not serialize raster");
6068 SET_VARSIZE(pgrtn, pgrtn->
size);
6069 PG_RETURN_POINTER(pgrtn);
6072 else if (strcmp(strFromText,
"NULL") == 0) {
6081 neighborData = (Datum *)palloc(winwidth * winheight *
sizeof(Datum));
6082 neighborNulls = (
bool *)palloc(winwidth * winheight *
sizeof(
bool));
6085 neighborDims[0] = winwidth;
6086 neighborDims[1] = winheight;
6093 get_typlenbyvalalign(FLOAT8OID, &typlen, &typbyval, &typalign);
6095 for (
x = 0 + ngbwidth;
x < width - ngbwidth;
x++) {
6096 for(
y = 0 + ngbheight;
y < height - ngbheight;
y++) {
6101 pixelreplace =
false;
6105 pixelreplace =
true;
6108 for (u =
x - ngbwidth; u <=
x + ngbwidth; u++) {
6109 for (v =
y - ngbheight; v <=
y + ngbheight; v++) {
6114 neighborData[nIndex] = Float8GetDatum((
double)
r);
6115 neighborNulls[nIndex] =
false;
6116 nNodataOnly =
false;
6120 if (valuereplace && pixelreplace) {
6122 neighborData[nIndex] = Float8GetDatum((
double)rpix);
6123 neighborNulls[nIndex] =
false;
6128 neighborData[nIndex] = PointerGetDatum(NULL);
6129 neighborNulls[nIndex] =
true;
6136 neighborData[nIndex] = PointerGetDatum(NULL);
6137 neighborNulls[nIndex] =
true;
6149 if (!(nNodataOnly ||
6150 (nNullSkip && nNullItems > 0) ||
6151 (valuereplace && nNullItems > 0))) {
6153 x,
y, winwidth, winheight);
6155 neighborDatum = construct_md_array((
void *)neighborData, neighborNulls, 2, neighborDims, neighborLbs,
6156 FLOAT8OID, typlen, typbyval, typalign);
6158 #if POSTGIS_PGSQL_VERSION < 120
6160 cbdata.arg[0] = PointerGetDatum(neighborDatum);
6163 tmpnewval = FunctionCallInvoke(&cbdata);
6166 if (cbdata.isnull) {
6167 newval = newnodatavalue;
6171 cbdata->args[0].value = PointerGetDatum(neighborDatum);
6174 tmpnewval = FunctionCallInvoke(cbdata);
6179 newval = newnodatavalue;
6183 newval = DatumGetFloat8(tmpnewval);
6199 pfree(neighborNulls);
6200 pfree(neighborData);
6202 pfree(txtCallbackParam);
6205 PG_FREE_IF_COPY(pgraster, 0);
6209 POSTGIS_RT_DEBUG(3,
"RASTER_mapAlgebraFctNgb: raster modified, serializing it.");
6220 SET_VARSIZE(pgrtn, pgrtn->
size);
6221 PG_RETURN_POINTER(pgrtn);
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_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.
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.
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.
uint16_t rt_raster_get_height(rt_raster raster)
void rt_raster_set_srid(rt_raster raster, int32_t srid)
Set raster's SRID.
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)
double rt_raster_get_y_scale(rt_raster raster)
Get scale Y in projection units.
double rt_raster_get_y_skew(rt_raster raster)
Get skew about the Y axis.
void rt_raster_set_offsets(rt_raster raster, double x, double y)
Set insertion points in projection units.
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.
int rt_raster_is_empty(rt_raster raster)
Return TRUE if the raster is empty.
rt_band rt_raster_get_band(rt_raster raster, int bandNum)
Return Nth band, or NULL if unavailable.
raster
Be careful!! Zeros function's input parameter can be a (height x width) array, not (width x height): ...
char * rtpg_strtoupper(char *str)
#define POSTGIS_RT_DEBUG(level, msg)
#define POSTGIS_RT_DEBUGF(level, msg,...)