PostGIS  2.5.0dev-r@@SVN_REVISION@@

◆ RASTER_nearestValue()

Datum RASTER_nearestValue ( PG_FUNCTION_ARGS  )

Definition at line 1836 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().

1837 {
1838  rt_pgraster *pgraster = NULL;
1839  rt_raster raster = NULL;
1840  rt_band band = NULL;
1841  int bandindex = 1;
1842  int num_bands = 0;
1843  GSERIALIZED *geom;
1844  bool exclude_nodata_value = TRUE;
1845  LWGEOM *lwgeom;
1846  LWPOINT *point = NULL;
1847  POINT2D p;
1848 
1849  double x;
1850  double y;
1851  int count;
1852  rt_pixel npixels = NULL;
1853  double value = 0;
1854  int hasvalue = 0;
1855  int isnodata = 0;
1856 
1857  if (PG_ARGISNULL(0))
1858  PG_RETURN_NULL();
1859  pgraster = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
1860  raster = rt_raster_deserialize(pgraster, FALSE);
1861  if (!raster) {
1862  PG_FREE_IF_COPY(pgraster, 0);
1863  elog(ERROR, "RASTER_nearestValue: Could not deserialize raster");
1864  PG_RETURN_NULL();
1865  }
1866 
1867  /* band index is 1-based */
1868  if (!PG_ARGISNULL(1))
1869  bandindex = PG_GETARG_INT32(1);
1870  num_bands = rt_raster_get_num_bands(raster);
1871  if (bandindex < 1 || bandindex > num_bands) {
1872  elog(NOTICE, "Invalid band index (must use 1-based). Returning NULL");
1873  rt_raster_destroy(raster);
1874  PG_FREE_IF_COPY(pgraster, 0);
1875  PG_RETURN_NULL();
1876  }
1877 
1878  /* point */
1879  geom = PG_GETARG_GSERIALIZED_P(2);
1880  if (gserialized_get_type(geom) != POINTTYPE) {
1881  elog(NOTICE, "Geometry provided must be a point");
1882  rt_raster_destroy(raster);
1883  PG_FREE_IF_COPY(pgraster, 0);
1884  PG_FREE_IF_COPY(geom, 2);
1885  PG_RETURN_NULL();
1886  }
1887 
1888  /* exclude_nodata_value flag */
1889  if (!PG_ARGISNULL(3))
1890  exclude_nodata_value = PG_GETARG_BOOL(3);
1891 
1892  /* SRIDs of raster and geometry must match */
1894  elog(NOTICE, "SRIDs of geometry and raster do not match");
1895  rt_raster_destroy(raster);
1896  PG_FREE_IF_COPY(pgraster, 0);
1897  PG_FREE_IF_COPY(geom, 2);
1898  PG_RETURN_NULL();
1899  }
1900 
1901  /* get band */
1902  band = rt_raster_get_band(raster, bandindex - 1);
1903  if (!band) {
1904  elog(NOTICE, "Could not find band at index %d. Returning NULL", bandindex);
1905  rt_raster_destroy(raster);
1906  PG_FREE_IF_COPY(pgraster, 0);
1907  PG_FREE_IF_COPY(geom, 2);
1908  PG_RETURN_NULL();
1909  }
1910 
1911  /* process geometry */
1912  lwgeom = lwgeom_from_gserialized(geom);
1913 
1914  if (lwgeom_is_empty(lwgeom)) {
1915  elog(NOTICE, "Geometry provided cannot be empty");
1916  rt_raster_destroy(raster);
1917  PG_FREE_IF_COPY(pgraster, 0);
1918  PG_FREE_IF_COPY(geom, 2);
1919  PG_RETURN_NULL();
1920  }
1921 
1922  /* Get a 2D version of the geometry if necessary */
1923  if (lwgeom_ndims(lwgeom) > 2) {
1924  LWGEOM *lwgeom2d = lwgeom_force_2d(lwgeom);
1925  lwgeom_free(lwgeom);
1926  lwgeom = lwgeom2d;
1927  }
1928 
1929  point = lwgeom_as_lwpoint(lwgeom);
1930  getPoint2d_p(point->point, 0, &p);
1931 
1933  raster,
1934  p.x, p.y,
1935  &x, &y,
1936  NULL
1937  ) != ES_NONE) {
1938  rt_raster_destroy(raster);
1939  PG_FREE_IF_COPY(pgraster, 0);
1940  lwgeom_free(lwgeom);
1941  PG_FREE_IF_COPY(geom, 2);
1942  elog(ERROR, "RASTER_nearestValue: Could not compute pixel coordinates from spatial coordinates");
1943  PG_RETURN_NULL();
1944  }
1945 
1946  /* get value at point */
1947  if (
1948  (x >= 0 && x < rt_raster_get_width(raster)) &&
1949  (y >= 0 && y < rt_raster_get_height(raster))
1950  ) {
1951  if (rt_band_get_pixel(band, x, y, &value, &isnodata) != ES_NONE) {
1952  rt_raster_destroy(raster);
1953  PG_FREE_IF_COPY(pgraster, 0);
1954  lwgeom_free(lwgeom);
1955  PG_FREE_IF_COPY(geom, 2);
1956  elog(ERROR, "RASTER_nearestValue: Could not get pixel value for band at index %d", bandindex);
1957  PG_RETURN_NULL();
1958  }
1959 
1960  /* value at point, return value */
1961  if (!exclude_nodata_value || !isnodata) {
1962  rt_raster_destroy(raster);
1963  PG_FREE_IF_COPY(pgraster, 0);
1964  lwgeom_free(lwgeom);
1965  PG_FREE_IF_COPY(geom, 2);
1966 
1967  PG_RETURN_FLOAT8(value);
1968  }
1969  }
1970 
1971  /* get neighborhood */
1972  count = rt_band_get_nearest_pixel(
1973  band,
1974  x, y,
1975  0, 0,
1976  exclude_nodata_value,
1977  &npixels
1978  );
1979  rt_band_destroy(band);
1980  /* error or no neighbors */
1981  if (count < 1) {
1982  /* error */
1983  if (count < 0)
1984  elog(NOTICE, "Could not get the nearest value for band at index %d", bandindex);
1985  /* no nearest pixel */
1986  else
1987  elog(NOTICE, "No nearest value found for band at index %d", bandindex);
1988 
1989  lwgeom_free(lwgeom);
1990  PG_FREE_IF_COPY(geom, 2);
1991  rt_raster_destroy(raster);
1992  PG_FREE_IF_COPY(pgraster, 0);
1993  PG_RETURN_NULL();
1994  }
1995 
1996  /* more than one nearest value, see which one is closest */
1997  if (count > 1) {
1998  int i = 0;
1999  LWPOLY *poly = NULL;
2000  double lastdist = -1;
2001  double dist;
2002 
2003  for (i = 0; i < count; i++) {
2004  /* convex-hull of pixel */
2005  poly = rt_raster_pixel_as_polygon(raster, npixels[i].x, npixels[i].y);
2006  if (!poly) {
2007  lwgeom_free(lwgeom);
2008  PG_FREE_IF_COPY(geom, 2);
2009  rt_raster_destroy(raster);
2010  PG_FREE_IF_COPY(pgraster, 0);
2011  elog(ERROR, "RASTER_nearestValue: Could not get polygon of neighboring pixel");
2012  PG_RETURN_NULL();
2013  }
2014 
2015  /* distance between convex-hull and point */
2016  dist = lwgeom_mindistance2d(lwpoly_as_lwgeom(poly), lwgeom);
2017  if (lastdist < 0 || dist < lastdist) {
2018  value = npixels[i].value;
2019  hasvalue = 1;
2020  }
2021  lastdist = dist;
2022 
2023  lwpoly_free(poly);
2024  }
2025  }
2026  else {
2027  value = npixels[0].value;
2028  hasvalue = 1;
2029  }
2030 
2031  pfree(npixels);
2032  lwgeom_free(lwgeom);
2033  PG_FREE_IF_COPY(geom, 2);
2034  rt_raster_destroy(raster);
2035  PG_FREE_IF_COPY(pgraster, 0);
2036 
2037  if (hasvalue)
2038  PG_RETURN_FLOAT8(value);
2039  else
2040  PG_RETURN_NULL();
2041 }
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:347
uint16_t 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:1144
band
Definition: ovdump.py:57
double value
Definition: librtcore.h:2318
void rt_band_destroy(rt_band band)
Destroy a raster band.
Definition: rt_band.c:336
LWPOINT * lwgeom_as_lwpoint(const LWGEOM *lwgeom)
Definition: lwgeom.c:161
POINTARRAY * point
Definition: liblwgeom.h:413
LWGEOM * lwpoly_as_lwgeom(const LWPOLY *obj)
Definition: lwgeom.c:320
int lwgeom_ndims(const LWGEOM *geom)
Return the number of dimensions (2, 3, 4) in a geometry.
Definition: lwgeom.c:944
double x
Definition: liblwgeom.h:330
void lwpoly_free(LWPOLY *poly)
Definition: lwpoly.c:175
rt_errorstate rt_band_get_pixel(rt_band band, int x, int y, double *value, int *nodata)
Get pixel value.
Definition: rt_band.c:1173
LWGEOM * lwgeom_force_2d(const LWGEOM *geom)
Strip out the Z/M components of an LWGEOM.
Definition: lwgeom.c:784
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:330
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
int getPoint2d_p(const POINTARRAY *pa, uint32_t n, POINT2D *point)
Definition: lwgeom_api.c:338
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:2230
#define POINTTYPE
LWTYPE numbers, used internally by PostGIS.
Definition: liblwgeom.h:84
double lwgeom_mindistance2d(const LWGEOM *lw1, const LWGEOM *lw2)
Function initializing 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:1393
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:1326
#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:99
Here is the call graph for this function:
Here is the caller graph for this function: