RASTER_clip()

Datum RASTER_clip

(

PG_FUNCTION_ARGS

)

Definition at line 2967 of file rtpg_mapalgebra.c.

{
        rt_pgraster *pgraster = NULL;
        LWGEOM *rastgeom = NULL;
        double gt[6] = {0};
        int32_t srid = SRID_UNKNOWN;
 
        rt_pgraster *pgrtn = NULL;
        rt_raster rtn = NULL;
 
        GSERIALIZED *gser = NULL;
        LWGEOM *geom = NULL;
        lwvarlena_t *wkb = NULL;
 
        ArrayType *array;
        Oid etype;
        Datum *e;
        bool *nulls;
 
        int16 typlen;
        bool typbyval;
        char typalign;
 
        int i = 0;
        int j = 0;
        int k = 0;
        rtpg_clip_arg arg = NULL;
        LWGEOM *tmpgeom = NULL;
 
        rt_pixtype pixtype;
        int hasnodata;
        double nodataval;
 
        double offset[4] = {0.};
        int input_x = 0;
        int input_y = 0;
        int mask_x = 0;
        int mask_y = 0;
        int x = 0;
        int y = 0;
        int width = 0;
        int height = 0;
        int mask_width = 0;
        int mask_height = 0;
        rt_band input_band = NULL;
        rt_band mask_band = NULL;
        rt_band output_band = NULL;
        double value;
        int isnodata;
 
        POSTGIS_RT_DEBUG(3, "Starting...");
 
        /* raster or geometry is NULL, return NULL */
        if (PG_ARGISNULL(0) || PG_ARGISNULL(2))
                PG_RETURN_NULL();
 
        /* init arg */
        arg = rtpg_clip_arg_init();
        if (arg == NULL) {
                elog(ERROR, "RASTER_clip: Could not initialize argument structure");
                PG_RETURN_NULL();
        }
 
        /* raster (0) */
        pgraster = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
 
        /* Get raster object */
        arg->raster = rt_raster_deserialize(pgraster, FALSE);
        if (arg->raster == NULL) {
                rtpg_clip_arg_destroy(arg);
                PG_FREE_IF_COPY(pgraster, 0);
                elog(ERROR, "RASTER_clip: Could not deserialize raster");
                PG_RETURN_NULL();
        }
 
        /* raster is empty, return empty raster */
        if (rt_raster_is_empty(arg->raster) || rt_raster_get_num_bands(arg->raster) == 0) {
                elog(NOTICE, "Input raster is empty or has no bands. Returning empty raster");
 
                rtpg_clip_arg_destroy(arg);
                PG_FREE_IF_COPY(pgraster, 0);
 
                rtn = rt_raster_new(0, 0);
                if (rtn == NULL) {
                        elog(ERROR, "RASTER_clip: Could not create empty raster");
                        PG_RETURN_NULL();
                }
 
                pgrtn = rt_raster_serialize(rtn);
                rt_raster_destroy(rtn);
                if (NULL == pgrtn)
                        PG_RETURN_NULL();
 
                SET_VARSIZE(pgrtn, pgrtn->size);
                PG_RETURN_POINTER(pgrtn);
        }
 
        /* metadata */
        rt_raster_get_geotransform_matrix(arg->raster, gt);
        srid = clamp_srid(rt_raster_get_srid(arg->raster));
 
        /* geometry (2) */
        gser = PG_GETARG_GSERIALIZED_P(2);
        geom = lwgeom_from_gserialized(gser);
 
        /* Get a 2D version of the geometry if necessary */
        if (lwgeom_ndims(geom) > 2) {
                LWGEOM *geom2d = lwgeom_force_2d(geom);
                lwgeom_free(geom);
                geom = geom2d;
        }
 
        /* check that SRIDs match */
        if (srid != clamp_srid(gserialized_get_srid(gser))) {
                elog(NOTICE, "Geometry provided does not have the same SRID as the raster. Returning NULL");
 
                rtpg_clip_arg_destroy(arg);
                PG_FREE_IF_COPY(pgraster, 0);
                lwgeom_free(geom);
                PG_FREE_IF_COPY(gser, 2);
 
                PG_RETURN_NULL();
        }
 
        /* crop (4) */
        if (!PG_ARGISNULL(4) && !PG_GETARG_BOOL(4))
                arg->extenttype = ET_FIRST;
 
        /* get intersection geometry of input raster and input geometry */
        if (rt_raster_get_convex_hull(arg->raster, &rastgeom) != ES_NONE) {
 
                rtpg_clip_arg_destroy(arg);
                PG_FREE_IF_COPY(pgraster, 0);
                lwgeom_free(geom);
                PG_FREE_IF_COPY(gser, 2);
 
                elog(ERROR, "RASTER_clip: Could not get convex hull of raster");
                PG_RETURN_NULL();
        }
 
        tmpgeom = lwgeom_intersection(rastgeom, geom);
        lwgeom_free(rastgeom);
        lwgeom_free(geom);
        PG_FREE_IF_COPY(gser, 2);
        geom = tmpgeom;
 
        /* intersection is empty AND extent type is INTERSECTION, return empty */
        if (lwgeom_is_empty(geom) && arg->extenttype == ET_INTERSECTION) {
                elog(NOTICE, "The input raster and input geometry do not intersect. Returning empty raster");
 
                rtpg_clip_arg_destroy(arg);
                PG_FREE_IF_COPY(pgraster, 0);
                lwgeom_free(geom);
 
                rtn = rt_raster_new(0, 0);
                if (rtn == NULL) {
                        elog(ERROR, "RASTER_clip: Could not create empty raster");
                        PG_RETURN_NULL();
                }
 
                pgrtn = rt_raster_serialize(rtn);
                rt_raster_destroy(rtn);
                if (NULL == pgrtn)
                        PG_RETURN_NULL();
 
                SET_VARSIZE(pgrtn, pgrtn->size);
                PG_RETURN_POINTER(pgrtn);
        }
 
        /* nband (1) */
        if (!PG_ARGISNULL(1)) {
                array = PG_GETARG_ARRAYTYPE_P(1);
                etype = ARR_ELEMTYPE(array);
                get_typlenbyvalalign(etype, &typlen, &typbyval, &typalign);
 
                switch (etype) {
                        case INT2OID:
                        case INT4OID:
                                break;
                        default:
                                rtpg_clip_arg_destroy(arg);
                                PG_FREE_IF_COPY(pgraster, 0);
                                lwgeom_free(geom);
                                elog(ERROR, "RASTER_clip: Invalid data type for band indexes");
                                PG_RETURN_NULL();
                                break;
                }
 
                deconstruct_array(
                        array, etype,
                        typlen, typbyval, typalign,
                        &e, &nulls, &(arg->numbands)
                );
 
                arg->band = palloc(sizeof(struct rtpg_clip_band_t) * arg->numbands);
                if (arg->band == NULL) {
                        rtpg_clip_arg_destroy(arg);
                        PG_FREE_IF_COPY(pgraster, 0);
                        lwgeom_free(geom);
                        elog(ERROR, "RASTER_clip: Could not allocate memory for band arguments");
                        PG_RETURN_NULL();
                }
 
                for (i = 0, j = 0; i < arg->numbands; i++) {
                        if (nulls[i]) continue;
 
                        switch (etype) {
                                case INT2OID:
                                        arg->band[j].nband = DatumGetInt16(e[i]) - 1;
                                        break;
                                case INT4OID:
                                        arg->band[j].nband = DatumGetInt32(e[i]) - 1;
                                        break;
                        }
 
                        j++;
                }
 
                if (j < arg->numbands) {
                        arg->band = repalloc(arg->band, sizeof(struct rtpg_clip_band_t) * j);
                        if (arg->band == NULL) {
                                rtpg_clip_arg_destroy(arg);
                                PG_FREE_IF_COPY(pgraster, 0);
                                lwgeom_free(geom);
                                elog(ERROR, "RASTER_clip: Could not reallocate memory for band arguments");
                                PG_RETURN_NULL();
                        }
 
                        arg->numbands = j;
                }
 
                /* validate band */
                for (i = 0; i < arg->numbands; i++) {
                        if (!rt_raster_has_band(arg->raster, arg->band[i].nband)) {
                                elog(NOTICE, "Band at index %d not found in raster", arg->band[i].nband + 1);
                                rtpg_clip_arg_destroy(arg);
                                PG_FREE_IF_COPY(pgraster, 0);
                                lwgeom_free(geom);
                                PG_RETURN_NULL();
                        }
 
                        arg->band[i].hasnodata = 0;
                        arg->band[i].nodataval = 0;
                }
        }
        else {
                arg->numbands = rt_raster_get_num_bands(arg->raster);
 
                /* raster may have no bands */
                if (arg->numbands) {
                        arg->band = palloc(sizeof(struct rtpg_clip_band_t) * arg->numbands);
                        if (arg->band == NULL) {
 
                                rtpg_clip_arg_destroy(arg);
                                PG_FREE_IF_COPY(pgraster, 0);
                                lwgeom_free(geom);
 
                                elog(ERROR, "RASTER_clip: Could not allocate memory for band arguments");
                                PG_RETURN_NULL();
                        }
 
                        for (i = 0; i < arg->numbands; i++) {
                                arg->band[i].nband = i;
                                arg->band[i].hasnodata = 0;
                                arg->band[i].nodataval = 0;
                        }
                }
        }
 
        /* nodataval (3) */
        if (!PG_ARGISNULL(3)) {
                array = PG_GETARG_ARRAYTYPE_P(3);
                etype = ARR_ELEMTYPE(array);
                get_typlenbyvalalign(etype, &typlen, &typbyval, &typalign);
 
                switch (etype) {
                        case FLOAT4OID:
                        case FLOAT8OID:
                                break;
                        default:
                                rtpg_clip_arg_destroy(arg);
                                PG_FREE_IF_COPY(pgraster, 0);
                                lwgeom_free(geom);
                                elog(ERROR, "RASTER_clip: Invalid data type for NODATA values");
                                PG_RETURN_NULL();
                                break;
                }
 
                deconstruct_array(
                        array, etype,
                        typlen, typbyval, typalign,
                        &e, &nulls, &k
                );
 
                /* it doesn't matter if there are more nodataval */
                for (i = 0, j = 0; i < arg->numbands; i++, j++) {
                        /* cap j to the last nodataval element */
                        if (j >= k)
                                j = k - 1;
 
                        if (nulls[j])
                                continue;
 
                        arg->band[i].hasnodata = 1;
                        switch (etype) {
                                case FLOAT4OID:
                                        arg->band[i].nodataval = DatumGetFloat4(e[j]);
                                        break;
                                case FLOAT8OID:
                                        arg->band[i].nodataval = DatumGetFloat8(e[j]);
                                        break;
                        }
                }
        }
 
        /* get wkb of geometry */
        POSTGIS_RT_DEBUG(3, "getting wkb of geometry");
        wkb = lwgeom_to_wkb_varlena(geom, WKB_SFSQL);
        lwgeom_free(geom);
 
        /* rasterize geometry */
        arg->mask = rt_raster_gdal_rasterize((unsigned char *)wkb->data,
                                             LWSIZE_GET(wkb->size) - LWVARHDRSZ,
                                             NULL,
                                             0,
                                             NULL,
                                             NULL,
                                             NULL,
                                             NULL,
                                             NULL,
                                             NULL,
                                             NULL,
                                             &(gt[1]),
                                             &(gt[5]),
                                             NULL,
                                             NULL,
                                             &(gt[0]),
                                             &(gt[3]),
                                             &(gt[2]),
                                             &(gt[4]),
                                             NULL);
 
        pfree(wkb);
        if (arg->mask == NULL) {
                rtpg_clip_arg_destroy(arg);
                PG_FREE_IF_COPY(pgraster, 0);
                elog(ERROR, "RASTER_clip: Could not rasterize intersection geometry");
                PG_RETURN_NULL();
        }
 
        /* set SRID */
        rt_raster_set_srid(arg->mask, srid);
 
        mask_width = rt_raster_get_width(arg->mask);
        mask_height = rt_raster_get_height(arg->mask);
 
        if (rt_raster_from_two_rasters(arg->raster, arg->mask, arg->extenttype, &rtn, offset) != ES_NONE) {
                rtpg_clip_arg_destroy(arg);
                PG_FREE_IF_COPY(pgraster, 0);
                elog(ERROR, "RASTER_clip: Could not compute extent of rasters");
                PG_RETURN_NULL();
        }
 
        width = rt_raster_get_width(rtn);
        height = rt_raster_get_height(rtn);
 
        mask_band = rt_raster_get_band(arg->mask, 0);
 
        for (i = 0; i < arg->numbands; i++) {
                input_band = rt_raster_get_band(arg->raster, arg->band[i].nband);
 
                /* band metadata */
                pixtype = rt_band_get_pixtype(input_band);
 
                if (arg->band[i].hasnodata) {
                        hasnodata = 1;
                        nodataval = arg->band[i].nodataval;
                }
                else if (rt_band_get_hasnodata_flag(input_band)) {
                        hasnodata = 1;
                        rt_band_get_nodata(input_band, &nodataval);
                }
                else {
                        hasnodata = 0;
                        nodataval = rt_band_get_min_value(input_band);
                }
 
                if (rt_raster_generate_new_band(rtn, pixtype, nodataval, hasnodata, nodataval, i) < 0) {
                        rtpg_clip_arg_destroy(arg);
                        PG_FREE_IF_COPY(pgraster, 0);
                        elog(ERROR, "RASTER_clip: Could not add new band in output raster");
                        PG_RETURN_NULL();
                }
 
                if (rt_band_get_isnodata_flag(input_band)) {
                        continue;
                }
 
                output_band = rt_raster_get_band(rtn, arg->band[i].nband);
 
                if (!mask_band) {
                        continue;
                }
 
                for (y = 0; y < height; y++) {
                        for (x = 0; x < width; x++) {
                                mask_x = x - (int)offset[2];
                                mask_y = y - (int)offset[3];
 
                                if (!(
                                        mask_x >= 0 &&
                                        mask_x < mask_width &&
                                        mask_y >= 0 &&
                                        mask_y < mask_height
                                )) {
                                        continue;
                                }
 
                                if (rt_band_get_pixel(mask_band, mask_x, mask_y, &value, &isnodata) != ES_NONE) {
                                        rtpg_clip_arg_destroy(arg);
                                        PG_FREE_IF_COPY(pgraster, 0);
                                        elog(ERROR, "RASTER_clip: Could not get pixel value");
                                        PG_RETURN_NULL();
                                }
 
                                if (isnodata) {
                                        continue;
                                }
 
                                input_x = x - (int)offset[0];
                                input_y = y - (int)offset[1];
 
                                if (rt_band_get_pixel(input_band, input_x, input_y, &value, &isnodata) != ES_NONE) {
                                        rtpg_clip_arg_destroy(arg);
                                        PG_FREE_IF_COPY(pgraster, 0);
                                        elog(ERROR, "RASTER_clip: Could not get pixel value");
                                        PG_RETURN_NULL();
                                }
 
                                if (isnodata) {
                                        continue;
                                }
 
                                if (rt_band_set_pixel(output_band, x, y, value, NULL)) {
                                        rtpg_clip_arg_destroy(arg);
                                        PG_FREE_IF_COPY(pgraster, 0);
                                        elog(ERROR, "RASTER_clip: Could not set pixel value");
                                        PG_RETURN_NULL();
                                }
                        }
                }
        }
 
        rtpg_clip_arg_destroy(arg);
        PG_FREE_IF_COPY(pgraster, 0);
 
        pgrtn = rt_raster_serialize(rtn);
        rt_raster_destroy(rtn);
 
        POSTGIS_RT_DEBUG(3, "Finished");
 
        if (!pgrtn)
                PG_RETURN_NULL();
 
        SET_VARSIZE(pgrtn, pgrtn->size);
        PG_RETURN_POINTER(pgrtn);
}

References rtpg_clip_arg_t::band, clamp_srid(), lwvarlena_t::data, ES_NONE, ET_FIRST, ET_INTERSECTION, rtpg_clip_arg_t::extenttype, FALSE, gserialized_get_srid(), window::gt, rtpg_clip_band_t::hasnodata, lwgeom_force_2d(), lwgeom_free(), lwgeom_from_gserialized(), lwgeom_intersection(), lwgeom_is_empty(), lwgeom_ndims(), lwgeom_to_wkb_varlena(), LWSIZE_GET, LWVARHDRSZ, rtpg_clip_arg_t::mask, rtpg_clip_band_t::nband, rtpg_clip_band_t::nodataval, rtpg_clip_arg_t::numbands, POSTGIS_RT_DEBUG, rtpg_clip_arg_t::raster, rt_band_get_hasnodata_flag(), rt_band_get_isnodata_flag(), rt_band_get_min_value(), rt_band_get_nodata(), rt_band_get_pixel(), rt_band_get_pixtype(), rt_band_set_pixel(), rt_raster_deserialize(), rt_raster_destroy(), rt_raster_from_two_rasters(), rt_raster_gdal_rasterize(), rt_raster_generate_new_band(), rt_raster_get_band(), rt_raster_get_convex_hull(), rt_raster_get_geotransform_matrix(), rt_raster_get_height(), rt_raster_get_num_bands(), rt_raster_get_srid(), rt_raster_get_width(), rt_raster_has_band(), rt_raster_is_empty(), rt_raster_new(), rt_raster_serialize(), rt_raster_set_srid(), rtpg_clip_arg_destroy(), rtpg_clip_arg_init(), lwvarlena_t::size, rt_raster_serialized_t::size, SRID_UNKNOWN, genraster::value, WKB_SFSQL, pixval::x, and pixval::y.

Here is the call graph for this function: