RASTER_clip()

Datum RASTER_clip

(

PG_FUNCTION_ARGS

)

Definition at line 3026 of file rtpg_mapalgebra.c.

{
        rt_pgraster *pgraster = NULL;
        LWGEOM *rastgeom = NULL;
        double gt[6] = {0};
        int srid = SRID_UNKNOWN;
 
        rt_pgraster *pgrtn = NULL;
        rt_raster rtn = NULL;
 
        GSERIALIZED *gser = NULL;
        LWGEOM *geom = NULL;
        unsigned char *wkb = NULL;
        size_t wkb_len;
 
        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_iterator itrset;
 
        rt_raster _raster = NULL;
        rt_band band = NULL;
        rt_pixtype pixtype;
        int hasnodata;
        double nodataval;
        int noerr = 0;
 
        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(geom, WKB_SFSQL, &wkb_len);
        lwgeom_free(geom);
 
        /* rasterize geometry */
        arg->mask = rt_raster_gdal_rasterize(
                wkb, wkb_len,
                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);
 
        /* run iterator */
 
        /* init itrset */
        itrset = palloc(sizeof(struct rt_iterator_t) * 2);
        if (itrset == NULL) {
                rtpg_clip_arg_destroy(arg);
                PG_FREE_IF_COPY(pgraster, 0);
                elog(ERROR, "RASTER_clip: Could not allocate memory for iterator arguments");
                PG_RETURN_NULL();
        }
 
        /* one band at a time */
        for (i = 0; i < arg->numbands; i++) {
                POSTGIS_RT_DEBUGF(4, "band arg %d (nband, hasnodata, nodataval) = (%d, %d, %f)",
                        i, arg->band[i].nband, arg->band[i].hasnodata, arg->band[i].nodataval);
 
                band = rt_raster_get_band(arg->raster, arg->band[i].nband);
 
                /* band metadata */
                pixtype = rt_band_get_pixtype(band);
 
                if (arg->band[i].hasnodata) {
                        hasnodata = 1;
                        nodataval = arg->band[i].nodataval;
                }
                else if (rt_band_get_hasnodata_flag(band)) {
                        hasnodata = 1;
                        rt_band_get_nodata(band, &nodataval);
                }
                else {
                        hasnodata = 0;
                        nodataval = rt_band_get_min_value(band);
                }
 
                /* band is NODATA, create NODATA band and continue */
                if (rt_band_get_isnodata_flag(band)) {
                        /* create raster */
                        if (rtn == NULL) {
                                noerr = rt_raster_from_two_rasters(arg->raster, arg->mask, arg->extenttype, &rtn, NULL);
                                if (noerr != ES_NONE) {
                                        rtpg_clip_arg_destroy(arg);
                                        PG_FREE_IF_COPY(pgraster, 0);
                                        elog(ERROR, "RASTER_clip: Could not create output raster");
                                        PG_RETURN_NULL();
                                }
                        }
 
                        /* create NODATA band */
                        if (rt_raster_generate_new_band(rtn, pixtype, nodataval, hasnodata, nodataval, i) < 0) {
                                rt_raster_destroy(rtn);
                                rtpg_clip_arg_destroy(arg);
                                PG_FREE_IF_COPY(pgraster, 0);
                                elog(ERROR, "RASTER_clip: Could not add NODATA band to output raster");
                                PG_RETURN_NULL();
                        }
 
                        continue;
                }
 
                /* raster */
                itrset[0].raster = arg->raster;
                itrset[0].nband = arg->band[i].nband;
                itrset[0].nbnodata = 1;
 
                /* mask */
                itrset[1].raster = arg->mask;
                itrset[1].nband = 0;
                itrset[1].nbnodata = 1;
 
                /* pass to iterator */
                noerr = rt_raster_iterator(
                        itrset, 2,
                        arg->extenttype, NULL,
                        pixtype,
                        hasnodata, nodataval,
                        0, 0,
                        NULL,
                        NULL,
                        rtpg_clip_callback,
                        &_raster
                );
 
                if (noerr != ES_NONE) {
                        pfree(itrset);
                        rtpg_clip_arg_destroy(arg);
                        if (rtn != NULL) rt_raster_destroy(rtn);
                        PG_FREE_IF_COPY(pgraster, 0);
                        elog(ERROR, "RASTER_clip: Could not run raster iterator function");
                        PG_RETURN_NULL();
                }
 
                /* new raster */
                if (rtn == NULL)
                        rtn = _raster;
                /* copy band */
                else {
                        band = rt_raster_get_band(_raster, 0);
                        if (band == NULL) {
                                pfree(itrset);
                                rtpg_clip_arg_destroy(arg);
                                rt_raster_destroy(_raster);
                                rt_raster_destroy(rtn);
                                PG_FREE_IF_COPY(pgraster, 0);
                                elog(NOTICE, "RASTER_clip: Could not get band from working raster");
                                PG_RETURN_NULL();
                        }
 
                        if (rt_raster_add_band(rtn, band, i) < 0) {
                                pfree(itrset);
                                rtpg_clip_arg_destroy(arg);
                                rt_raster_destroy(_raster);
                                rt_raster_destroy(rtn);
                                PG_FREE_IF_COPY(pgraster, 0);
                                elog(ERROR, "RASTER_clip: Could not add new band to output raster");
                                PG_RETURN_NULL();
                        }
 
                        rt_raster_destroy(_raster);
                }
        }
 
        pfree(itrset);
        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, ovdump::band, clamp_srid(), 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(), rtpg_clip_arg_t::mask, rt_iterator_t::nband, rtpg_clip_band_t::nband, rt_iterator_t::nbnodata, rtpg_clip_band_t::nodataval, rtpg_clip_arg_t::numbands, POSTGIS_RT_DEBUG, POSTGIS_RT_DEBUGF, rt_iterator_t::raster, 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_pixtype(), rt_raster_add_band(), 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_num_bands(), rt_raster_get_srid(), rt_raster_has_band(), rt_raster_is_empty(), rt_raster_iterator(), rt_raster_new(), rt_raster_serialize(), rt_raster_set_srid(), rtpg_clip_arg_destroy(), rtpg_clip_arg_init(), rtpg_clip_callback(), rt_raster_serialized_t::size, SRID_UNKNOWN, and WKB_SFSQL.

Here is the call graph for this function: