PostGIS  2.5.1dev-r@@SVN_REVISION@@

◆ RASTER_neighborhood()

Datum RASTER_neighborhood ( PG_FUNCTION_ARGS  )

Definition at line 2050 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.

Referenced by RASTER_nearestValue().

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