PostGIS  2.4.9dev-r@@SVN_REVISION@@

◆ RASTER_nearestValue()

Datum RASTER_nearestValue ( PG_FUNCTION_ARGS  )

Definition at line 1839 of file rtpg_pixel.c.

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(), PG_FUNCTION_INFO_V1(), LWPOINT::point, POINTTYPE, rtrowdump::raster, RASTER_neighborhood(), 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, genraster::value, rt_pixel_t::value, pixval::x, POINT2D::x, pixval::y, and POINT2D::y.

Referenced by RASTER_pixelOfValue().

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