PostGIS  2.1.10dev-r@@SVN_REVISION@@
Datum RASTER_neighborhood ( PG_FUNCTION_ARGS  )

Definition at line 5069 of file rt_pg.c.

References ovdump::band, genraster::count, distance(), ES_NONE, FALSE, rt_pixel_t::nodata, pixval::pixval, POSTGIS_RT_DEBUGF, rtrowdump::raster, rt_band_destroy(), rt_band_get_hasnodata_flag(), rt_band_get_height(), rt_band_get_min_value(), rt_band_get_nearest_pixel(), rt_band_get_nodata(), rt_band_get_pixel(), rt_band_get_width(), rt_pixel_set_to_array(), rt_raster_deserialize(), rt_raster_destroy(), rt_raster_get_band(), rt_raster_get_num_bands(), TRUE, rt_pixel_t::value, pixval::x, rt_pixel_t::x, pixval::y, and rt_pixel_t::y.

5070 {
5071  rt_pgraster *pgraster = NULL;
5072  rt_raster raster = NULL;
5073  rt_band band = NULL;
5074  int bandindex = 1;
5075  int num_bands = 0;
5076  int x = 0;
5077  int y = 0;
5078  int _x = 0;
5079  int _y = 0;
5080  int distance[2] = {0};
5081  bool exclude_nodata_value = TRUE;
5082  double pixval;
5083  int isnodata = 0;
5084 
5085  rt_pixel npixels = NULL;
5086  int count;
5087  double **value2D = NULL;
5088  int **nodata2D = NULL;
5089 
5090  int i = 0;
5091  int j = 0;
5092  int k = 0;
5093  Datum *value1D = NULL;
5094  bool *nodata1D = NULL;
5095  int dim[2] = {0};
5096  int lbound[2] = {1, 1};
5097  ArrayType *mdArray = NULL;
5098 
5099  int16 typlen;
5100  bool typbyval;
5101  char typalign;
5102 
5103  /* pgraster is null, return nothing */
5104  if (PG_ARGISNULL(0))
5105  PG_RETURN_NULL();
5106  pgraster = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
5107 
5108  raster = rt_raster_deserialize(pgraster, FALSE);
5109  if (!raster) {
5110  PG_FREE_IF_COPY(pgraster, 0);
5111  elog(ERROR, "RASTER_neighborhood: Could not deserialize raster");
5112  PG_RETURN_NULL();
5113  }
5114 
5115  /* band index is 1-based */
5116  if (!PG_ARGISNULL(1))
5117  bandindex = PG_GETARG_INT32(1);
5118  num_bands = rt_raster_get_num_bands(raster);
5119  if (bandindex < 1 || bandindex > num_bands) {
5120  elog(NOTICE, "Invalid band index (must use 1-based). Returning NULL");
5121  rt_raster_destroy(raster);
5122  PG_FREE_IF_COPY(pgraster, 0);
5123  PG_RETURN_NULL();
5124  }
5125 
5126  /* pixel column, 1-based */
5127  x = PG_GETARG_INT32(2);
5128  _x = x - 1;
5129 
5130  /* pixel row, 1-based */
5131  y = PG_GETARG_INT32(3);
5132  _y = y - 1;
5133 
5134  /* distance X axis */
5135  distance[0] = PG_GETARG_INT32(4);
5136  if (distance[0] < 0) {
5137  elog(NOTICE, "Invalid value for distancex (must be >= zero). Returning NULL");
5138  rt_raster_destroy(raster);
5139  PG_FREE_IF_COPY(pgraster, 0);
5140  PG_RETURN_NULL();
5141  }
5142  distance[0] = (uint16_t) distance[0];
5143 
5144  /* distance Y axis */
5145  distance[1] = PG_GETARG_INT32(5);
5146  if (distance[1] < 0) {
5147  elog(NOTICE, "Invalid value for distancey (must be >= zero). Returning NULL");
5148  rt_raster_destroy(raster);
5149  PG_FREE_IF_COPY(pgraster, 0);
5150  PG_RETURN_NULL();
5151  }
5152  distance[1] = (uint16_t) distance[1];
5153 
5154  /* exclude_nodata_value flag */
5155  if (!PG_ARGISNULL(6))
5156  exclude_nodata_value = PG_GETARG_BOOL(6);
5157 
5158  /* get band */
5159  band = rt_raster_get_band(raster, bandindex - 1);
5160  if (!band) {
5161  elog(NOTICE, "Could not find band at index %d. Returning NULL", bandindex);
5162  rt_raster_destroy(raster);
5163  PG_FREE_IF_COPY(pgraster, 0);
5164  PG_RETURN_NULL();
5165  }
5166 
5167  /* get neighborhood */
5168  count = 0;
5169  npixels = NULL;
5170  if (distance[0] > 0 || distance[1] > 0) {
5171  count = rt_band_get_nearest_pixel(
5172  band,
5173  _x, _y,
5174  distance[0], distance[1],
5175  exclude_nodata_value,
5176  &npixels
5177  );
5178  /* error */
5179  if (count < 0) {
5180  elog(NOTICE, "Could not get the pixel's neighborhood for band at index %d", bandindex);
5181 
5182  rt_band_destroy(band);
5183  rt_raster_destroy(raster);
5184  PG_FREE_IF_COPY(pgraster, 0);
5185 
5186  PG_RETURN_NULL();
5187  }
5188  }
5189 
5190  /* get pixel's value */
5191  if (
5192  (_x >= 0 && _x < rt_band_get_width(band)) &&
5193  (_y >= 0 && _y < rt_band_get_height(band))
5194  ) {
5195  if (rt_band_get_pixel(
5196  band,
5197  _x, _y,
5198  &pixval,
5199  &isnodata
5200  ) != ES_NONE) {
5201  elog(NOTICE, "Could not get the pixel of band at index %d. Returning NULL", bandindex);
5202  rt_band_destroy(band);
5203  rt_raster_destroy(raster);
5204  PG_FREE_IF_COPY(pgraster, 0);
5205  PG_RETURN_NULL();
5206  }
5207  }
5208  /* outside band extent, set to NODATA */
5209  else {
5210  /* has NODATA, use NODATA */
5211  if (rt_band_get_hasnodata_flag(band))
5212  rt_band_get_nodata(band, &pixval);
5213  /* no NODATA, use min possible value */
5214  else
5215  pixval = rt_band_get_min_value(band);
5216  isnodata = 1;
5217  }
5218  POSTGIS_RT_DEBUGF(4, "pixval: %f", pixval);
5219 
5220 
5221  /* add pixel to neighborhood */
5222  count++;
5223  if (count > 1)
5224  npixels = (rt_pixel) repalloc(npixels, sizeof(struct rt_pixel_t) * count);
5225  else
5226  npixels = (rt_pixel) palloc(sizeof(struct rt_pixel_t));
5227  if (npixels == NULL) {
5228 
5229  rt_band_destroy(band);
5230  rt_raster_destroy(raster);
5231  PG_FREE_IF_COPY(pgraster, 0);
5232 
5233  elog(ERROR, "RASTER_neighborhood: Could not reallocate memory for neighborhood");
5234  PG_RETURN_NULL();
5235  }
5236  npixels[count - 1].x = _x;
5237  npixels[count - 1].y = _y;
5238  npixels[count - 1].nodata = 1;
5239  npixels[count - 1].value = pixval;
5240 
5241  /* set NODATA */
5242  if (!exclude_nodata_value || !isnodata) {
5243  npixels[count - 1].nodata = 0;
5244  }
5245 
5246  /* free unnecessary stuff */
5247  rt_band_destroy(band);
5248  rt_raster_destroy(raster);
5249  PG_FREE_IF_COPY(pgraster, 0);
5250 
5251  /* convert set of rt_pixel to 2D array */
5252  /* dim is passed with element 0 being Y-axis and element 1 being X-axis */
5253  count = rt_pixel_set_to_array(
5254  npixels, count,
5255  _x, _y,
5256  distance[0], distance[1],
5257  &value2D,
5258  &nodata2D,
5259  &(dim[1]), &(dim[0])
5260  );
5261  pfree(npixels);
5262  if (count != ES_NONE) {
5263  elog(NOTICE, "Could not create 2D array of neighborhood");
5264  PG_RETURN_NULL();
5265  }
5266 
5267  /* 1D arrays for values and nodata from 2D arrays */
5268  value1D = palloc(sizeof(Datum) * dim[0] * dim[1]);
5269  nodata1D = palloc(sizeof(bool) * dim[0] * dim[1]);
5270 
5271  if (value1D == NULL || nodata1D == NULL) {
5272 
5273  for (i = 0; i < dim[0]; i++) {
5274  pfree(value2D[i]);
5275  pfree(nodata2D[i]);
5276  }
5277  pfree(value2D);
5278  pfree(nodata2D);
5279 
5280  elog(ERROR, "RASTER_neighborhood: Could not allocate memory for return 2D array");
5281  PG_RETURN_NULL();
5282  }
5283 
5284  /* copy values from 2D arrays to 1D arrays */
5285  k = 0;
5286  /* Y-axis */
5287  for (i = 0; i < dim[0]; i++) {
5288  /* X-axis */
5289  for (j = 0; j < dim[1]; j++) {
5290  nodata1D[k] = (bool) nodata2D[i][j];
5291  if (!nodata1D[k])
5292  value1D[k] = Float8GetDatum(value2D[i][j]);
5293  else
5294  value1D[k] = PointerGetDatum(NULL);
5295 
5296  k++;
5297  }
5298  }
5299 
5300  /* no more need for 2D arrays */
5301  for (i = 0; i < dim[0]; i++) {
5302  pfree(value2D[i]);
5303  pfree(nodata2D[i]);
5304  }
5305  pfree(value2D);
5306  pfree(nodata2D);
5307 
5308  /* info about the type of item in the multi-dimensional array (float8). */
5309  get_typlenbyvalalign(FLOAT8OID, &typlen, &typbyval, &typalign);
5310 
5311  mdArray = construct_md_array(
5312  value1D, nodata1D,
5313  2, dim, lbound,
5314  FLOAT8OID,
5315  typlen, typbyval, typalign
5316  );
5317 
5318  pfree(value1D);
5319  pfree(nodata1D);
5320 
5321  PG_RETURN_ARRAYTYPE_P(mdArray);
5322 }
rt_errorstate rt_pixel_set_to_array(rt_pixel npixel, int count, int x, int y, uint16_t distancex, uint16_t distancey, double ***value, int ***nodata, int *dimx, int *dimy)
Definition: rt_api.c:1341
int rt_raster_get_num_bands(rt_raster raster)
Definition: rt_api.c:5677
void rt_raster_destroy(rt_raster raster)
Release memory associated to a raster.
Definition: rt_api.c:5387
int rt_band_get_nearest_pixel(rt_band band, int x, int y, uint16_t distancex, uint16_t distancey, int exclude_nodata_value, rt_pixel *npixels)
Get nearest pixel(s) with value (not NODATA) to specified pixel.
Definition: rt_api.c:2702
tuple band
Definition: ovdump.py:57
double value
Definition: rt_api.h:2263
tuple raster
Be careful!! Zeros function's input parameter can be a (height x width) array, not (width x height): ...
Definition: rtrowdump.py:123
struct rt_pixel_t * rt_pixel
Definition: rt_api.h:135
tuple pixval
Definition: pixval.py:93
double rt_band_get_min_value(rt_band band)
Returns the minimal possible value for the band according to the pixel type.
Definition: rt_api.c:3073
rt_errorstate rt_band_get_nodata(rt_band band, double *nodata)
Get NODATA value.
Definition: rt_api.c:3058
rt_errorstate rt_band_get_pixel(rt_band band, int x, int y, double *value, int *nodata)
Get pixel value.
Definition: rt_api.c:2549
int count
Definition: genraster.py:57
uint16_t rt_band_get_width(rt_band band)
Return width of this band.
Definition: rt_api.c:1909
uint8_t nodata
Definition: rt_api.h:2262
tuple x
Definition: pixval.py:53
Datum distance(PG_FUNCTION_ARGS)
void rt_band_destroy(rt_band band)
Destroy a raster band.
Definition: rt_api.c:1650
rt_band rt_raster_get_band(rt_raster raster, int n)
Return Nth band, or NULL if unavailable.
Definition: rt_api.c:5686
#define FALSE
Definition: dbfopen.c:169
Struct definitions.
Definition: rt_api.h:2175
#define POSTGIS_RT_DEBUGF(level, msg,...)
Definition: rt_pg.h:62
uint16_t rt_band_get_height(rt_band band)
Return height of this band.
Definition: rt_api.c:1918
rt_raster rt_raster_deserialize(void *serialized, int header_only)
Return a raster from a serialized form.
Definition: rt_api.c:8350
#define TRUE
Definition: dbfopen.c:170
int rt_band_get_hasnodata_flag(rt_band band)
Get hasnodata flag value.
Definition: rt_api.c:2002
tuple y
Definition: pixval.py:54

Here is the call graph for this function: