RASTER_neighborhood()

Datum RASTER_neighborhood

(

PG_FUNCTION_ARGS

)

Definition at line 2050 of file rtpg_pixel.c.

References ovdump::band, genraster::count, distance(), ES_NONE, FALSE, rt_pixel_t::nodata, pixval::pixval, POSTGIS_RT_DEBUGF, rtrowdump::raster, rt_band_destroy(), rt_band_get_hasnodata_flag(), rt_band_get_height(), rt_band_get_min_value(), rt_band_get_nearest_pixel(), rt_band_get_nodata(), rt_band_get_pixel(), rt_band_get_width(), rt_pixel_set_to_array(), rt_raster_deserialize(), rt_raster_destroy(), rt_raster_get_band(), rt_raster_get_num_bands(), TRUE, rt_pixel_t::value, pixval::x, rt_pixel_t::x, pixval::y, and rt_pixel_t::y.

Referenced by RASTER_nearestValue().

{
        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 _x = 0;
        int _y = 0;
        int distance[2] = {0};
        bool exclude_nodata_value = TRUE;
        double pixval;
        int isnodata = 0;

        rt_pixel npixels = NULL;
        int count;
        double **value2D = NULL;
        int **nodata2D = NULL;

        int i = 0;
        int j = 0;
        int k = 0;
        Datum *value1D = NULL;
        bool *nodata1D = NULL;
        int dim[2] = {0};
        int lbound[2] = {1, 1};
        ArrayType *mdArray = NULL;

        int16 typlen;
        bool typbyval;
        char typalign;

        /* pgraster is null, return nothing */
        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_neighborhood: 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();
        }

        /* pixel column, 1-based */
        x = PG_GETARG_INT32(2);
        _x = x - 1;

        /* pixel row, 1-based */
        y = PG_GETARG_INT32(3);
        _y = y - 1;

        /* distance X axis */
        distance[0] = PG_GETARG_INT32(4);
        if (distance[0] < 0) {
                elog(NOTICE, "Invalid value for distancex (must be >= zero). Returning NULL");
                rt_raster_destroy(raster);
                PG_FREE_IF_COPY(pgraster, 0);
                PG_RETURN_NULL();
        }
        distance[0] = (uint16_t) distance[0];

        /* distance Y axis */
        distance[1] = PG_GETARG_INT32(5);
        if (distance[1] < 0) {
                elog(NOTICE, "Invalid value for distancey (must be >= zero). Returning NULL");
                rt_raster_destroy(raster);
                PG_FREE_IF_COPY(pgraster, 0);
                PG_RETURN_NULL();
        }
        distance[1] = (uint16_t) distance[1];

        /* exclude_nodata_value flag */
        if (!PG_ARGISNULL(6))
                exclude_nodata_value = PG_GETARG_BOOL(6);

        /* 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_RETURN_NULL();
        }

        /* get neighborhood */
        count = 0;
        npixels = NULL;
        if (distance[0] > 0 || distance[1] > 0) {
                count = rt_band_get_nearest_pixel(
                        band,
                        _x, _y,
                        distance[0], distance[1],
                        exclude_nodata_value,
                        &npixels
                );
                /* error */
                if (count < 0) {
                        elog(NOTICE, "Could not get the pixel's neighborhood for band at index %d", bandindex);

                        rt_band_destroy(band);
                        rt_raster_destroy(raster);
                        PG_FREE_IF_COPY(pgraster, 0);

                        PG_RETURN_NULL();
                }
        }

        /* get pixel's value */
        if (
                (_x >= 0 && _x < rt_band_get_width(band)) &&
                (_y >= 0 && _y < rt_band_get_height(band))
        ) {
                if (rt_band_get_pixel(
                        band,
                        _x, _y,
                        &pixval,
                        &isnodata
                ) != ES_NONE) {
                        elog(NOTICE, "Could not get the pixel of band at index %d. Returning NULL", bandindex);
                        rt_band_destroy(band);
                        rt_raster_destroy(raster);
                        PG_FREE_IF_COPY(pgraster, 0);
                        PG_RETURN_NULL();
                }
        }
        /* outside band extent, set to NODATA */
        else {
                /* has NODATA, use NODATA */
                if (rt_band_get_hasnodata_flag(band))
                        rt_band_get_nodata(band, &pixval);
                /* no NODATA, use min possible value */
                else
                        pixval = rt_band_get_min_value(band);
                isnodata = 1;
        }
        POSTGIS_RT_DEBUGF(4, "pixval: %f", pixval);


        /* add pixel to neighborhood */
        count++;
        if (count > 1)
                npixels = (rt_pixel) repalloc(npixels, sizeof(struct rt_pixel_t) * count);
        else
                npixels = (rt_pixel) palloc(sizeof(struct rt_pixel_t));
        if (npixels == NULL) {

                rt_band_destroy(band);
                rt_raster_destroy(raster);
                PG_FREE_IF_COPY(pgraster, 0);

                elog(ERROR, "RASTER_neighborhood: Could not reallocate memory for neighborhood");
                PG_RETURN_NULL();
        }
        npixels[count - 1].x = _x;
        npixels[count - 1].y = _y;
        npixels[count - 1].nodata = 1;
        npixels[count - 1].value = pixval;

        /* set NODATA */
        if (!exclude_nodata_value || !isnodata) {
                npixels[count - 1].nodata = 0;
        }

        /* free unnecessary stuff */
        rt_band_destroy(band);
        rt_raster_destroy(raster);
        PG_FREE_IF_COPY(pgraster, 0);

        /* convert set of rt_pixel to 2D array */
        /* dim is passed with element 0 being Y-axis and element 1 being X-axis */
        count = rt_pixel_set_to_array(
                npixels, count, NULL,
                _x, _y,
                distance[0], distance[1],
                &value2D,
                &nodata2D,
                &(dim[1]), &(dim[0])
        );
        pfree(npixels);
        if (count != ES_NONE) {
                elog(NOTICE, "Could not create 2D array of neighborhood");
                PG_RETURN_NULL();
        }

        /* 1D arrays for values and nodata from 2D arrays */
        value1D = palloc(sizeof(Datum) * dim[0] * dim[1]);
        nodata1D = palloc(sizeof(bool) * dim[0] * dim[1]);

        if (value1D == NULL || nodata1D == NULL) {

                for (i = 0; i < dim[0]; i++) {
                        pfree(value2D[i]);
                        pfree(nodata2D[i]);
                }
                pfree(value2D);
                pfree(nodata2D);

                elog(ERROR, "RASTER_neighborhood: Could not allocate memory for return 2D array");
                PG_RETURN_NULL();
        }

        /* copy values from 2D arrays to 1D arrays */
        k = 0;
        /* Y-axis */
        for (i = 0; i < dim[0]; i++) {
                /* X-axis */
                for (j = 0; j < dim[1]; j++) {
                        nodata1D[k] = (bool) nodata2D[i][j];
                        if (!nodata1D[k])
                                value1D[k] = Float8GetDatum(value2D[i][j]);
                        else
                                value1D[k] = PointerGetDatum(NULL);

                        k++;
                }
        }

        /* no more need for 2D arrays */
        for (i = 0; i < dim[0]; i++) {
                pfree(value2D[i]);
                pfree(nodata2D[i]);
        }
        pfree(value2D);
        pfree(nodata2D);

        /* info about the type of item in the multi-dimensional array (float8). */
        get_typlenbyvalalign(FLOAT8OID, &typlen, &typbyval, &typalign);

        mdArray = construct_md_array(
                value1D, nodata1D,
                2, dim, lbound,
                FLOAT8OID,
                typlen, typbyval, typalign
        );

        pfree(value1D);
        pfree(nodata1D);

        PG_RETURN_ARRAYTYPE_P(mdArray);
}

Here is the call graph for this function:

Here is the caller graph for this function: