PostGIS  3.1.6dev-r@@SVN_REVISION@@

◆ RASTER_nearestValue()

Datum RASTER_nearestValue ( PG_FUNCTION_ARGS  )

Definition at line 1845 of file rtpg_pixel.c.

1846 {
1847  rt_pgraster *pgraster = NULL;
1848  rt_raster raster = NULL;
1849  rt_band band = NULL;
1850  int bandindex = 1;
1851  int num_bands = 0;
1852  GSERIALIZED *geom;
1853  bool exclude_nodata_value = TRUE;
1854  LWGEOM *lwgeom;
1855  LWPOINT *point = NULL;
1856  POINT2D p;
1857 
1858  double x;
1859  double y;
1860  int count;
1861  rt_pixel npixels = NULL;
1862  double value = 0;
1863  int hasvalue = 0;
1864  int isnodata = 0;
1865 
1866  if (PG_ARGISNULL(0))
1867  PG_RETURN_NULL();
1868  pgraster = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
1869  raster = rt_raster_deserialize(pgraster, FALSE);
1870  if (!raster) {
1871  PG_FREE_IF_COPY(pgraster, 0);
1872  elog(ERROR, "RASTER_nearestValue: Could not deserialize raster");
1873  PG_RETURN_NULL();
1874  }
1875 
1876  /* band index is 1-based */
1877  if (!PG_ARGISNULL(1))
1878  bandindex = PG_GETARG_INT32(1);
1879  num_bands = rt_raster_get_num_bands(raster);
1880  if (bandindex < 1 || bandindex > num_bands) {
1881  elog(NOTICE, "Invalid band index (must use 1-based). Returning NULL");
1883  PG_FREE_IF_COPY(pgraster, 0);
1884  PG_RETURN_NULL();
1885  }
1886 
1887  /* point */
1888  geom = PG_GETARG_GSERIALIZED_P(2);
1889  if (gserialized_get_type(geom) != POINTTYPE) {
1890  elog(NOTICE, "Geometry provided must be a point");
1892  PG_FREE_IF_COPY(pgraster, 0);
1893  PG_FREE_IF_COPY(geom, 2);
1894  PG_RETURN_NULL();
1895  }
1896 
1897  /* exclude_nodata_value flag */
1898  if (!PG_ARGISNULL(3))
1899  exclude_nodata_value = PG_GETARG_BOOL(3);
1900 
1901  /* SRIDs of raster and geometry must match */
1903  elog(NOTICE, "SRIDs of geometry and raster do not match");
1905  PG_FREE_IF_COPY(pgraster, 0);
1906  PG_FREE_IF_COPY(geom, 2);
1907  PG_RETURN_NULL();
1908  }
1909 
1910  /* get band */
1911  band = rt_raster_get_band(raster, bandindex - 1);
1912  if (!band) {
1913  elog(NOTICE, "Could not find band at index %d. Returning NULL", bandindex);
1915  PG_FREE_IF_COPY(pgraster, 0);
1916  PG_FREE_IF_COPY(geom, 2);
1917  PG_RETURN_NULL();
1918  }
1919 
1920  /* process geometry */
1921  lwgeom = lwgeom_from_gserialized(geom);
1922 
1923  if (lwgeom_is_empty(lwgeom)) {
1924  elog(NOTICE, "Geometry provided cannot be empty");
1926  PG_FREE_IF_COPY(pgraster, 0);
1927  PG_FREE_IF_COPY(geom, 2);
1928  PG_RETURN_NULL();
1929  }
1930 
1931  /* Get a 2D version of the geometry if necessary */
1932  if (lwgeom_ndims(lwgeom) > 2) {
1933  LWGEOM *lwgeom2d = lwgeom_force_2d(lwgeom);
1934  lwgeom_free(lwgeom);
1935  lwgeom = lwgeom2d;
1936  }
1937 
1938  point = lwgeom_as_lwpoint(lwgeom);
1939  getPoint2d_p(point->point, 0, &p);
1940 
1942  raster,
1943  p.x, p.y,
1944  &x, &y,
1945  NULL
1946  ) != ES_NONE) {
1948  PG_FREE_IF_COPY(pgraster, 0);
1949  lwgeom_free(lwgeom);
1950  PG_FREE_IF_COPY(geom, 2);
1951  elog(ERROR, "RASTER_nearestValue: Could not compute pixel coordinates from spatial coordinates");
1952  PG_RETURN_NULL();
1953  }
1954 
1955  /* get value at point */
1956  if (
1957  (x >= 0 && x < rt_raster_get_width(raster)) &&
1958  (y >= 0 && y < rt_raster_get_height(raster))
1959  ) {
1960  if (rt_band_get_pixel(band, x, y, &value, &isnodata) != ES_NONE) {
1962  PG_FREE_IF_COPY(pgraster, 0);
1963  lwgeom_free(lwgeom);
1964  PG_FREE_IF_COPY(geom, 2);
1965  elog(ERROR, "RASTER_nearestValue: Could not get pixel value for band at index %d", bandindex);
1966  PG_RETURN_NULL();
1967  }
1968 
1969  /* value at point, return value */
1970  if (!exclude_nodata_value || !isnodata) {
1972  PG_FREE_IF_COPY(pgraster, 0);
1973  lwgeom_free(lwgeom);
1974  PG_FREE_IF_COPY(geom, 2);
1975 
1976  PG_RETURN_FLOAT8(value);
1977  }
1978  }
1979 
1980  /* get neighborhood */
1982  band,
1983  x, y,
1984  0, 0,
1985  exclude_nodata_value,
1986  &npixels
1987  );
1989  /* error or no neighbors */
1990  if (count < 1) {
1991  /* error */
1992  if (count < 0)
1993  elog(NOTICE, "Could not get the nearest value for band at index %d", bandindex);
1994  /* no nearest pixel */
1995  else
1996  elog(NOTICE, "No nearest value found for band at index %d", bandindex);
1997 
1998  lwgeom_free(lwgeom);
1999  PG_FREE_IF_COPY(geom, 2);
2001  PG_FREE_IF_COPY(pgraster, 0);
2002  PG_RETURN_NULL();
2003  }
2004 
2005  /* more than one nearest value, see which one is closest */
2006  if (count > 1) {
2007  int i = 0;
2008  LWPOLY *poly = NULL;
2009  double lastdist = -1;
2010  double dist;
2011 
2012  for (i = 0; i < count; i++) {
2013  /* convex-hull of pixel */
2014  poly = rt_raster_pixel_as_polygon(raster, npixels[i].x, npixels[i].y);
2015  if (!poly) {
2016  lwgeom_free(lwgeom);
2017  PG_FREE_IF_COPY(geom, 2);
2019  PG_FREE_IF_COPY(pgraster, 0);
2020  elog(ERROR, "RASTER_nearestValue: Could not get polygon of neighboring pixel");
2021  PG_RETURN_NULL();
2022  }
2023 
2024  /* distance between convex-hull and point */
2025  dist = lwgeom_mindistance2d(lwpoly_as_lwgeom(poly), lwgeom);
2026  if (lastdist < 0 || dist < lastdist) {
2027  value = npixels[i].value;
2028  hasvalue = 1;
2029  }
2030  lastdist = dist;
2031 
2032  lwpoly_free(poly);
2033  }
2034  }
2035  else {
2036  value = npixels[0].value;
2037  hasvalue = 1;
2038  }
2039 
2040  pfree(npixels);
2041  lwgeom_free(lwgeom);
2042  PG_FREE_IF_COPY(geom, 2);
2044  PG_FREE_IF_COPY(pgraster, 0);
2045 
2046  if (hasvalue)
2047  PG_RETURN_FLOAT8(value);
2048  else
2049  PG_RETURN_NULL();
2050 }
#define TRUE
Definition: dbfopen.c:73
#define FALSE
Definition: dbfopen.c:72
int32_t gserialized_get_srid(const GSERIALIZED *g)
Extract the SRID from the serialized form (it is packed into three bytes so this is a handy function)...
Definition: gserialized.c:126
LWGEOM * lwgeom_from_gserialized(const GSERIALIZED *g)
Allocate a new LWGEOM from a GSERIALIZED.
Definition: gserialized.c:239
uint32_t gserialized_get_type(const GSERIALIZED *g)
Extract the geometry type from the serialized form (it hides in the anonymous data area,...
Definition: gserialized.c:89
int lwgeom_ndims(const LWGEOM *geom)
Return the number of dimensions (2, 3, 4) in a geometry.
Definition: lwgeom.c:938
void lwgeom_free(LWGEOM *geom)
Definition: lwgeom.c:1138
LWGEOM * lwpoly_as_lwgeom(const LWPOLY *obj)
Definition: lwgeom.c:312
int getPoint2d_p(const POINTARRAY *pa, uint32_t n, POINT2D *point)
Definition: lwgeom_api.c:343
#define POINTTYPE
LWTYPE numbers, used internally by PostGIS.
Definition: liblwgeom.h:116
double lwgeom_mindistance2d(const LWGEOM *lw1, const LWGEOM *lw2)
Function initializing min distance calculation.
Definition: measures.c:197
void lwpoly_free(LWPOLY *poly)
Definition: lwpoly.c:175
int32_t clamp_srid(int32_t srid)
Return a valid SRID from an arbitrary integer Raises a notice if what comes out is different from wha...
Definition: lwutil.c:333
LWGEOM * lwgeom_force_2d(const LWGEOM *geom)
Strip out the Z/M components of an LWGEOM.
Definition: lwgeom.c:776
int32_t rt_raster_get_srid(rt_raster raster)
Get raster's SRID.
Definition: rt_raster.c:356
rt_errorstate rt_raster_geopoint_to_cell(rt_raster raster, double xw, double yw, double *xr, double *yr, double *igt)
Convert an xw, yw map point to a xr, yr raster point.
Definition: rt_raster.c:804
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
void rt_raster_destroy(rt_raster raster)
Release memory associated to a raster.
Definition: rt_raster.c:82
LWPOLY * rt_raster_pixel_as_polygon(rt_raster raster, int x, int y)
Get a raster pixel as a polygon.
Definition: rt_geometry.c:610
@ ES_NONE
Definition: librtcore.h:180
void rt_band_destroy(rt_band band)
Destroy a raster band.
Definition: rt_band.c:340
uint16_t rt_raster_get_num_bands(rt_raster raster)
Definition: rt_raster.c:372
uint16_t rt_raster_get_height(rt_raster raster)
Definition: rt_raster.c:129
uint16_t rt_raster_get_width(rt_raster raster)
Definition: rt_raster.c:121
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
rt_raster rt_raster_deserialize(void *serialized, int header_only)
Return a raster from a serialized form.
Definition: rt_serialize.c:725
rt_band rt_raster_get_band(rt_raster raster, int bandNum)
Return Nth band, or NULL if unavailable.
Definition: rt_raster.c:381
static int lwgeom_is_empty(const LWGEOM *geom)
Return true or false depending on whether a geometry is an "empty" geometry (no vertices members)
Definition: lwinline.h:203
static LWPOINT * lwgeom_as_lwpoint(const LWGEOM *lwgeom)
Definition: lwinline.h:131
int value
Definition: genraster.py:62
int count
Definition: genraster.py:57
band
Definition: ovdump.py:58
raster
Be careful!! Zeros function's input parameter can be a (height x width) array, not (width x height): ...
Definition: rtrowdump.py:121
POINTARRAY * point
Definition: liblwgeom.h:485
double y
Definition: liblwgeom.h:404
double x
Definition: liblwgeom.h:404
double value
Definition: librtcore.h:2339
Struct definitions.
Definition: librtcore.h:2251

References ovdump::band, clamp_srid(), genraster::count, ES_NONE, FALSE, getPoint2d_p(), gserialized_get_srid(), gserialized_get_type(), lwgeom_as_lwpoint(), lwgeom_force_2d(), lwgeom_free(), lwgeom_from_gserialized(), lwgeom_is_empty(), lwgeom_mindistance2d(), lwgeom_ndims(), lwpoly_as_lwgeom(), lwpoly_free(), LWPOINT::point, POINTTYPE, rtrowdump::raster, rt_band_destroy(), rt_band_get_nearest_pixel(), rt_band_get_pixel(), rt_raster_deserialize(), rt_raster_destroy(), rt_raster_geopoint_to_cell(), rt_raster_get_band(), rt_raster_get_height(), rt_raster_get_num_bands(), rt_raster_get_srid(), rt_raster_get_width(), rt_raster_pixel_as_polygon(), TRUE, rt_pixel_t::value, genraster::value, POINT2D::x, pixval::x, POINT2D::y, and pixval::y.

Here is the call graph for this function: