}
}
+/**
+ * Get nearest pixel(s) with value (not NODATA) to specified pixel
+ *
+ * @param band: the band to get nearest pixel(s) from
+ * @param x: the column of the pixel (0-based)
+ * @param y: the line of the pixel (0-based)
+ * @param distance: the number of pixels around the specified pixel
+ * @param exclude_nodata_value: if non-zero, ignore nodata values
+ * to check for pixels with value
+ * @param npixels: return set of rt_pixel object or NULL
+ *
+ * @return -1 on error, otherwise the number of rt_pixel objects
+ * in npixels
+ */
+int rt_band_get_nearest_pixel(
+ rt_band band,
+ int x, int y,
+ uint16_t distance,
+ int exclude_nodata_value,
+ rt_pixel *npixels
+) {
+ rt_pixel npixel = NULL;
+ int extent[4] = {0};
+ int d0 = 0;
+ uint32_t _distance;
+ uint32_t i = 0;
+ uint32_t j = 0;
+ uint32_t k = 0;
+ uint32_t max = 0;
+ int _x = 0;
+ int _y = 0;
+ int *_inc = NULL;
+ double pixval = 0;
+ double minval = 0;
+ int err;
+ uint32_t count = 0;
+
+ assert(NULL != band);
+ assert(NULL != npixels);
+
+ RASTER_DEBUG(3, "Starting");
+
+ minval = rt_pixtype_get_min_value(band->pixtype);
+
+ if (!distance)
+ d0 = 1;
+
+ RASTER_DEBUGF(4, "Selected pixel: %d x %d", x, y);
+ RASTER_DEBUGF(4, "Distance: %d", distance);
+
+ /* shortcuts if outside band extent */
+ if (
+ exclude_nodata_value && (
+ (x < 0 || x > band->width) ||
+ (y < 0 || y > band->height)
+ )
+ ) {
+ /* no distance specified, jump to pixel close to extent */
+ if (!distance) {
+ if (x < 0)
+ x = -1;
+ else if (x > band->width)
+ x = band->width;
+
+ if (y < 0)
+ y = -1;
+ else if (y > band->height)
+ y = band->height;
+
+ RASTER_DEBUGF(4, "Moved selected pixel: %d x %d", x, y);
+ }
+ /*
+ distance specified
+ if distance won't capture extent of band, return 0
+ */
+ else if (
+ ((x < 0 && abs(x) > distance) || (x - band->width >= distance)) ||
+ ((x < 0 && abs(y) > distance) || (y - band->height >= distance))
+ ) {
+ RASTER_DEBUG(4, "No nearest pixels possible for provided pixel and distance");
+ return 0;
+ }
+ }
+
+ /* determine the maximum distance to prevent an infinite loop */
+ if (!distance) {
+ int a, b;
+ a = abs(x);
+ b = abs(x - band->width);
+
+ if (a > b)
+ distance = a;
+ else
+ distance = b;
+
+ a = abs(y);
+ b = abs(y - band->height);
+ if (a > distance)
+ distance = a;
+ if (b > distance)
+ distance = b;
+
+ RASTER_DEBUGF(4, "Maximum distance: %d", distance);
+ }
+
+ count = 0;
+ _distance = 0;
+ *npixels = NULL;
+ do {
+ _distance++;
+ extent[0] = x - _distance; /* min x */
+ extent[1] = y - _distance; /* min y */
+ extent[2] = x + _distance; /* max x */
+ extent[3] = y + _distance; /* max y */
+
+ RASTER_DEBUGF(4, "Processing distance: %d", _distance);
+ RASTER_DEBUGF(4, "Extent: (%d, %d, %d, %d)",
+ extent[0], extent[1], extent[2], extent[3]);
+
+ for (i = 0; i < 2; i++) {
+
+ /* by row */
+ if (i < 1)
+ max = extent[2] - extent[0] + 1;
+ /* by column */
+ else
+ max = extent[3] - extent[1] + 1;
+ max = abs(max);
+
+ for (j = 0; j < 2; j++) {
+ /* by row */
+ if (i < 1) {
+ _x = extent[0];
+ _inc = &_x;
+
+ /* top row */
+ if (j < 1)
+ _y = extent[1];
+ /* bottom row */
+ else
+ _y = extent[3];
+ }
+ /* by column */
+ else {
+ _y = extent[1] + 1;
+ _inc = &_y;
+
+ /* left column */
+ if (j < 1) {
+ _x = extent[0];
+ max -= 2;
+ }
+ /* right column */
+ else
+ _x = extent[2];
+ }
+
+ RASTER_DEBUGF(4, "max: %d", max);
+ for (k = 0; k < max; k++) {
+ /* outside band extent, set to NODATA */
+ if (
+ (_x < 0 || _x >= band->width) ||
+ (_y < 0 || _y >= band->height)
+ ) {
+ /* no NODATA, set to minimum possible value */
+ if (!band->hasnodata)
+ pixval = minval;
+ else
+ pixval = band->nodataval;
+ RASTER_DEBUGF(4, "NODATA pixel outside band extent: (x, y, val) = (%d, %d, %f)", _x, _y, pixval);
+ }
+ else {
+ err = rt_band_get_pixel(
+ band,
+ _x, _y,
+ &pixval
+ );
+ if (err < 0) {
+ rterror("rt_band_get_nearest_pixel: Unable to get pixel value");
+ if (count) rtdealloc(*npixels);
+ return -1;
+ }
+ RASTER_DEBUGF(4, "Pixel: (x, y, val) = (%d, %d, %f)", _x, _y, pixval);
+ }
+
+ /* use pixval? */
+ if (
+ !exclude_nodata_value || (
+ exclude_nodata_value &&
+ (band->hasnodata != FALSE) && (
+ FLT_NEQ(pixval, band->nodataval) &&
+ (rt_band_clamped_value_is_nodata(band, pixval) != 1)
+ )
+ )
+ ) {
+ /* add pixel to result set */
+ RASTER_DEBUGF(4, "Adding pixel to set of nearest pixels: (x, y, val) = (%d, %d, %f)", _x, _y, pixval);
+ count++;
+
+ if (*npixels == NULL)
+ *npixels = (rt_pixel) rtalloc(sizeof(struct rt_pixel_t) * count);
+ else
+ *npixels = (rt_pixel) rtrealloc(*npixels, sizeof(struct rt_pixel_t) * count);
+ if (*npixels == NULL) {
+ rterror("rt_band_get_nearest_pixel: Unable to allocate memory for nearest pixel(s)");
+ return -1;
+ }
+
+ npixel = &((*npixels)[count - 1]);
+ npixel->x = _x;
+ npixel->y = _y;
+ npixel->value = pixval;
+ }
+
+ (*_inc)++;
+ }
+ }
+ }
+
+ /* distance threshhold met */
+ if (_distance >= distance)
+ break;
+ else if (d0 && count)
+ break;
+ }
+ while (1);
+
+ RASTER_DEBUGF(3, "Nearest pixels in return: %d", count);
+
+ return count;
+}
+
double
rt_band_get_nodata(rt_band band) {
*/
typedef struct rt_raster_t* rt_raster;
typedef struct rt_band_t* rt_band;
+typedef struct rt_pixel_t* rt_pixel;
typedef struct rt_geomval_t* rt_geomval;
typedef struct rt_bandstats_t* rt_bandstats;
typedef struct rt_histogram_t* rt_histogram;
double *result
);
+/**
+ * Get nearest pixel(s) with value (not NODATA) to specified pixel
+ *
+ * @param band: the band to get nearest pixel(s) from
+ * @param x: the column of the pixel (0-based)
+ * @param y: the line of the pixel (0-based)
+ * @param distance: the number of pixels around the specified pixel
+ * @param exclude_nodata_value: if non-zero, ignore nodata values
+ * to check for pixels with value
+ * @param npixels: return set of rt_pixel object or NULL
+ *
+ * @return -1 on error, otherwise the number of rt_pixel objects
+ * in npixels
+ */
+int rt_band_get_nearest_pixel(
+ rt_band band,
+ int x, int y,
+ uint16_t distance,
+ int exclude_nodata_value,
+ rt_pixel *npixels
+);
+
/**
* Returns the minimal possible value for the band according to the pixel type.
* @param band: the band to get info from
};
+struct rt_pixel_t {
+ int x; /* column */
+ int y; /* line */
+ double value;
+};
+
/* polygon as LWPOLY with associated value */
struct rt_geomval_t {
LWPOLY *geom;
/* Get pixel geographical shape */
Datum RASTER_getPixelPolygon(PG_FUNCTION_ARGS);
+/* Get nearest value to a point */
+Datum RASTER_nearestValue(PG_FUNCTION_ARGS);
+
+/* Get the neighborhood around a pixel */
+Datum RASTER_neighborhood(PG_FUNCTION_ARGS);
+
/* Raster and band creation */
Datum RASTER_makeEmpty(PG_FUNCTION_ARGS);
Datum RASTER_addband(PG_FUNCTION_ARGS);
PG_RETURN_POINTER(gser);
}
+/**
+ * Return nearest value to a point
+ */
+PG_FUNCTION_INFO_V1(RASTER_nearestValue);
+Datum RASTER_nearestValue(PG_FUNCTION_ARGS)
+{
+ rt_pgraster *pgraster = NULL;
+ rt_raster raster = NULL;
+ rt_band band = NULL;
+ int bandindex = 1;
+ int num_bands = 0;
+ GSERIALIZED *geom;
+ bool exclude_nodata_value = TRUE;
+ LWGEOM *lwgeom;
+ LWPOINT *point = NULL;
+ POINT2D p;
+
+ double x;
+ double y;
+ int count;
+ rt_pixel npixels = NULL;
+ double value = 0;
+ int hasvalue = 0;
+
+ if (PG_ARGISNULL(0))
+ PG_RETURN_NULL();
+ pgraster = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
+ raster = rt_raster_deserialize(pgraster, FALSE);
+ if (!raster) {
+ elog(ERROR, "RASTER_nearestValue: Could not deserialize raster");
+ PG_RETURN_NULL();
+ }
+
+ /* band index is 1-based */
+ if (!PG_ARGISNULL(1))
+ bandindex = PG_GETARG_INT32(1);
+ num_bands = rt_raster_get_num_bands(raster);
+ if (bandindex < 1 || bandindex > num_bands) {
+ elog(NOTICE, "Invalid band index (must use 1-based). Returning NULL");
+ rt_raster_destroy(raster);
+ PG_RETURN_NULL();
+ }
+
+ /* point */
+ geom = (GSERIALIZED *) PG_DETOAST_DATUM(PG_GETARG_DATUM(2));
+ if (gserialized_get_type(geom) != POINTTYPE) {
+ elog(NOTICE, "Geometry provided must be a point");
+ rt_raster_destroy(raster);
+ PG_RETURN_NULL();
+ }
+
+ /* exclude_nodata_value flag */
+ if (!PG_ARGISNULL(3))
+ exclude_nodata_value = PG_GETARG_BOOL(3);
+
+ /* get band */
+ band = rt_raster_get_band(raster, bandindex - 1);
+ if (!band) {
+ elog(NOTICE, "Could not find band at index %d. Returning NULL", bandindex);
+ rt_raster_destroy(raster);
+ PG_RETURN_NULL();
+ }
+
+ /* process geometry */
+ lwgeom = lwgeom_from_gserialized(geom);
+
+ if (lwgeom_is_empty(lwgeom)) {
+ elog(NOTICE, "Geometry provided cannot be empty");
+ rt_raster_destroy(raster);
+ PG_RETURN_NULL();
+ }
+
+ point = lwgeom_as_lwpoint(lwgeom);
+ getPoint2d_p(point->point, 0, &p);
+
+ if (!rt_raster_geopoint_to_cell(
+ raster,
+ p.x, p.y,
+ &x, &y,
+ NULL
+ )) {
+ elog(ERROR, "RASTER_nearestValue: Unable to compute pixel coordinates from spatial coordinates");
+ rt_raster_destroy(raster);
+ lwgeom_free(lwgeom);
+ PG_RETURN_NULL();
+ }
+
+ /* get value at point */
+ if (
+ (x >= 0 && x < rt_raster_get_width(raster)) &&
+ (y >= 0 && y < rt_raster_get_height(raster))
+ ) {
+ if (rt_band_get_pixel(band, x, y, &value) < 0) {
+ elog(ERROR, "RASTER_nearestValue: Unable to get pixel value for band at index %d", bandindex);
+ rt_raster_destroy(raster);
+ lwgeom_free(lwgeom);
+ PG_RETURN_NULL();
+ }
+
+ /* value at point, return value */
+ if (
+ !exclude_nodata_value || (
+ exclude_nodata_value &&
+ rt_band_get_hasnodata_flag(band) && (
+ FLT_NEQ(value, rt_band_get_nodata(band)) &&
+ (rt_band_clamped_value_is_nodata(band, value) != 1)
+ )
+ )
+ ) {
+ rt_raster_destroy(raster);
+ lwgeom_free(lwgeom);
+
+ PG_RETURN_FLOAT8(value);
+ }
+ }
+
+ /* get neighborhood */
+ count = rt_band_get_nearest_pixel(
+ band,
+ x, y,
+ 0,
+ exclude_nodata_value,
+ &npixels
+ );
+ rt_band_destroy(band);
+ /* error or no neighbors */
+ if (count < 1) {
+ /* error */
+ if (count < 0)
+ elog(NOTICE, "Unable to get the nearest value for band at index %d", bandindex);
+ /* no nearest pixel */
+ else
+ elog(NOTICE, "No nearest value found for band at index %d", bandindex);
+
+ lwgeom_free(lwgeom);
+ rt_raster_destroy(raster);
+ PG_RETURN_NULL();
+ }
+
+ /* more than one nearest value, see which one is closest */
+ if (count > 1) {
+ int i = 0;
+ LWPOLY *poly = NULL;
+ double lastdist = -1;
+ double dist;
+
+ for (i = 0; i < count; i++) {
+ /* convex-hull of pixel */
+ poly = rt_raster_pixel_as_polygon(raster, npixels[i].x, npixels[i].y);
+ if (!poly) {
+ elog(ERROR, "RASTER_nearestValue: Unable to get polygon of neighboring pixel");
+ lwgeom_free(lwgeom);
+ rt_raster_destroy(raster);
+ PG_RETURN_NULL();
+ }
+
+ /* distance between convex-hull and point */
+ dist = lwgeom_mindistance2d(lwpoly_as_lwgeom(poly), lwgeom);
+ if (lastdist < 0 || dist < lastdist) {
+ value = npixels[i].value;
+ hasvalue = 1;
+ }
+ lastdist = dist;
+
+ lwpoly_free(poly);
+ }
+ }
+ else {
+ value = npixels[0].value;
+ hasvalue = 1;
+ }
+
+ pfree(npixels);
+ lwgeom_free(lwgeom);
+ rt_raster_destroy(raster);
+
+ if (hasvalue)
+ PG_RETURN_FLOAT8(value);
+ else
+ PG_RETURN_NULL();
+}
+
+/**
+ * Return the neighborhood around a pixel
+ */
+PG_FUNCTION_INFO_V1(RASTER_neighborhood);
+Datum RASTER_neighborhood(PG_FUNCTION_ARGS)
+{
+ FuncCallContext *funcctx;
+ TupleDesc tupdesc;
+
+ int call_cntr;
+ int max_calls;
+
+ rt_pixel npixel1 = NULL;
+ rt_pixel npixel2 = NULL;
+
+ if (SRF_IS_FIRSTCALL()) {
+ MemoryContext oldcontext;
+
+ rt_pgraster *pgraster = NULL;
+ rt_raster raster = NULL;
+ rt_band band = NULL;
+ int bandindex = 1;
+ int num_bands = 0;
+ int x = 0;
+ int y = 0;
+ int distance = 0;
+ bool exclude_nodata_value = TRUE;
+
+ int count;
+
+ /* create a function context for cross-call persistence */
+ funcctx = SRF_FIRSTCALL_INIT();
+
+ /* switch to memory context appropriate for multiple function calls */
+ oldcontext = MemoryContextSwitchTo(funcctx->multi_call_memory_ctx);
+
+ /* pgraster is null, return nothing */
+ if (PG_ARGISNULL(0)) {
+ MemoryContextSwitchTo(oldcontext);
+ SRF_RETURN_DONE(funcctx);
+ }
+ pgraster = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
+
+ raster = rt_raster_deserialize(pgraster, FALSE);
+ if (!raster) {
+ elog(ERROR, "RASTER_neighborhood: Could not deserialize raster");
+ MemoryContextSwitchTo(oldcontext);
+ SRF_RETURN_DONE(funcctx);
+ }
+
+ /* band index is 1-based */
+ if (!PG_ARGISNULL(1))
+ bandindex = PG_GETARG_INT32(1);
+ num_bands = rt_raster_get_num_bands(raster);
+ if (bandindex < 1 || bandindex > num_bands) {
+ elog(NOTICE, "Invalid band index (must use 1-based). Returning NULL");
+ rt_raster_destroy(raster);
+ MemoryContextSwitchTo(oldcontext);
+ SRF_RETURN_DONE(funcctx);
+ }
+
+ /* pixel column, 1-based */
+ x = PG_GETARG_INT32(2);
+
+ /* pixel row, 1-based */
+ y = PG_GETARG_INT32(3);
+
+ /* distance */
+ distance = PG_GETARG_INT32(4);
+ if (distance < 1) {
+ elog(NOTICE, "Invalid value for distance (must be greater than zero). Returning NULL");
+ rt_raster_destroy(raster);
+ MemoryContextSwitchTo(oldcontext);
+ SRF_RETURN_DONE(funcctx);
+ }
+
+ /* exclude_nodata_value flag */
+ if (!PG_ARGISNULL(5))
+ exclude_nodata_value = PG_GETARG_BOOL(5);
+
+ /* get band */
+ band = rt_raster_get_band(raster, bandindex - 1);
+ if (!band) {
+ elog(NOTICE, "Could not find band at index %d. Returning NULL", bandindex);
+ rt_raster_destroy(raster);
+ MemoryContextSwitchTo(oldcontext);
+ SRF_RETURN_DONE(funcctx);
+ }
+
+ /* get neighborhood */
+ count = rt_band_get_nearest_pixel(
+ band,
+ x - 1, y - 1,
+ (uint16_t) distance,
+ exclude_nodata_value,
+ &npixel1
+ );
+ rt_band_destroy(band);
+ rt_raster_destroy(raster);
+ /* error or no neighbors */
+ if (count < 1) {
+ /* error */
+ if (count < 0)
+ elog(NOTICE, "Unable to get the pixel's neighborhood for band at index %d", bandindex);
+ /* no neighbors */
+ else
+ elog(NOTICE, "Pixel has no neighbors for band at distance %d", distance);
+
+ MemoryContextSwitchTo(oldcontext);
+ SRF_RETURN_DONE(funcctx);
+ }
+
+ /* Store needed information */
+ funcctx->user_fctx = npixel1;
+
+ /* total number of tuples to be returned */
+ funcctx->max_calls = count;
+
+ /* Build a tuple descriptor for our result type */
+ if (get_call_result_type(fcinfo, NULL, &tupdesc) != TYPEFUNC_COMPOSITE) {
+ ereport(ERROR, (
+ errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
+ errmsg(
+ "function returning record called in context "
+ "that cannot accept type record"
+ )
+ ));
+ }
+
+ BlessTupleDesc(tupdesc);
+ funcctx->tuple_desc = tupdesc;
+
+ MemoryContextSwitchTo(oldcontext);
+ }
+
+ /* stuff done on every call of the function */
+ funcctx = SRF_PERCALL_SETUP();
+
+ call_cntr = funcctx->call_cntr;
+ max_calls = funcctx->max_calls;
+ tupdesc = funcctx->tuple_desc;
+ npixel2 = funcctx->user_fctx;
+
+ /* do when there is more left to send */
+ if (call_cntr < max_calls) {
+ int values_length = 3;
+ Datum values[values_length];
+ bool *nulls = NULL;
+ HeapTuple tuple;
+ Datum result;
+
+ POSTGIS_RT_DEBUGF(3, "Result %d", call_cntr);
+
+ nulls = palloc(sizeof(bool) * values_length);
+ memset(nulls, FALSE, values_length);
+
+ /* x,y are 0-based, make 1-based for end users */
+ values[0] = Int64GetDatum(npixel2[call_cntr].x + 1);
+ values[1] = Int64GetDatum(npixel2[call_cntr].y + 1);
+ values[2] = Float8GetDatum(npixel2[call_cntr].value);
+
+ /* build a tuple */
+ tuple = heap_form_tuple(tupdesc, values, nulls);
+
+ /* make the tuple into a datum */
+ result = HeapTupleGetDatum(tuple);
+
+ /* clean up */
+ pfree(nulls);
+
+ SRF_RETURN_NEXT(funcctx, result);
+ }
+ /* do when there is no more left */
+ else {
+ pfree(npixel2);
+ SRF_RETURN_DONE(funcctx);
+ }
+
+}
+
/**
* Add a band to the given raster at the given position.
*/
$$ SELECT ST_Clip($1, NULL, $2, null::float8[], $3) $$
LANGUAGE 'sql' STABLE;
+-----------------------------------------------------------------------
+-- ST_NearestValue
+-----------------------------------------------------------------------
+
+CREATE OR REPLACE FUNCTION st_nearestvalue(
+ rast raster, band integer,
+ pt geometry,
+ exclude_nodata_value boolean DEFAULT TRUE
+)
+ RETURNS double precision
+ AS 'MODULE_PATHNAME', 'RASTER_nearestValue'
+ LANGUAGE 'C' IMMUTABLE STRICT;
+
+CREATE OR REPLACE FUNCTION st_nearestvalue(
+ rast raster,
+ pt geometry,
+ exclude_nodata_value boolean DEFAULT TRUE
+)
+ RETURNS double precision
+ AS $$ SELECT st_nearestvalue($1, 1, $2, $3) $$
+ LANGUAGE 'SQL' IMMUTABLE STRICT;
+
+CREATE OR REPLACE FUNCTION st_nearestvalue(
+ rast raster, band integer,
+ x integer, y integer,
+ exclude_nodata_value boolean DEFAULT TRUE
+)
+ RETURNS double precision
+ AS $$ SELECT st_nearestvalue($1, $2, st_makepoint(st_raster2worldcoordx($1, $3, $4), st_raster2worldcoordy($1, $3, $4)), $5) $$
+ LANGUAGE 'SQL' IMMUTABLE STRICT;
+
+CREATE OR REPLACE FUNCTION st_nearestvalue(
+ rast raster,
+ x integer, y integer,
+ exclude_nodata_value boolean DEFAULT TRUE
+)
+ RETURNS double precision
+ AS $$ SELECT st_nearestvalue($1, 1, st_makepoint(st_raster2worldcoordx($1, $2, $3), st_raster2worldcoordy($1, $2, $3)), $4) $$
+ LANGUAGE 'SQL' IMMUTABLE STRICT;
+
+-----------------------------------------------------------------------
+-- ST_Neighborhood
+-----------------------------------------------------------------------
+
+CREATE OR REPLACE FUNCTION st_neighborhood(
+ rast raster, band integer,
+ ix integer, iy integer,
+ distance integer,
+ exclude_nodata_value boolean DEFAULT TRUE,
+ OUT x integer, OUT y integer,
+ OUT val double precision
+)
+ RETURNS SETOF record
+ AS 'MODULE_PATHNAME', 'RASTER_neighborhood'
+ LANGUAGE 'C' IMMUTABLE STRICT;
+
+CREATE OR REPLACE FUNCTION st_neighborhood(
+ rast raster,
+ ix integer, iy integer,
+ distance integer,
+ exclude_nodata_value boolean DEFAULT TRUE,
+ OUT x integer, OUT y integer,
+ OUT val double precision
+)
+ RETURNS SETOF record
+ AS $$ SELECT x, y, val FROM st_neighborhood($1, 1, $2, $3, $4, $5) $$
+ LANGUAGE 'SQL' IMMUTABLE STRICT;
+
+CREATE OR REPLACE FUNCTION st_neighborhood(
+ rast raster, band integer,
+ pt geometry,
+ distance integer,
+ exclude_nodata_value boolean DEFAULT TRUE,
+ OUT x integer, OUT y integer,
+ OUT val double precision
+)
+ RETURNS SETOF record
+ AS $$
+ DECLARE
+ wx int;
+ wy int;
+ BEGIN
+ IF (st_geometrytype($3) != 'ST_Point') THEN
+ RAISE EXCEPTION 'Attempting to get the neighbor of a pixel with a non-point geometry';
+ END IF;
+ wx := st_x($3);
+ wy := st_y($3);
+
+ RETURN QUERY
+ SELECT x, y, val
+ FROM st_neighborhood(
+ $1, $2,
+ st_world2rastercoordx(rast, wx, wy),
+ st_world2rastercoordy(rast, wx, wy),
+ $4,
+ $5
+ );
+ END;
+ $$ LANGUAGE 'plpgsql' IMMUTABLE STRICT;
+
+CREATE OR REPLACE FUNCTION st_neighborhood(
+ rast raster,
+ pt geometry,
+ distance integer,
+ exclude_nodata_value boolean DEFAULT TRUE,
+ OUT x integer, OUT y integer,
+ OUT val double precision
+)
+ RETURNS SETOF record
+ AS $$ SELECT x, y, val FROM st_neighborhood($1, 1, $2, $3, $4) $$
+ LANGUAGE 'SQL' IMMUTABLE STRICT;
+
------------------------------------------------------------------------------
-- raster constraint functions
-------------------------------------------------------------------------------
raster = rt_raster_new(xmax, ymax);
assert(raster);
- band = addBand(raster, PT_32BUI, 0, 0);
+ band = addBand(raster, PT_32BUI, 1, 0);
CHECK(band);
- rt_band_set_nodata(band, 0);
for (x = 0; x < xmax; x++) {
for (y = 0; y < ymax; y++) {
deepRelease(raster);
}
+static void testNearestPixel() {
+ rt_raster rast;
+ rt_band band;
+ uint32_t x, y;
+ int rtn;
+ const int maxX = 10;
+ const int maxY = 10;
+ rt_pixel npixels = NULL;
+
+ rast = rt_raster_new(maxX, maxY);
+ assert(rast);
+
+ band = addBand(rast, PT_32BUI, 1, 0);
+ CHECK(band);
+
+ for (x = 0; x < maxX; x++) {
+ for (y = 0; y < maxY; y++) {
+ rtn = rt_band_set_pixel(band, x, y, 1);
+ CHECK((rtn != -1));
+ }
+ }
+
+ rt_band_set_pixel(band, 0, 0, 0);
+ rt_band_set_pixel(band, 3, 0, 0);
+ rt_band_set_pixel(band, 6, 0, 0);
+ rt_band_set_pixel(band, 9, 0, 0);
+ rt_band_set_pixel(band, 1, 2, 0);
+ rt_band_set_pixel(band, 4, 2, 0);
+ rt_band_set_pixel(band, 7, 2, 0);
+ rt_band_set_pixel(band, 2, 4, 0);
+ rt_band_set_pixel(band, 5, 4, 0);
+ rt_band_set_pixel(band, 8, 4, 0);
+ rt_band_set_pixel(band, 0, 6, 0);
+ rt_band_set_pixel(band, 3, 6, 0);
+ rt_band_set_pixel(band, 6, 6, 0);
+ rt_band_set_pixel(band, 9, 6, 0);
+ rt_band_set_pixel(band, 1, 8, 0);
+ rt_band_set_pixel(band, 4, 8, 0);
+ rt_band_set_pixel(band, 7, 8, 0);
+
+ /* 0,0 */
+ rtn = rt_band_get_nearest_pixel(
+ band,
+ 0, 0,
+ 0,
+ 1,
+ &npixels
+ );
+ CHECK((rtn == 3));
+ if (rtn)
+ rtdealloc(npixels);
+
+ /* 1,1 */
+ rtn = rt_band_get_nearest_pixel(
+ band,
+ 1, 1,
+ 0,
+ 1,
+ &npixels
+ );
+ CHECK((rtn == 6));
+ if (rtn)
+ rtdealloc(npixels);
+
+ /* 4,4 */
+ rtn = rt_band_get_nearest_pixel(
+ band,
+ 4, 4,
+ 0,
+ 1,
+ &npixels
+ );
+ CHECK((rtn == 7));
+ if (rtn)
+ rtdealloc(npixels);
+
+ /* 4,4 distance 2 */
+ rtn = rt_band_get_nearest_pixel(
+ band,
+ 4, 4,
+ 2,
+ 1,
+ &npixels
+ );
+ CHECK((rtn == 19));
+ if (rtn)
+ rtdealloc(npixels);
+
+ /* 10,10 */
+ rtn = rt_band_get_nearest_pixel(
+ band,
+ 10, 10,
+ 0,
+ 1,
+ &npixels
+ );
+ CHECK((rtn == 1));
+ if (rtn)
+ rtdealloc(npixels);
+
+ /* 11,11 distance 1 */
+ rtn = rt_band_get_nearest_pixel(
+ band,
+ 11, 11,
+ 1,
+ 1,
+ &npixels
+ );
+ CHECK((rtn == 0));
+ if (rtn)
+ rtdealloc(npixels);
+
+ /* -1,-1 */
+ rtn = rt_band_get_nearest_pixel(
+ band,
+ -1, -1,
+ 0,
+ 1,
+ &npixels
+ );
+ CHECK((rtn == 3));
+ if (rtn)
+ rtdealloc(npixels);
+
+ /* -1,-1 distance 1 */
+ rtn = rt_band_get_nearest_pixel(
+ band,
+ -1, -1,
+ 1,
+ 1,
+ &npixels
+ );
+ CHECK((rtn == 0));
+ if (rtn)
+ rtdealloc(npixels);
+
+ /* -1,1 distance 1 */
+ rtn = rt_band_get_nearest_pixel(
+ band,
+ -1, 1,
+ 1,
+ 1,
+ &npixels
+ );
+ CHECK((rtn == 2));
+ if (rtn)
+ rtdealloc(npixels);
+
+ /* -2,2 distance 1 */
+ rtn = rt_band_get_nearest_pixel(
+ band,
+ -2, 2,
+ 1,
+ 1,
+ &npixels
+ );
+ CHECK((rtn == 0));
+ if (rtn)
+ rtdealloc(npixels);
+
+ /* -10,2 distance 3 */
+ rtn = rt_band_get_nearest_pixel(
+ band,
+ -10, 2,
+ 3,
+ 1,
+ &npixels
+ );
+ CHECK((rtn == 0));
+ if (rtn)
+ rtdealloc(npixels);
+
+ /* -10,2 distance 3 include NODATA */
+ rtn = rt_band_get_nearest_pixel(
+ band,
+ -10, 2,
+ 3,
+ 0,
+ &npixels
+ );
+ CHECK((rtn == 48));
+ if (rtn)
+ rtdealloc(npixels);
+
+ deepRelease(rast);
+}
+
int
main()
{
testCellGeoPoint();
printf("OK\n");
+ printf("Testing rt_band_get_nearest_pixel... ");
+ testNearestPixel();
+ printf("OK\n");
+
deepRelease(raster);
return EXIT_SUCCESS;
rt_valuepercent \
rt_bandmetadata \
rt_pixelvalue \
- drop_rt_band_properties_test
+ drop_rt_band_properties_test \
+ rt_neighborhood \
+ rt_nearestvalue
TEST_UTILITY = \
rt_utility \
--- /dev/null
+DROP TABLE IF EXISTS raster_nearestvalue;
+CREATE TABLE raster_nearestvalue (
+ rast raster
+);
+CREATE OR REPLACE FUNCTION make_test_raster()
+ RETURNS void
+ AS $$
+ DECLARE
+ width int := 10;
+ height int := 10;
+ x int;
+ y int;
+ rast raster;
+ BEGIN
+ rast := ST_MakeEmptyRaster(width, height, 0, 0, 1, -1, 0, 0, 0);
+ rast := ST_AddBand(rast, 1, '8BUI', 1, 0);
+
+ FOR y IN 1..height LOOP
+ FOR x IN 1..width LOOP
+ rast := ST_SetValue(rast, x, y, 2 * x + (1/3) * y);
+ END LOOP;
+ END LOOP;
+
+ rast := ST_SetValue(rast, 1, 1, 0);
+ rast := ST_SetValue(rast, 4, 1, 0);
+ rast := ST_SetValue(rast, 7, 1, 0);
+ rast := ST_SetValue(rast, 10, 1, 0);
+ rast := ST_SetValue(rast, 2, 3, 0);
+ rast := ST_SetValue(rast, 5, 3, 0);
+ rast := ST_SetValue(rast, 8, 3, 0);
+ rast := ST_SetValue(rast, 3, 5, 0);
+ rast := ST_SetValue(rast, 6, 5, 0);
+ rast := ST_SetValue(rast, 9, 5, 0);
+ rast := ST_SetValue(rast, 1, 7, 0);
+ rast := ST_SetValue(rast, 4, 7, 0);
+ rast := ST_SetValue(rast, 7, 7, 0);
+ rast := ST_SetValue(rast, 10, 7, 0);
+ rast := ST_SetValue(rast, 2, 9, 0);
+ rast := ST_SetValue(rast, 5, 9, 0);
+ rast := ST_SetValue(rast, 8, 9, 0);
+
+ INSERT INTO raster_nearestvalue VALUES (rast);
+
+ RETURN;
+ END;
+ $$ LANGUAGE 'plpgsql';
+SELECT make_test_raster();
+DROP FUNCTION IF EXISTS make_test_raster();
+
+SELECT
+ ST_NearestValue(rast, 1, 1, 1)
+FROM raster_nearestvalue;
+SELECT
+ ST_NearestValue(rast, 1, 2, 2)
+FROM raster_nearestvalue;
+SELECT
+ ST_NearestValue(rast, 1, 5, 5)
+FROM raster_nearestvalue;
+SELECT
+ ST_NearestValue(rast, 1, 5, 5)
+FROM raster_nearestvalue;
+SELECT
+ ST_NearestValue(rast, 1, 11, 11)
+FROM raster_nearestvalue;
+SELECT
+ ST_NearestValue(rast, 1, 12, 12)
+FROM raster_nearestvalue;
+SELECT
+ ST_NearestValue(rast, 1, 0, 0)
+FROM raster_nearestvalue;
+SELECT
+ ST_NearestValue(rast, 1, 0, 2)
+FROM raster_nearestvalue;
+SELECT
+ ST_NearestValue(rast, 1, -1, 3)
+FROM raster_nearestvalue;
+SELECT
+ ST_NearestValue(rast, 1, -9, 3)
+FROM raster_nearestvalue;
+SELECT
+ ST_NearestValue(rast, 1, -9, 3, FALSE)
+FROM raster_nearestvalue;
+
+SELECT
+ ST_NearestValue(rast, 1, ST_MakePoint(1, 1))
+FROM raster_nearestvalue;
+SELECT
+ ST_NearestValue(rast, 1, ST_MakePoint(2, 2))
+FROM raster_nearestvalue;
+SELECT
+ ST_NearestValue(rast, 1, ST_MakePoint(5, 5))
+FROM raster_nearestvalue;
+SELECT
+ ST_NearestValue(rast, 1, ST_MakePoint(5, 5))
+FROM raster_nearestvalue;
+SELECT
+ ST_NearestValue(rast, 1, ST_MakePoint(11, 11))
+FROM raster_nearestvalue;
+SELECT
+ ST_NearestValue(rast, 1, ST_MakePoint(12, 12))
+FROM raster_nearestvalue;
+SELECT
+ ST_NearestValue(rast, 1, ST_MakePoint(0, 0))
+FROM raster_nearestvalue;
+SELECT
+ ST_NearestValue(rast, 1, ST_MakePoint(0, 2))
+FROM raster_nearestvalue;
+SELECT
+ ST_NearestValue(rast, 1, ST_MakePoint(-1, 3))
+FROM raster_nearestvalue;
+SELECT
+ ST_NearestValue(rast, 1, ST_MakePoint(-9, 3))
+FROM raster_nearestvalue;
+SELECT
+ ST_NearestValue(rast, 1, ST_MakePoint(-9, 3), FALSE)
+FROM raster_nearestvalue;
+
+DROP TABLE IF EXISTS raster_nearestvalue;
--- /dev/null
+NOTICE: table "raster_nearestvalue" does not exist, skipping
+4
+4
+10
+10
+20
+20
+4
+2
+2
+2
+0
+4
+4
+10
+10
+18
+18
+4
+4
+4
+4
+0
--- /dev/null
+DROP TABLE IF EXISTS raster_neighborhood;
+CREATE TABLE raster_neighborhood (
+ rast raster
+);
+CREATE OR REPLACE FUNCTION make_test_raster()
+ RETURNS void
+ AS $$
+ DECLARE
+ rast raster;
+ BEGIN
+ rast := ST_MakeEmptyRaster(10, 10, 0, 0, 1, -1, 0, 0, 0);
+ rast := ST_AddBand(rast, 1, '8BUI', 1, 0);
+
+ rast := ST_SetValue(rast, 1, 1, 0);
+ rast := ST_SetValue(rast, 4, 1, 0);
+ rast := ST_SetValue(rast, 7, 1, 0);
+ rast := ST_SetValue(rast, 10, 1, 0);
+ rast := ST_SetValue(rast, 2, 3, 0);
+ rast := ST_SetValue(rast, 5, 3, 0);
+ rast := ST_SetValue(rast, 8, 3, 0);
+ rast := ST_SetValue(rast, 3, 5, 0);
+ rast := ST_SetValue(rast, 6, 5, 0);
+ rast := ST_SetValue(rast, 9, 5, 0);
+ rast := ST_SetValue(rast, 1, 7, 0);
+ rast := ST_SetValue(rast, 4, 7, 0);
+ rast := ST_SetValue(rast, 7, 7, 0);
+ rast := ST_SetValue(rast, 10, 7, 0);
+ rast := ST_SetValue(rast, 2, 9, 0);
+ rast := ST_SetValue(rast, 5, 9, 0);
+ rast := ST_SetValue(rast, 8, 9, 0);
+
+ INSERT INTO raster_neighborhood VALUES (rast);
+
+ RETURN;
+ END;
+ $$ LANGUAGE 'plpgsql';
+SELECT make_test_raster();
+DROP FUNCTION IF EXISTS make_test_raster();
+
+SELECT
+ ST_Neighborhood(rast, 1, 1, 1, 1)
+FROM raster_neighborhood;
+SELECT
+ ST_Neighborhood(rast, 1, 2, 2, 1)
+FROM raster_neighborhood;
+SELECT
+ ST_Neighborhood(rast, 1, 5, 5, 1)
+FROM raster_neighborhood;
+SELECT
+ ST_Neighborhood(rast, 1, 5, 5, 2)
+FROM raster_neighborhood;
+SELECT
+ ST_Neighborhood(rast, 1, 11, 11, 1)
+FROM raster_neighborhood;
+SELECT
+ ST_Neighborhood(rast, 1, 12, 12, 1)
+FROM raster_neighborhood;
+SELECT
+ ST_Neighborhood(rast, 1, 0, 0, 1)
+FROM raster_neighborhood;
+SELECT
+ ST_Neighborhood(rast, 1, 0, 2, 1)
+FROM raster_neighborhood;
+SELECT
+ ST_Neighborhood(rast, 1, -1, 3, 1)
+FROM raster_neighborhood;
+SELECT
+ ST_Neighborhood(rast, 1, -9, 3, 3)
+FROM raster_neighborhood;
+SELECT
+ ST_Neighborhood(rast, 1, -9, 3, 3, FALSE)
+FROM raster_neighborhood;
+
+DROP TABLE IF EXISTS raster_neighborhood;
--- /dev/null
+NOTICE: table "raster_neighborhood" does not exist, skipping
+(1,2,1)
+(2,2,1)
+(2,1,1)
+(2,1,1)
+(3,1,1)
+(1,3,1)
+(3,3,1)
+(1,2,1)
+(3,2,1)
+(4,4,1)
+(5,4,1)
+(6,4,1)
+(4,6,1)
+(5,6,1)
+(6,6,1)
+(4,5,1)
+(4,4,1)
+(5,4,1)
+(6,4,1)
+(4,6,1)
+(5,6,1)
+(6,6,1)
+(4,5,1)
+(3,3,1)
+(4,3,1)
+(6,3,1)
+(7,3,1)
+(3,7,1)
+(5,7,1)
+(6,7,1)
+(3,4,1)
+(3,6,1)
+(7,4,1)
+(7,5,1)
+(7,6,1)
+(10,10,1)
+NOTICE: Pixel has no neighbors for band at distance 1
+NOTICE: Pixel has no neighbors for band at distance 1
+(1,3,1)
+(1,2,1)
+NOTICE: Pixel has no neighbors for band at distance 1
+NOTICE: Pixel has no neighbors for band at distance 3
+(-10,2,0)
+(-9,2,0)
+(-8,2,0)
+(-10,4,0)
+(-9,4,0)
+(-8,4,0)
+(-10,3,0)
+(-8,3,0)
+(-11,1,0)
+(-10,1,0)
+(-9,1,0)
+(-8,1,0)
+(-7,1,0)
+(-11,5,0)
+(-10,5,0)
+(-9,5,0)
+(-8,5,0)
+(-7,5,0)
+(-11,2,0)
+(-11,3,0)
+(-11,4,0)
+(-7,2,0)
+(-7,3,0)
+(-7,4,0)
+(-12,0,0)
+(-11,0,0)
+(-10,0,0)
+(-9,0,0)
+(-8,0,0)
+(-7,0,0)
+(-6,0,0)
+(-12,6,0)
+(-11,6,0)
+(-10,6,0)
+(-9,6,0)
+(-8,6,0)
+(-7,6,0)
+(-6,6,0)
+(-12,1,0)
+(-12,2,0)
+(-12,3,0)
+(-12,4,0)
+(-12,5,0)
+(-6,1,0)
+(-6,2,0)
+(-6,3,0)
+(-6,4,0)
+(-6,5,0)
projects / postgis / commitdiff