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

Definition at line 2047 of file rtpg_pixel.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.

2048 {
2049  rt_pgraster *pgraster = NULL;
2050  rt_raster raster = NULL;
2051  rt_band band = NULL;
2052  int bandindex = 1;
2053  int num_bands = 0;
2054  int x = 0;
2055  int y = 0;
2056  int _x = 0;
2057  int _y = 0;
2058  int distance[2] = {0};
2059  bool exclude_nodata_value = TRUE;
2060  double pixval;
2061  int isnodata = 0;
2062 
2063  rt_pixel npixels = NULL;
2064  int count;
2065  double **value2D = NULL;
2066  int **nodata2D = NULL;
2067 
2068  int i = 0;
2069  int j = 0;
2070  int k = 0;
2071  Datum *value1D = NULL;
2072  bool *nodata1D = NULL;
2073  int dim[2] = {0};
2074  int lbound[2] = {1, 1};
2075  ArrayType *mdArray = NULL;
2076 
2077  int16 typlen;
2078  bool typbyval;
2079  char typalign;
2080 
2081  /* pgraster is null, return nothing */
2082  if (PG_ARGISNULL(0))
2083  PG_RETURN_NULL();
2084  pgraster = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
2085 
2086  raster = rt_raster_deserialize(pgraster, FALSE);
2087  if (!raster) {
2088  PG_FREE_IF_COPY(pgraster, 0);
2089  elog(ERROR, "RASTER_neighborhood: Could not deserialize raster");
2090  PG_RETURN_NULL();
2091  }
2092 
2093  /* band index is 1-based */
2094  if (!PG_ARGISNULL(1))
2095  bandindex = PG_GETARG_INT32(1);
2096  num_bands = rt_raster_get_num_bands(raster);
2097  if (bandindex < 1 || bandindex > num_bands) {
2098  elog(NOTICE, "Invalid band index (must use 1-based). Returning NULL");
2099  rt_raster_destroy(raster);
2100  PG_FREE_IF_COPY(pgraster, 0);
2101  PG_RETURN_NULL();
2102  }
2103 
2104  /* pixel column, 1-based */
2105  x = PG_GETARG_INT32(2);
2106  _x = x - 1;
2107 
2108  /* pixel row, 1-based */
2109  y = PG_GETARG_INT32(3);
2110  _y = y - 1;
2111 
2112  /* distance X axis */
2113  distance[0] = PG_GETARG_INT32(4);
2114  if (distance[0] < 0) {
2115  elog(NOTICE, "Invalid value for distancex (must be >= zero). Returning NULL");
2116  rt_raster_destroy(raster);
2117  PG_FREE_IF_COPY(pgraster, 0);
2118  PG_RETURN_NULL();
2119  }
2120  distance[0] = (uint16_t) distance[0];
2121 
2122  /* distance Y axis */
2123  distance[1] = PG_GETARG_INT32(5);
2124  if (distance[1] < 0) {
2125  elog(NOTICE, "Invalid value for distancey (must be >= zero). Returning NULL");
2126  rt_raster_destroy(raster);
2127  PG_FREE_IF_COPY(pgraster, 0);
2128  PG_RETURN_NULL();
2129  }
2130  distance[1] = (uint16_t) distance[1];
2131 
2132  /* exclude_nodata_value flag */
2133  if (!PG_ARGISNULL(6))
2134  exclude_nodata_value = PG_GETARG_BOOL(6);
2135 
2136  /* get band */
2137  band = rt_raster_get_band(raster, bandindex - 1);
2138  if (!band) {
2139  elog(NOTICE, "Could not find band at index %d. Returning NULL", bandindex);
2140  rt_raster_destroy(raster);
2141  PG_FREE_IF_COPY(pgraster, 0);
2142  PG_RETURN_NULL();
2143  }
2144 
2145  /* get neighborhood */
2146  count = 0;
2147  npixels = NULL;
2148  if (distance[0] > 0 || distance[1] > 0) {
2149  count = rt_band_get_nearest_pixel(
2150  band,
2151  _x, _y,
2152  distance[0], distance[1],
2153  exclude_nodata_value,
2154  &npixels
2155  );
2156  /* error */
2157  if (count < 0) {
2158  elog(NOTICE, "Could not get the pixel's neighborhood for band at index %d", bandindex);
2159 
2160  rt_band_destroy(band);
2161  rt_raster_destroy(raster);
2162  PG_FREE_IF_COPY(pgraster, 0);
2163 
2164  PG_RETURN_NULL();
2165  }
2166  }
2167 
2168  /* get pixel's value */
2169  if (
2170  (_x >= 0 && _x < rt_band_get_width(band)) &&
2171  (_y >= 0 && _y < rt_band_get_height(band))
2172  ) {
2173  if (rt_band_get_pixel(
2174  band,
2175  _x, _y,
2176  &pixval,
2177  &isnodata
2178  ) != ES_NONE) {
2179  elog(NOTICE, "Could not get the pixel of band at index %d. Returning NULL", bandindex);
2180  rt_band_destroy(band);
2181  rt_raster_destroy(raster);
2182  PG_FREE_IF_COPY(pgraster, 0);
2183  PG_RETURN_NULL();
2184  }
2185  }
2186  /* outside band extent, set to NODATA */
2187  else {
2188  /* has NODATA, use NODATA */
2189  if (rt_band_get_hasnodata_flag(band))
2190  rt_band_get_nodata(band, &pixval);
2191  /* no NODATA, use min possible value */
2192  else
2193  pixval = rt_band_get_min_value(band);
2194  isnodata = 1;
2195  }
2196  POSTGIS_RT_DEBUGF(4, "pixval: %f", pixval);
2197 
2198 
2199  /* add pixel to neighborhood */
2200  count++;
2201  if (count > 1)
2202  npixels = (rt_pixel) repalloc(npixels, sizeof(struct rt_pixel_t) * count);
2203  else
2204  npixels = (rt_pixel) palloc(sizeof(struct rt_pixel_t));
2205  if (npixels == NULL) {
2206 
2207  rt_band_destroy(band);
2208  rt_raster_destroy(raster);
2209  PG_FREE_IF_COPY(pgraster, 0);
2210 
2211  elog(ERROR, "RASTER_neighborhood: Could not reallocate memory for neighborhood");
2212  PG_RETURN_NULL();
2213  }
2214  npixels[count - 1].x = _x;
2215  npixels[count - 1].y = _y;
2216  npixels[count - 1].nodata = 1;
2217  npixels[count - 1].value = pixval;
2218 
2219  /* set NODATA */
2220  if (!exclude_nodata_value || !isnodata) {
2221  npixels[count - 1].nodata = 0;
2222  }
2223 
2224  /* free unnecessary stuff */
2225  rt_band_destroy(band);
2226  rt_raster_destroy(raster);
2227  PG_FREE_IF_COPY(pgraster, 0);
2228 
2229  /* convert set of rt_pixel to 2D array */
2230  /* dim is passed with element 0 being Y-axis and element 1 being X-axis */
2231  count = rt_pixel_set_to_array(
2232  npixels, count, NULL,
2233  _x, _y,
2234  distance[0], distance[1],
2235  &value2D,
2236  &nodata2D,
2237  &(dim[1]), &(dim[0])
2238  );
2239  pfree(npixels);
2240  if (count != ES_NONE) {
2241  elog(NOTICE, "Could not create 2D array of neighborhood");
2242  PG_RETURN_NULL();
2243  }
2244 
2245  /* 1D arrays for values and nodata from 2D arrays */
2246  value1D = palloc(sizeof(Datum) * dim[0] * dim[1]);
2247  nodata1D = palloc(sizeof(bool) * dim[0] * dim[1]);
2248 
2249  if (value1D == NULL || nodata1D == NULL) {
2250 
2251  for (i = 0; i < dim[0]; i++) {
2252  pfree(value2D[i]);
2253  pfree(nodata2D[i]);
2254  }
2255  pfree(value2D);
2256  pfree(nodata2D);
2257 
2258  elog(ERROR, "RASTER_neighborhood: Could not allocate memory for return 2D array");
2259  PG_RETURN_NULL();
2260  }
2261 
2262  /* copy values from 2D arrays to 1D arrays */
2263  k = 0;
2264  /* Y-axis */
2265  for (i = 0; i < dim[0]; i++) {
2266  /* X-axis */
2267  for (j = 0; j < dim[1]; j++) {
2268  nodata1D[k] = (bool) nodata2D[i][j];
2269  if (!nodata1D[k])
2270  value1D[k] = Float8GetDatum(value2D[i][j]);
2271  else
2272  value1D[k] = PointerGetDatum(NULL);
2273 
2274  k++;
2275  }
2276  }
2277 
2278  /* no more need for 2D arrays */
2279  for (i = 0; i < dim[0]; i++) {
2280  pfree(value2D[i]);
2281  pfree(nodata2D[i]);
2282  }
2283  pfree(value2D);
2284  pfree(nodata2D);
2285 
2286  /* info about the type of item in the multi-dimensional array (float8). */
2287  get_typlenbyvalalign(FLOAT8OID, &typlen, &typbyval, &typalign);
2288 
2289  mdArray = construct_md_array(
2290  value1D, nodata1D,
2291  2, dim, lbound,
2292  FLOAT8OID,
2293  typlen, typbyval, typalign
2294  );
2295 
2296  pfree(value1D);
2297  pfree(nodata1D);
2298 
2299  PG_RETURN_ARRAYTYPE_P(mdArray);
2300 }
uint16_t rt_raster_get_num_bands(rt_raster raster)
Definition: rt_raster.c:372
struct rt_pixel_t * rt_pixel
Definition: librtcore.h:147
tuple band
Definition: ovdump.py:57
double value
Definition: librtcore.h:2289
tuple raster
Be careful!! Zeros function's input parameter can be a (height x width) array, not (width x height): ...
Definition: rtrowdump.py:121
void rt_band_destroy(rt_band band)
Destroy a raster band.
Definition: rt_band.c:242
tuple pixval
Definition: pixval.py:93
#define POSTGIS_RT_DEBUGF(level, msg,...)
Definition: rtpostgis.h:65
rt_errorstate rt_band_get_nodata(rt_band band, double *nodata)
Get NODATA value.
Definition: rt_band.c:1597
rt_errorstate rt_band_get_pixel(rt_band band, int x, int y, double *value, int *nodata)
Get pixel value.
Definition: rt_band.c:1088
rt_band rt_raster_get_band(rt_raster raster, int bandNum)
Return Nth band, or NULL if unavailable.
Definition: rt_raster.c:381
int rt_band_get_hasnodata_flag(rt_band band)
Get hasnodata flag value.
Definition: rt_band.c:541
int count
Definition: genraster.py:56
uint16_t rt_band_get_width(rt_band band)
Return width of this band.
Definition: rt_band.c:507
uint16_t rt_band_get_height(rt_band band)
Return height of this band.
Definition: rt_band.c:516
void rt_raster_destroy(rt_raster raster)
Release memory associated to a raster.
Definition: rt_raster.c:82
uint8_t nodata
Definition: librtcore.h:2288
tuple x
Definition: pixval.py:53
Datum distance(PG_FUNCTION_ARGS)
rt_errorstate rt_pixel_set_to_array(rt_pixel npixel, uint32_t count, rt_mask mask, int x, int y, uint16_t distancex, uint16_t distancey, double ***value, int ***nodata, int *dimx, int *dimy)
Definition: rt_pixel.c:286
#define FALSE
Definition: dbfopen.c:168
Struct definitions.
Definition: librtcore.h:2201
double rt_band_get_min_value(rt_band band)
Returns the minimal possible value for the band according to the pixel type.
Definition: rt_band.c:1612
uint32_t 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_band.c:1241
#define TRUE
Definition: dbfopen.c:169
rt_raster rt_raster_deserialize(void *serialized, int header_only)
Return a raster from a serialized form.
Definition: rt_serialize.c:717
tuple y
Definition: pixval.py:54

Here is the call graph for this function: