PostGIS  2.1.10dev-r@@SVN_REVISION@@
Datum RASTER_nearestValue ( PG_FUNCTION_ARGS  )

Definition at line 4858 of file rt_pg.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(), 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, genraster::value, rt_pixel_t::value, pixval::x, POINT2D::x, pixval::y, and POINT2D::y.

4859 {
4860  rt_pgraster *pgraster = NULL;
4861  rt_raster raster = NULL;
4862  rt_band band = NULL;
4863  int bandindex = 1;
4864  int num_bands = 0;
4865  GSERIALIZED *geom;
4866  bool exclude_nodata_value = TRUE;
4867  LWGEOM *lwgeom;
4868  LWPOINT *point = NULL;
4869  POINT2D p;
4870 
4871  double x;
4872  double y;
4873  int count;
4874  rt_pixel npixels = NULL;
4875  double value = 0;
4876  int hasvalue = 0;
4877  int isnodata = 0;
4878 
4879  if (PG_ARGISNULL(0))
4880  PG_RETURN_NULL();
4881  pgraster = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
4882  raster = rt_raster_deserialize(pgraster, FALSE);
4883  if (!raster) {
4884  PG_FREE_IF_COPY(pgraster, 0);
4885  elog(ERROR, "RASTER_nearestValue: Could not deserialize raster");
4886  PG_RETURN_NULL();
4887  }
4888 
4889  /* band index is 1-based */
4890  if (!PG_ARGISNULL(1))
4891  bandindex = PG_GETARG_INT32(1);
4892  num_bands = rt_raster_get_num_bands(raster);
4893  if (bandindex < 1 || bandindex > num_bands) {
4894  elog(NOTICE, "Invalid band index (must use 1-based). Returning NULL");
4895  rt_raster_destroy(raster);
4896  PG_FREE_IF_COPY(pgraster, 0);
4897  PG_RETURN_NULL();
4898  }
4899 
4900  /* point */
4901  geom = (GSERIALIZED *) PG_DETOAST_DATUM(PG_GETARG_DATUM(2));
4902  if (gserialized_get_type(geom) != POINTTYPE) {
4903  elog(NOTICE, "Geometry provided must be a point");
4904  rt_raster_destroy(raster);
4905  PG_FREE_IF_COPY(pgraster, 0);
4906  PG_FREE_IF_COPY(geom, 2);
4907  PG_RETURN_NULL();
4908  }
4909 
4910  /* exclude_nodata_value flag */
4911  if (!PG_ARGISNULL(3))
4912  exclude_nodata_value = PG_GETARG_BOOL(3);
4913 
4914  /* SRIDs of raster and geometry must match */
4916  elog(NOTICE, "SRIDs of geometry and raster do not match");
4917  rt_raster_destroy(raster);
4918  PG_FREE_IF_COPY(pgraster, 0);
4919  PG_FREE_IF_COPY(geom, 2);
4920  PG_RETURN_NULL();
4921  }
4922 
4923  /* get band */
4924  band = rt_raster_get_band(raster, bandindex - 1);
4925  if (!band) {
4926  elog(NOTICE, "Could not find band at index %d. Returning NULL", bandindex);
4927  rt_raster_destroy(raster);
4928  PG_FREE_IF_COPY(pgraster, 0);
4929  PG_FREE_IF_COPY(geom, 2);
4930  PG_RETURN_NULL();
4931  }
4932 
4933  /* process geometry */
4934  lwgeom = lwgeom_from_gserialized(geom);
4935 
4936  if (lwgeom_is_empty(lwgeom)) {
4937  elog(NOTICE, "Geometry provided cannot be empty");
4938  rt_raster_destroy(raster);
4939  PG_FREE_IF_COPY(pgraster, 0);
4940  PG_FREE_IF_COPY(geom, 2);
4941  PG_RETURN_NULL();
4942  }
4943 
4944  /* Get a 2D version of the geometry if necessary */
4945  if (lwgeom_ndims(lwgeom) > 2) {
4946  LWGEOM *lwgeom2d = lwgeom_force_2d(lwgeom);
4947  lwgeom_free(lwgeom);
4948  lwgeom = lwgeom2d;
4949  }
4950 
4951  point = lwgeom_as_lwpoint(lwgeom);
4952  getPoint2d_p(point->point, 0, &p);
4953 
4955  raster,
4956  p.x, p.y,
4957  &x, &y,
4958  NULL
4959  ) != ES_NONE) {
4960  rt_raster_destroy(raster);
4961  PG_FREE_IF_COPY(pgraster, 0);
4962  lwgeom_free(lwgeom);
4963  PG_FREE_IF_COPY(geom, 2);
4964  elog(ERROR, "RASTER_nearestValue: Could not compute pixel coordinates from spatial coordinates");
4965  PG_RETURN_NULL();
4966  }
4967 
4968  /* get value at point */
4969  if (
4970  (x >= 0 && x < rt_raster_get_width(raster)) &&
4971  (y >= 0 && y < rt_raster_get_height(raster))
4972  ) {
4973  if (rt_band_get_pixel(band, x, y, &value, &isnodata) != ES_NONE) {
4974  rt_raster_destroy(raster);
4975  PG_FREE_IF_COPY(pgraster, 0);
4976  lwgeom_free(lwgeom);
4977  PG_FREE_IF_COPY(geom, 2);
4978  elog(ERROR, "RASTER_nearestValue: Could not get pixel value for band at index %d", bandindex);
4979  PG_RETURN_NULL();
4980  }
4981 
4982  /* value at point, return value */
4983  if (!exclude_nodata_value || !isnodata) {
4984  rt_raster_destroy(raster);
4985  PG_FREE_IF_COPY(pgraster, 0);
4986  lwgeom_free(lwgeom);
4987  PG_FREE_IF_COPY(geom, 2);
4988 
4989  PG_RETURN_FLOAT8(value);
4990  }
4991  }
4992 
4993  /* get neighborhood */
4994  count = rt_band_get_nearest_pixel(
4995  band,
4996  x, y,
4997  0, 0,
4998  exclude_nodata_value,
4999  &npixels
5000  );
5001  rt_band_destroy(band);
5002  /* error or no neighbors */
5003  if (count < 1) {
5004  /* error */
5005  if (count < 0)
5006  elog(NOTICE, "Could not get the nearest value for band at index %d", bandindex);
5007  /* no nearest pixel */
5008  else
5009  elog(NOTICE, "No nearest value found for band at index %d", bandindex);
5010 
5011  lwgeom_free(lwgeom);
5012  PG_FREE_IF_COPY(geom, 2);
5013  rt_raster_destroy(raster);
5014  PG_FREE_IF_COPY(pgraster, 0);
5015  PG_RETURN_NULL();
5016  }
5017 
5018  /* more than one nearest value, see which one is closest */
5019  if (count > 1) {
5020  int i = 0;
5021  LWPOLY *poly = NULL;
5022  double lastdist = -1;
5023  double dist;
5024 
5025  for (i = 0; i < count; i++) {
5026  /* convex-hull of pixel */
5027  poly = rt_raster_pixel_as_polygon(raster, npixels[i].x, npixels[i].y);
5028  if (!poly) {
5029  lwgeom_free(lwgeom);
5030  PG_FREE_IF_COPY(geom, 2);
5031  rt_raster_destroy(raster);
5032  PG_FREE_IF_COPY(pgraster, 0);
5033  elog(ERROR, "RASTER_nearestValue: Could not get polygon of neighboring pixel");
5034  PG_RETURN_NULL();
5035  }
5036 
5037  /* distance between convex-hull and point */
5038  dist = lwgeom_mindistance2d(lwpoly_as_lwgeom(poly), lwgeom);
5039  if (lastdist < 0 || dist < lastdist) {
5040  value = npixels[i].value;
5041  hasvalue = 1;
5042  }
5043  lastdist = dist;
5044 
5045  lwpoly_free(poly);
5046  }
5047  }
5048  else {
5049  value = npixels[0].value;
5050  hasvalue = 1;
5051  }
5052 
5053  pfree(npixels);
5054  lwgeom_free(lwgeom);
5055  PG_FREE_IF_COPY(geom, 2);
5056  rt_raster_destroy(raster);
5057  PG_FREE_IF_COPY(pgraster, 0);
5058 
5059  if (hasvalue)
5060  PG_RETURN_FLOAT8(value);
5061  else
5062  PG_RETURN_NULL();
5063 }
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:56
LWPOLY * rt_raster_pixel_as_polygon(rt_raster rast, int x, int y)
Get a raster pixel as a polygon.
Definition: rt_api.c:13182
int rt_raster_get_num_bands(rt_raster raster)
Definition: rt_api.c:5677
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:326
void rt_raster_destroy(rt_raster raster)
Release memory associated to a raster.
Definition: rt_api.c:5387
LWGEOM * lwgeom_from_gserialized(const GSERIALIZED *g)
Allocate a new LWGEOM from a GSERIALIZED.
void lwgeom_free(LWGEOM *geom)
Definition: lwgeom.c:1006
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_api.c:2702
tuple band
Definition: ovdump.py:57
double value
Definition: rt_api.h:2263
tuple raster
Be careful!! Zeros function's input parameter can be a (height x width) array, not (width x height): ...
Definition: rtrowdump.py:123
int32_t rt_raster_get_srid(rt_raster raster)
Get raster's SRID.
Definition: rt_api.c:5661
LWPOINT * lwgeom_as_lwpoint(const LWGEOM *lwgeom)
Definition: lwgeom.c:80
POINTARRAY * point
Definition: liblwgeom.h:367
LWGEOM * lwpoly_as_lwgeom(const LWPOLY *obj)
Definition: lwgeom.c:239
int lwgeom_ndims(const LWGEOM *geom)
Return the number of dimensions (2, 3, 4) in a geometry.
Definition: lwgeom.c:806
double x
Definition: liblwgeom.h:284
double lwgeom_mindistance2d(LWGEOM *lw1, LWGEOM *lw2)
Function initialazing min distance calculation.
Definition: measures.c:162
void lwpoly_free(LWPOLY *poly)
Definition: lwpoly.c:79
rt_errorstate rt_band_get_pixel(rt_band band, int x, int y, double *value, int *nodata)
Get pixel value.
Definition: rt_api.c:2549
LWGEOM * lwgeom_force_2d(const LWGEOM *geom)
Strip out the Z/M components of an LWGEOM.
Definition: lwgeom.c:646
int count
Definition: genraster.py:57
double y
Definition: liblwgeom.h:284
int getPoint2d_p(const POINTARRAY *pa, int n, POINT2D *point)
Definition: lwgeom_api.c:434
tuple x
Definition: pixval.py:53
void rt_band_destroy(rt_band band)
Destroy a raster band.
Definition: rt_api.c:1650
uint16_t rt_raster_get_height(rt_raster raster)
Definition: rt_api.c:5434
rt_band rt_raster_get_band(rt_raster raster, int n)
Return Nth band, or NULL if unavailable.
Definition: rt_api.c:5686
#define FALSE
Definition: dbfopen.c:169
Struct definitions.
Definition: rt_api.h:2175
#define POINTTYPE
LWTYPE numbers, used internally by PostGIS.
Definition: liblwgeom.h:60
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_api.c:6105
rt_raster rt_raster_deserialize(void *serialized, int header_only)
Return a raster from a serialized form.
Definition: rt_api.c:8350
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:1229
#define TRUE
Definition: dbfopen.c:170
tuple y
Definition: pixval.py:54
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:70
uint16_t rt_raster_get_width(rt_raster raster)
Definition: rt_api.c:5426

Here is the call graph for this function: