RASTER_nearestValue()

Datum RASTER_nearestValue

(

PG_FUNCTION_ARGS

)

Definition at line 2044 of file rtpg_pixel.c.

{
        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;
        int isnodata = 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) {
                PG_FREE_IF_COPY(pgraster, 0);
                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_FREE_IF_COPY(pgraster, 0);
                PG_RETURN_NULL();
        }
 
        /* point */
        geom = PG_GETARG_GSERIALIZED_P(2);
        if (gserialized_get_type(geom) != POINTTYPE) {
                elog(NOTICE, "Geometry provided must be a point");
                rt_raster_destroy(raster);
                PG_FREE_IF_COPY(pgraster, 0);
                PG_FREE_IF_COPY(geom, 2);
                PG_RETURN_NULL();
        }
 
        /* exclude_nodata_value flag */
        if (!PG_ARGISNULL(3))
                exclude_nodata_value = PG_GETARG_BOOL(3);
 
        /* SRIDs of raster and geometry must match  */
        if (clamp_srid(gserialized_get_srid(geom)) != clamp_srid(rt_raster_get_srid(raster))) {
                elog(NOTICE, "SRIDs of geometry and raster do not match");
                rt_raster_destroy(raster);
                PG_FREE_IF_COPY(pgraster, 0);
                PG_FREE_IF_COPY(geom, 2);
                PG_RETURN_NULL();
        }
 
        /* 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_FREE_IF_COPY(pgraster, 0);
                PG_FREE_IF_COPY(geom, 2);
                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_FREE_IF_COPY(pgraster, 0);
                PG_FREE_IF_COPY(geom, 2);
                PG_RETURN_NULL();
        }
 
        /* Get a 2D version of the geometry if necessary */
        if (lwgeom_ndims(lwgeom) > 2) {
                LWGEOM *lwgeom2d = lwgeom_force_2d(lwgeom);
                lwgeom_free(lwgeom);
                lwgeom = lwgeom2d;
        }
 
        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
        ) != ES_NONE) {
                rt_raster_destroy(raster);
                PG_FREE_IF_COPY(pgraster, 0);
                lwgeom_free(lwgeom);
                PG_FREE_IF_COPY(geom, 2);
                elog(ERROR, "RASTER_nearestValue: Could not compute pixel coordinates from spatial coordinates");
                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, &isnodata) != ES_NONE) {
                        rt_raster_destroy(raster);
                        PG_FREE_IF_COPY(pgraster, 0);
                        lwgeom_free(lwgeom);
                        PG_FREE_IF_COPY(geom, 2);
                        elog(ERROR, "RASTER_nearestValue: Could not get pixel value for band at index %d", bandindex);
                        PG_RETURN_NULL();
                }
 
                /* value at point, return value */
                if (!exclude_nodata_value || !isnodata) {
                        rt_raster_destroy(raster);
                        PG_FREE_IF_COPY(pgraster, 0);
                        lwgeom_free(lwgeom);
                        PG_FREE_IF_COPY(geom, 2);
 
                        PG_RETURN_FLOAT8(value);
                }
        }
 
        /* get neighborhood */
        count = rt_band_get_nearest_pixel(
                band,
                x, y,
                0, 0,
                exclude_nodata_value,
                &npixels
        );
        rt_band_destroy(band);
        /* error or no neighbors */
        if (count < 1) {
                /* error */
                if (count < 0)
                        elog(NOTICE, "Could not 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);
                PG_FREE_IF_COPY(geom, 2);
                rt_raster_destroy(raster);
                PG_FREE_IF_COPY(pgraster, 0);
                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) {
                                lwgeom_free(lwgeom);
                                PG_FREE_IF_COPY(geom, 2);
                                rt_raster_destroy(raster);
                                PG_FREE_IF_COPY(pgraster, 0);
                                elog(ERROR, "RASTER_nearestValue: Could not get polygon of neighboring pixel");
                                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);
        PG_FREE_IF_COPY(geom, 2);
        rt_raster_destroy(raster);
        PG_FREE_IF_COPY(pgraster, 0);
 
        if (hasvalue)
                PG_RETURN_FLOAT8(value);
        else
                PG_RETURN_NULL();
}

References clamp_srid(), 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, 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, POINT2D::x, and POINT2D::y.

Here is the call graph for this function: