RASTER_setPixelValuesGeomval()

Datum RASTER_setPixelValuesGeomval

(

PG_FUNCTION_ARGS

)

Definition at line 1377 of file rtpg_pixel.c.

{
        rt_pgraster *pgraster = NULL;
        rt_pgraster *pgrtn = NULL;
        rt_raster raster = NULL;
        rt_band band = NULL;
        rt_raster _raster = NULL;
        rt_band _band = NULL;
        int nband = 0; /* 1-based */
 
        int numbands = 0;
        int width = 0;
        int height = 0;
        int32_t srid = 0;
        double gt[6] = {0};
 
        rt_pixtype pixtype = PT_END;
        int hasnodata = 0;
        double nodataval = 0;
 
        rtpg_setvaluesgv_arg arg = NULL;
        int allpoint = 0;
 
        ArrayType *array;
        Oid etype;
        Datum *e;
        bool *nulls;
        int16 typlen;
        bool typbyval;
        char typalign;
        int n = 0;
 
        HeapTupleHeader tup;
        bool isnull;
        Datum tupv;
 
        GSERIALIZED *gser = NULL;
        uint8_t gtype;
        lwvarlena_t *wkb = NULL;
 
        int i = 0;
        uint32_t j = 0;
        int noerr = 1;
 
        /* pgraster is null, return null */
        if (PG_ARGISNULL(0))
                PG_RETURN_NULL();
        pgraster = (rt_pgraster *) PG_DETOAST_DATUM_COPY(PG_GETARG_DATUM(0));
 
        /* raster */
        raster = rt_raster_deserialize(pgraster, FALSE);
        if (!raster) {
                PG_FREE_IF_COPY(pgraster, 0);
                elog(ERROR, "RASTER_setPixelValuesGeomval: Could not deserialize raster");
                PG_RETURN_NULL();
        }
 
        /* raster attributes */
        numbands = rt_raster_get_num_bands(raster);
        width = rt_raster_get_width(raster);
        height = rt_raster_get_height(raster);
        srid = clamp_srid(rt_raster_get_srid(raster));
        rt_raster_get_geotransform_matrix(raster, gt);
 
        /* nband */
        if (PG_ARGISNULL(1)) {
                elog(NOTICE, "Band index cannot be NULL.  Value must be 1-based.  Returning original raster");
                rt_raster_destroy(raster);
                PG_RETURN_POINTER(pgraster);
        }
 
        nband = PG_GETARG_INT32(1);
        if (nband < 1 || nband > numbands) {
                elog(NOTICE, "Band index is invalid.  Value must be 1-based.  Returning original raster");
                rt_raster_destroy(raster);
                PG_RETURN_POINTER(pgraster);
        }
 
        /* get band attributes */
        band = rt_raster_get_band(raster, nband - 1);
        pixtype = rt_band_get_pixtype(band);
        hasnodata = rt_band_get_hasnodata_flag(band);
        if (hasnodata)
                rt_band_get_nodata(band, &nodataval);
 
        /* array of geomval (2) */
        if (PG_ARGISNULL(2)) {
                elog(NOTICE, "No values to set.  Returning original raster");
                rt_raster_destroy(raster);
                PG_RETURN_POINTER(pgraster);
        }
 
        array = PG_GETARG_ARRAYTYPE_P(2);
        etype = ARR_ELEMTYPE(array);
        get_typlenbyvalalign(etype, &typlen, &typbyval, &typalign);
 
        deconstruct_array(
                array,
                etype,
                typlen, typbyval, typalign,
                &e, &nulls, &n
        );
 
        if (!n) {
                elog(NOTICE, "No values to set.  Returning original raster");
                rt_raster_destroy(raster);
                PG_RETURN_POINTER(pgraster);
        }
 
        /* init arg */
        arg = rtpg_setvaluesgv_arg_init();
        if (arg == NULL) {
                rt_raster_destroy(raster);
                PG_FREE_IF_COPY(pgraster, 0);
                elog(ERROR, "RASTER_setPixelValuesGeomval: Could not intialize argument structure");
                PG_RETURN_NULL();
        }
 
        arg->gv = palloc(sizeof(struct rtpg_setvaluesgv_geomval_t) * n);
        if (arg->gv == NULL) {
                rtpg_setvaluesgv_arg_destroy(arg);
                rt_raster_destroy(raster);
                PG_FREE_IF_COPY(pgraster, 0);
                elog(ERROR, "RASTER_setPixelValuesGeomval: Could not allocate memory for geomval array");
                PG_RETURN_NULL();
        }
 
        /* process each element */
        arg->ngv = 0;
        for (i = 0; i < n; i++) {
                if (nulls[i])
                        continue;
 
                arg->gv[arg->ngv].pixval.nodata = 0;
                arg->gv[arg->ngv].pixval.value = 0;
                arg->gv[arg->ngv].geom = NULL;
                arg->gv[arg->ngv].mask = NULL;
 
                /* each element is a tuple */
                tup = (HeapTupleHeader) DatumGetPointer(e[i]);
                if (NULL == tup) {
                        rtpg_setvaluesgv_arg_destroy(arg);
                        rt_raster_destroy(raster);
                        PG_FREE_IF_COPY(pgraster, 0);
                        elog(ERROR, "RASTER_setPixelValuesGeomval: Invalid argument for geomval at index %d", i);
                        PG_RETURN_NULL();
                }
 
                /* first element, geometry */
                POSTGIS_RT_DEBUG(4, "Processing first element (geometry)");
                tupv = GetAttributeByName(tup, "geom", &isnull);
                if (isnull) {
                        elog(NOTICE, "First argument (geom) of geomval at index %d is NULL. Skipping", i);
                        continue;
                }
 
                gser = (GSERIALIZED *) PG_DETOAST_DATUM(tupv);
                arg->gv[arg->ngv].geom = lwgeom_from_gserialized(gser);
                if (arg->gv[arg->ngv].geom == NULL) {
                        rtpg_setvaluesgv_arg_destroy(arg);
                        rt_raster_destroy(raster);
                        PG_FREE_IF_COPY(pgraster, 0);
                        elog(ERROR, "RASTER_setPixelValuesGeomval: Could not deserialize geometry of geomval at index %d", i);
                        PG_RETURN_NULL();
                }
 
                /* empty geometry */
                if (lwgeom_is_empty(arg->gv[arg->ngv].geom)) {
                        elog(NOTICE, "First argument (geom) of geomval at index %d is an empty geometry. Skipping", i);
                        continue;
                }
 
                /* check SRID */
                if (clamp_srid(gserialized_get_srid(gser)) != srid) {
                        elog(NOTICE, "Geometry provided for geomval at index %d does not have the same SRID as the raster: %d. Returning original raster", i, srid);
                        rtpg_setvaluesgv_arg_destroy(arg);
                        rt_raster_destroy(raster);
                        PG_RETURN_POINTER(pgraster);
                }
 
                /* Get a 2D version of the geometry if necessary */
                if (lwgeom_ndims(arg->gv[arg->ngv].geom) > 2) {
                        LWGEOM *geom2d = lwgeom_force_2d(arg->gv[arg->ngv].geom);
                        lwgeom_free(arg->gv[arg->ngv].geom);
                        arg->gv[arg->ngv].geom = geom2d;
                }
 
                /* filter for types */
                gtype = gserialized_get_type(gser);
 
                /* shortcuts for POINT and MULTIPOINT */
                if (gtype == POINTTYPE || gtype == MULTIPOINTTYPE)
                        allpoint++;
 
                /* get wkb of geometry */
                POSTGIS_RT_DEBUG(3, "getting wkb of geometry");
                wkb = lwgeom_to_wkb_varlena(arg->gv[arg->ngv].geom, WKB_SFSQL);
 
                /* rasterize geometry */
                arg->gv[arg->ngv].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 (gtype != POINTTYPE && gtype != MULTIPOINTTYPE) {
                        lwgeom_free(arg->gv[arg->ngv].geom);
                        arg->gv[arg->ngv].geom = NULL;
                }
 
                if (arg->gv[arg->ngv].mask == NULL) {
                        rtpg_setvaluesgv_arg_destroy(arg);
                        rt_raster_destroy(raster);
                        PG_FREE_IF_COPY(pgraster, 0);
                        elog(ERROR, "RASTER_setPixelValuesGeomval: Could not rasterize geometry of geomval at index %d", i);
                        PG_RETURN_NULL();
                }
 
                /* set SRID */
                rt_raster_set_srid(arg->gv[arg->ngv].mask, srid);
 
                /* second element, value */
                POSTGIS_RT_DEBUG(4, "Processing second element (val)");
                tupv = GetAttributeByName(tup, "val", &isnull);
                if (isnull) {
                        elog(NOTICE, "Second argument (val) of geomval at index %d is NULL. Treating as NODATA", i);
                        arg->gv[arg->ngv].pixval.nodata = 1;
                }
                else
                        arg->gv[arg->ngv].pixval.value = DatumGetFloat8(tupv);
 
                (arg->ngv)++;
        }
 
        /* redim arg->gv if needed */
        if (arg->ngv < n) {
                arg->gv = repalloc(arg->gv, sizeof(struct rtpg_setvaluesgv_geomval_t) * arg->ngv);
                if (arg->gv == NULL) {
                        rtpg_setvaluesgv_arg_destroy(arg);
                        rt_raster_destroy(raster);
                        PG_FREE_IF_COPY(pgraster, 0);
                        elog(ERROR, "RASTER_setPixelValuesGeomval: Could not reallocate memory for geomval array");
                        PG_RETURN_NULL();
                }
        }
 
        /* keepnodata */
        if (!PG_ARGISNULL(3))
                arg->keepnodata = PG_GETARG_BOOL(3);
        POSTGIS_RT_DEBUGF(3, "keepnodata = %d", arg->keepnodata);
 
        /* keepnodata = TRUE and band is NODATA */
        if (arg->keepnodata && rt_band_get_isnodata_flag(band)) {
                POSTGIS_RT_DEBUG(3, "keepnodata = TRUE and band is NODATA. Not doing anything");
        }
        /* all elements are points */
        else if (allpoint == arg->ngv) {
                double igt[6] = {0};
                double xy[2] = {0};
                double value = 0;
                int isnodata = 0;
 
                LWCOLLECTION *coll = NULL;
                LWPOINT *point = NULL;
                POINT2D p;
 
                POSTGIS_RT_DEBUG(3, "all geometries are points, using direct to pixel method");
 
                /* cache inverse gretransform matrix */
                rt_raster_get_inverse_geotransform_matrix(NULL, gt, igt);
 
                /* process each geometry */
                for (i = 0; i < arg->ngv; i++) {
                        /* convert geometry to collection */
                        coll = lwgeom_as_lwcollection(lwgeom_as_multi(arg->gv[i].geom));
 
                        /* process each point in collection */
                        for (j = 0; j < coll->ngeoms; j++) {
                                point = lwgeom_as_lwpoint(coll->geoms[j]);
                                getPoint2d_p(point->point, 0, &p);
 
                                if (rt_raster_geopoint_to_cell(raster, p.x, p.y, &(xy[0]), &(xy[1]), igt) != ES_NONE) {
                                        rtpg_setvaluesgv_arg_destroy(arg);
                                        rt_raster_destroy(raster);
                                        PG_FREE_IF_COPY(pgraster, 0);
                                        elog(ERROR, "RASTER_setPixelValuesGeomval: Could not process coordinates of point");
                                        PG_RETURN_NULL();
                                }
 
                                /* skip point if outside raster */
                                if (
                                        (xy[0] < 0 || xy[0] >= width) ||
                                        (xy[1] < 0 || xy[1] >= height)
                                ) {
                                        elog(NOTICE, "Point is outside raster extent. Skipping");
                                        continue;
                                }
 
                                /* get pixel value */
                                if (rt_band_get_pixel(band, xy[0], xy[1], &value, &isnodata) != ES_NONE) {
                                        rtpg_setvaluesgv_arg_destroy(arg);
                                        rt_raster_destroy(raster);
                                        PG_FREE_IF_COPY(pgraster, 0);
                                        elog(ERROR, "RASTER_setPixelValuesGeomval: Could not get pixel value");
                                        PG_RETURN_NULL();
                                }
 
                                /* keepnodata = TRUE AND pixel value is NODATA */
                                if (arg->keepnodata && isnodata)
                                        continue;
 
                                /* set pixel */
                                if (arg->gv[i].pixval.nodata)
                                        noerr = rt_band_set_pixel(band, xy[0], xy[1], nodataval, NULL);
                                else
                                        noerr = rt_band_set_pixel(band, xy[0], xy[1], arg->gv[i].pixval.value, NULL);
 
                                if (noerr != ES_NONE) {
                                        rtpg_setvaluesgv_arg_destroy(arg);
                                        rt_raster_destroy(raster);
                                        PG_FREE_IF_COPY(pgraster, 0);
                                        elog(ERROR, "RASTER_setPixelValuesGeomval: Could not set pixel value");
                                        PG_RETURN_NULL();
                                }
                        }
                }
        }
        /* run iterator otherwise */
        else {
                rt_iterator itrset;
 
                POSTGIS_RT_DEBUG(3, "a mix of geometries, using iterator method");
 
                /* init itrset */
                itrset = palloc(sizeof(struct rt_iterator_t) * (arg->ngv + 1));
                if (itrset == NULL) {
                        rtpg_setvaluesgv_arg_destroy(arg);
                        rt_raster_destroy(raster);
                        PG_FREE_IF_COPY(pgraster, 0);
                        elog(ERROR, "RASTER_setPixelValuesGeomval: Could not allocate memory for iterator arguments");
                        PG_RETURN_NULL();
                }
 
                /* set first raster's info */
                itrset[0].raster = raster;
                itrset[0].nband = nband - 1;
                itrset[0].nbnodata = 1;
 
                /* set other raster's info */
                for (i = 0, j = 1; i < arg->ngv; i++, j++) {
                        itrset[j].raster = arg->gv[i].mask;
                        itrset[j].nband = 0;
                        itrset[j].nbnodata = 1;
                }
 
                /* pass to iterator */
                noerr = rt_raster_iterator(
                        itrset, arg->ngv + 1,
                        ET_FIRST, NULL,
                        pixtype,
                        hasnodata, nodataval,
                        0, 0,
                        NULL,
                        arg,
                        rtpg_setvalues_geomval_callback,
                        &_raster
                );
                pfree(itrset);
 
                if (noerr != ES_NONE) {
                        rtpg_setvaluesgv_arg_destroy(arg);
                        rt_raster_destroy(raster);
                        PG_FREE_IF_COPY(pgraster, 0);
                        elog(ERROR, "RASTER_setPixelValuesGeomval: Could not run raster iterator function");
                        PG_RETURN_NULL();
                }
 
                /* copy band from _raster to raster */
                _band = rt_raster_get_band(_raster, 0);
                if (_band == NULL) {
                        rtpg_setvaluesgv_arg_destroy(arg);
                        rt_raster_destroy(_raster);
                        rt_raster_destroy(raster);
                        PG_FREE_IF_COPY(pgraster, 0);
                        elog(ERROR, "RASTER_setPixelValuesGeomval: Could not get band from working raster");
                        PG_RETURN_NULL();
                }
 
                _band = rt_raster_replace_band(raster, _band, nband - 1);
                if (_band == NULL) {
                        rtpg_setvaluesgv_arg_destroy(arg);
                        rt_raster_destroy(_raster);
                        rt_raster_destroy(raster);
                        PG_FREE_IF_COPY(pgraster, 0);
                        elog(ERROR, "RASTER_setPixelValuesGeomval: Could not replace band in output raster");
                        PG_RETURN_NULL();
                }
 
                rt_band_destroy(_band);
                rt_raster_destroy(_raster);
        }
 
        rtpg_setvaluesgv_arg_destroy(arg);
 
        pgrtn = rt_raster_serialize(raster);
        rt_raster_destroy(raster);
        PG_FREE_IF_COPY(pgraster, 0);
 
        POSTGIS_RT_DEBUG(3, "Finished");
 
        if (!pgrtn)
                PG_RETURN_NULL();
 
        SET_VARSIZE(pgrtn, pgrtn->size);
        PG_RETURN_POINTER(pgrtn);
}

References ovdump::band, clamp_srid(), lwvarlena_t::data, ES_NONE, ET_FIRST, FALSE, rtpg_setvaluesgv_geomval_t::geom, LWCOLLECTION::geoms, getPoint2d_p(), gserialized_get_srid(), gserialized_get_type(), window::gt, rtpg_setvaluesgv_arg_t::gv, rtpg_setvaluesgv_arg_t::keepnodata, lwgeom_as_lwcollection(), lwgeom_as_lwpoint(), lwgeom_as_multi(), lwgeom_force_2d(), lwgeom_free(), lwgeom_from_gserialized(), lwgeom_is_empty(), lwgeom_ndims(), lwgeom_to_wkb_varlena(), LWSIZE_GET, LWVARHDRSZ, rtpg_setvaluesgv_geomval_t::mask, MULTIPOINTTYPE, rt_iterator_t::nband, pixval::nband, rt_iterator_t::nbnodata, LWCOLLECTION::ngeoms, rtpg_setvaluesgv_arg_t::ngv, rtpg_setvaluesgv_geomval_t::nodata, rtpg_setvaluesgv_geomval_t::pixval, LWPOINT::point, POINTTYPE, POSTGIS_RT_DEBUG, POSTGIS_RT_DEBUGF, PT_END, rt_iterator_t::raster, rtrowdump::raster, rt_band_destroy(), rt_band_get_hasnodata_flag(), rt_band_get_isnodata_flag(), rt_band_get_nodata(), rt_band_get_pixel(), rt_band_get_pixtype(), rt_band_set_pixel(), rt_raster_deserialize(), rt_raster_destroy(), rt_raster_gdal_rasterize(), rt_raster_geopoint_to_cell(), rt_raster_get_band(), rt_raster_get_geotransform_matrix(), rt_raster_get_height(), rt_raster_get_inverse_geotransform_matrix(), rt_raster_get_num_bands(), rt_raster_get_srid(), rt_raster_get_width(), rt_raster_iterator(), rt_raster_replace_band(), rt_raster_serialize(), rt_raster_set_srid(), rtpg_setvalues_geomval_callback(), rtpg_setvaluesgv_arg_destroy(), rtpg_setvaluesgv_arg_init(), lwvarlena_t::size, rt_raster_serialized_t::size, rtpg_setvaluesgv_geomval_t::value, genraster::value, WKB_SFSQL, POINT2D::x, and POINT2D::y.

Here is the call graph for this function: