From: Bborie Park Date: Thu, 18 Oct 2012 21:45:06 +0000 (+0000) Subject: Addition of geomval variants of ST_SetValues() and regression tests. X-Git-Tag: 2.1.0beta2~511 X-Git-Url: https://granicus.if.org/sourcecode?a=commitdiff_plain;h=8da4897b8cadd0c6bc1c3ceebc79d27c91f6dae8;p=postgis Addition of geomval variants of ST_SetValues() and regression tests. Added helper function rt_raster_get_inverse_geotransform_matrix(). Additional code cleanup for rt_raster_geopoint_to_cell() and rt_raster_cell_to_geopoint(). git-svn-id: http://svn.osgeo.org/postgis/trunk@10465 b70326c6-7e19-0410-871a-916f4a2858ee --- diff --git a/raster/rt_core/rt_api.c b/raster/rt_core/rt_api.c index 3f359c5d2..a87512888 100644 --- a/raster/rt_core/rt_api.c +++ b/raster/rt_core/rt_api.c @@ -5583,6 +5583,35 @@ rt_raster_generate_new_band(rt_raster raster, rt_pixtype pixtype, return index; } +/** + * Get 6-element array of raster inverse geotransform matrix + * + * @param raster : the raster to get matrix of + * @param gt : optional input parameter, 6-element geotransform matrix + * @param igt : output parameter, 6-element inverse geotransform matrix + * + * @return if zero, error occurred in function + */ +int rt_raster_get_inverse_geotransform_matrix(rt_raster raster, + double *gt, double *igt) { + double _gt[6] = {0}; + + assert((raster != NULL || gt != NULL)); + assert(igt != NULL); + + if (gt == NULL) + rt_raster_get_geotransform_matrix(raster, _gt); + else + memcpy(_gt, gt, sizeof(double) * 6); + + if (!GDALInvGeoTransform(_gt, igt)) { + rterror("rt_raster_get_inverse_geotransform_matrix: Unable to compute inverse geotransform matrix"); + return 0; + } + + return 1; +} + /** * Get 6-element array of raster geotransform matrix * @@ -5643,28 +5672,14 @@ rt_raster_cell_to_geopoint(rt_raster raster, double *xw, double *yw, double *gt ) { - double *_gt = NULL; - int init_gt = 0; - int i = 0; + double _gt[6] = {0}; assert(NULL != raster); assert(NULL != xw); assert(NULL != yw); - if (NULL == gt) { - _gt = rtalloc(sizeof(double) * 6); - if (NULL == _gt) { - rterror("rt_raster_cell_to_geopoint: Unable to allocate memory for geotransform matrix"); - return 0; - } - init_gt = 1; - - for (i = 0; i < 6; i++) _gt[i] = 0; - } - else { - _gt = gt; - init_gt = 0; - } + if (NULL != gt) + memcpy(_gt, gt, sizeof(double) * 6); /* scale of matrix is not set */ if ( @@ -5687,7 +5702,6 @@ rt_raster_cell_to_geopoint(rt_raster raster, RASTER_DEBUGF(4, "GDALApplyGeoTransform (c -> g) for (%f, %f) = (%f, %f)", xr, yr, *xw, *yw); - if (init_gt) rtdealloc(_gt); return 1; } @@ -5699,7 +5713,7 @@ rt_raster_cell_to_geopoint(rt_raster raster, * @param yw : Y ordinate of the geographical point * @param xr : output parameter, the pixel's column * @param yr : output parameter, the pixel's row - * @param igt : input/output parameter, 3x2 inverse geotransform matrix + * @param igt : input/output parameter, inverse geotransform matrix * * @return if zero, error occurred in function */ @@ -5709,29 +5723,15 @@ rt_raster_geopoint_to_cell(rt_raster raster, double *xr, double *yr, double *igt ) { - double *_igt = NULL; - int i = 0; - int init_igt = 0; + double _igt[6] = {0}; double rnd = 0; assert(NULL != raster); assert(NULL != xr); assert(NULL != yr); - if (igt == NULL) { - _igt = rtalloc(sizeof(double) * 6); - if (_igt == NULL) { - rterror("rt_raster_geopoint_to_cell: Unable to allocate memory for inverse geotransform matrix"); - return 0; - } - init_igt = 1; - - for (i = 0; i < 6; i++) _igt[i] = 0; - } - else { - _igt = igt; - init_igt = 0; - } + if (igt != NULL) + memcpy(_igt, igt, sizeof(double) * 6); /* matrix is not set */ if ( @@ -5742,12 +5742,8 @@ rt_raster_geopoint_to_cell(rt_raster raster, FLT_EQ(_igt[4], 0.) && FLT_EQ(_igt[5], 0.) ) { - double gt[6] = {0.0}; - rt_raster_get_geotransform_matrix(raster, gt); - - if (!GDALInvGeoTransform(gt, _igt)) { - rterror("rt_raster_geopoint_to_cell: Unable to compute inverse geotransform matrix"); - if (init_igt) rtdealloc(_igt); + if (!rt_raster_get_inverse_geotransform_matrix(raster, NULL, _igt)) { + rterror("rt_raster_geopoint_to_cell: Unable to get inverse geotransform matrix"); return 0; } } @@ -5771,7 +5767,6 @@ rt_raster_geopoint_to_cell(rt_raster raster, RASTER_DEBUGF(4, "Corrected GDALApplyGeoTransform (g -> c) for (%f, %f) = (%f, %f)", xw, yw, *xr, *yr); - if (init_igt) rtdealloc(_igt); return 1; } diff --git a/raster/rt_core/rt_api.h b/raster/rt_core/rt_api.h index ac624d6f4..83bd94fc4 100644 --- a/raster/rt_core/rt_api.h +++ b/raster/rt_core/rt_api.h @@ -1114,6 +1114,18 @@ void rt_raster_set_srid(rt_raster raster, int32_t srid); */ int32_t rt_raster_get_srid(rt_raster raster); +/** + * Get 6-element array of raster inverse geotransform matrix + * + * @param raster : the raster to get matrix of + * @param gt : optional input parameter, 6-element geotransform matrix + * @param igt : output parameter, 6-element inverse geotransform matrix + * + * @return if zero, error occurred in function + */ +int rt_raster_get_inverse_geotransform_matrix(rt_raster raster, + double *gt, double *igt); + /** * Get 6-element array of raster geotransform matrix * diff --git a/raster/rt_pg/rt_pg.c b/raster/rt_pg/rt_pg.c index 129ef6223..2c8ffa180 100644 --- a/raster/rt_pg/rt_pg.c +++ b/raster/rt_pg/rt_pg.c @@ -247,6 +247,7 @@ Datum RASTER_dumpValues(PG_FUNCTION_ARGS); /* Set pixel value(s) */ Datum RASTER_setPixelValue(PG_FUNCTION_ARGS); Datum RASTER_setPixelValuesArray(PG_FUNCTION_ARGS); +Datum RASTER_setPixelValuesGeomval(PG_FUNCTION_ARGS); /* Get pixel geographical shape */ Datum RASTER_getPixelPolygons(PG_FUNCTION_ARGS); @@ -2435,7 +2436,7 @@ Datum RASTER_getPixelValue(PG_FUNCTION_ARGS) } /* ---------------------------------------------------------------- */ -/* ST_DumpValue function */ +/* ST_DumpValues function */ /* ---------------------------------------------------------------- */ typedef struct rtpg_dumpvalues_arg_t *rtpg_dumpvalues_arg; @@ -2947,6 +2948,9 @@ Datum RASTER_setPixelValue(PG_FUNCTION_ARGS) PG_RETURN_POINTER(pgrtn); } +/** + * Set pixels to value from array + */ PG_FUNCTION_INFO_V1(RASTER_setPixelValuesArray); Datum RASTER_setPixelValuesArray(PG_FUNCTION_ARGS) { @@ -3368,6 +3372,511 @@ Datum RASTER_setPixelValuesArray(PG_FUNCTION_ARGS) PG_RETURN_POINTER(pgrtn); } +/* ---------------------------------------------------------------- */ +/* ST_SetValues using geomval array */ +/* ---------------------------------------------------------------- */ + +typedef struct rtpg_setvaluesgv_arg_t *rtpg_setvaluesgv_arg; +typedef struct rtpg_setvaluesgv_geomval_t *rtpg_setvaluesgv_geomval; + +struct rtpg_setvaluesgv_arg_t { + int ngv; + rtpg_setvaluesgv_geomval gv; + + bool keepnodata; +}; + +struct rtpg_setvaluesgv_geomval_t { + struct { + int nodata; + double value; + } pixval; + + LWGEOM *geom; + rt_raster mask; +}; + +static rtpg_setvaluesgv_arg rtpg_setvaluesgv_arg_init() { + rtpg_setvaluesgv_arg arg = palloc(sizeof(struct rtpg_setvaluesgv_arg_t)); + if (arg == NULL) { + elog(ERROR, "rtpg_setvaluesgv_arg_init: Unable to allocate memory for function arguments"); + return NULL; + } + + arg->ngv = 0; + arg->gv = NULL; + arg->keepnodata = 0; + + return arg; +} + +static void rtpg_setvaluesgv_arg_destroy(rtpg_setvaluesgv_arg arg) { + int i = 0; + + if (arg->gv != NULL) { + for (i = 0; i < arg->ngv; i++) { + if (arg->gv[i].geom != NULL) + lwgeom_free(arg->gv[i].geom); + if (arg->gv[i].mask != NULL) + rt_raster_destroy(arg->gv[i].mask); + } + + pfree(arg->gv); + } + + pfree(arg); +} + +static int rtpg_setvalues_geomval_callback( + rt_iterator_arg arg, void *userarg, + double *value, int *nodata +) { + rtpg_setvaluesgv_arg funcarg = (rtpg_setvaluesgv_arg) userarg; + int i = 0; + int j = 0; + + *value = 0; + *nodata = 0; + + POSTGIS_RT_DEBUGF(4, "keepnodata = %d", funcarg->keepnodata); + + /* keepnodata = TRUE */ + if (funcarg->keepnodata && arg->nodata[0][0][0]) { + POSTGIS_RT_DEBUG(3, "keepnodata = 1 AND source raster pixel is NODATA"); + *nodata = 1; + return 1; + } + + for (i = arg->rasters - 1, j = funcarg->ngv - 1; i > 0; i--, j--) { + POSTGIS_RT_DEBUGF(4, "checking raster %d", i); + + /* mask is NODATA */ + if (arg->nodata[i][0][0]) + continue; + /* mask is NOT NODATA */ + else { + POSTGIS_RT_DEBUGF(4, "Using information from geometry %d", j); + + if (funcarg->gv[j].pixval.nodata) + *nodata = 1; + else + *value = funcarg->gv[j].pixval.value; + + return 1; + } + } + + POSTGIS_RT_DEBUG(4, "Using information from source raster"); + + /* still here */ + if (arg->nodata[0][0][0]) + *nodata = 1; + else + *value = arg->values[0][0][0]; + + return 1; +} + +PG_FUNCTION_INFO_V1(RASTER_setPixelValuesGeomval); +Datum RASTER_setPixelValuesGeomval(PG_FUNCTION_ARGS) +{ + 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 */ + rt_iterator itrset; + + int numbands = 0; + int width = 0; + int height = 0; + int 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; + unsigned char *wkb = NULL; + size_t wkb_len; + + int i = 0; + int 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) { + elog(ERROR, "RASTER_setPixelValuesGeomval: Could not deserialize raster"); + PG_FREE_IF_COPY(pgraster, 0); + 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 = 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) + nodataval = rt_band_get_nodata(band); + + /* 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) { + elog(ERROR, "RASTER_setPixelValuesGeomval: Unable to intialize argument structure"); + rt_raster_destroy(raster); + PG_FREE_IF_COPY(pgraster, 0); + PG_RETURN_NULL(); + } + + arg->gv = palloc(sizeof(struct rtpg_setvaluesgv_geomval_t) * n); + if (arg->gv == NULL) { + elog(ERROR, "RASTER_setPixelValuesGeomval: Unable to allocate memory for geomval array"); + rtpg_setvaluesgv_arg_destroy(arg); + rt_raster_destroy(raster); + PG_FREE_IF_COPY(pgraster, 0); + 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) { + elog(ERROR, "RASTER_setPixelValuesGeomval: Invalid argument for geomval at index %d", i); + rtpg_setvaluesgv_arg_destroy(arg); + rt_raster_destroy(raster); + PG_FREE_IF_COPY(pgraster, 0); + 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) { + elog(ERROR, "RASTER_setPixelValuesGeomval: Unable to deserialize geometry of geomval at index %d", i); + rtpg_setvaluesgv_arg_destroy(arg); + rt_raster_destroy(raster); + PG_FREE_IF_COPY(pgraster, 0); + 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 (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); + } + + /* 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(arg->gv[arg->ngv].geom, WKB_SFSQL, &wkb_len); + + /* rasterize geometry */ + arg->gv[arg->ngv].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 (gtype != POINTTYPE && gtype != MULTIPOINTTYPE) { + lwgeom_free(arg->gv[arg->ngv].geom); + arg->gv[arg->ngv].geom = NULL; + } + + if (arg->gv[arg->ngv].mask == NULL) { + elog(ERROR, "RASTER_setPixelValuesGeomval: Unable to rasterize geometry of geomval at index %d", i); + rtpg_setvaluesgv_arg_destroy(arg); + rt_raster_destroy(raster); + PG_FREE_IF_COPY(pgraster, 0); + PG_RETURN_NULL(); + } + + /* 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) { + elog(ERROR, "RASTER_setPixelValuesGeomval: Unable to reallocate memory for geomval array"); + rtpg_setvaluesgv_arg_destroy(arg); + rt_raster_destroy(raster); + PG_FREE_IF_COPY(pgraster, 0); + PG_RETURN_NULL(); + } + } + + /* keepnodata */ + if (!PG_ARGISNULL(3)) + arg->keepnodata = PG_GETARG_BOOL(3); + POSTGIS_RT_DEBUGF(3, "keepnodata = %d", arg->keepnodata); + + /* all elements are points */ + if (allpoint == arg->ngv) { + POSTGIS_RT_DEBUG(3, "all geometries are points, using direct to pixel method"); + double igt[6] = {0}; + double xy[2] = {0}; + double value = 0; + + LWCOLLECTION *coll = NULL; + LWPOINT *point = NULL; + POINT2D p; + + /* 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)) { + elog(ERROR, "RASTER_setPixelValuesGeomval: Unable to process coordinates of point"); + rtpg_setvaluesgv_arg_destroy(arg); + rt_raster_destroy(raster); + PG_FREE_IF_COPY(pgraster, 0); + 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) != 0) { + elog(ERROR, "RASTER_setPixelValuesGeomval: Unable to get pixel value"); + rtpg_setvaluesgv_arg_destroy(arg); + rt_raster_destroy(raster); + PG_FREE_IF_COPY(pgraster, 0); + PG_RETURN_NULL(); + } + + /* keepnodata = TRUE AND pixel value is NODATA */ + if ( + arg->keepnodata && + hasnodata && ( + FLT_EQ(value, nodataval) || + rt_band_clamped_value_is_nodata(band, value) == 1 + ) + ) { + continue; + } + + /* set pixel */ + if (arg->gv[i].pixval.nodata) + noerr = rt_band_set_pixel(band, xy[0], xy[1], nodataval); + else + noerr = rt_band_set_pixel(band, xy[0], xy[1], arg->gv[i].pixval.value); + + if (noerr < 0) { + elog(ERROR, "RASTER_setPixelValuesGeomval: Unable to set pixel value"); + rtpg_setvaluesgv_arg_destroy(arg); + rt_raster_destroy(raster); + PG_FREE_IF_COPY(pgraster, 0); + PG_RETURN_NULL(); + } + } + } + } + /* run iterator otherwise */ + else { + 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) { + elog(ERROR, "RASTER_setPixelValuesGeomval: Unable to allocate memory for iterator arguments"); + rtpg_setvaluesgv_arg_destroy(arg); + rt_raster_destroy(raster); + PG_FREE_IF_COPY(pgraster, 0); + 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 */ + _raster = rt_raster_iterator( + itrset, arg->ngv + 1, + ET_FIRST, NULL, + pixtype, + hasnodata, nodataval, + 0, 0, + arg, + rtpg_setvalues_geomval_callback, + &noerr + ); + pfree(itrset); + + if (!noerr) { + elog(ERROR, "RASTER_setPixelValuesGeomval: Unable to run raster iterator function"); + rtpg_setvaluesgv_arg_destroy(arg); + rt_raster_destroy(raster); + PG_FREE_IF_COPY(pgraster, 0); + PG_RETURN_NULL(); + } + + /* copy band from _raster to raster */ + _band = rt_raster_get_band(_raster, 0); + _band = rt_raster_replace_band(raster, _band, nband -1); + + 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); +} + /** * Return the geographical shape of all pixels */ @@ -13171,7 +13680,7 @@ Datum RASTER_mapAlgebra2(PG_FUNCTION_ARGS) /* expression has argument(s) */ if (spi_argcount[i]) { /* reset values to (Datum) NULL */ - memset(values, NULL, sizeof(Datum) * argkwcount); + memset(values, (Datum) NULL, sizeof(Datum) * argkwcount); /* reset nulls to FALSE */ memset(nulls, FALSE, sizeof(bool) * argkwcount); @@ -14155,26 +14664,25 @@ static int rtpg_nmapalgebra_rastbandarg_process(rtpg_nmapalgebra_arg arg, ArrayT continue; } - else { - arg->pgraster[i] = (rt_pgraster *) PG_DETOAST_DATUM(tupv); - - /* see if this is a copy of an existing pgraster */ - for (j = 0; j < i; j++) { - if (arg->pgraster[i] == arg->pgraster[j]) { - POSTGIS_RT_DEBUG(4, "raster matching existing same raster found"); - arg->raster[i] = arg->raster[j]; - arg->ownsdata[i] = 0; - break; - } + + arg->pgraster[i] = (rt_pgraster *) PG_DETOAST_DATUM(tupv); + + /* see if this is a copy of an existing pgraster */ + for (j = 0; j < i; j++) { + if (arg->pgraster[i] == arg->pgraster[j]) { + POSTGIS_RT_DEBUG(4, "raster matching existing same raster found"); + arg->raster[i] = arg->raster[j]; + arg->ownsdata[i] = 0; + break; } + } - if (arg->ownsdata[i]) { - POSTGIS_RT_DEBUG(4, "deserializing raster"); - arg->raster[i] = rt_raster_deserialize(arg->pgraster[i], FALSE); - if (arg->raster[i] == NULL) { - elog(ERROR, "rtpg_nmapalgebra_rastbandarg_process: Could not deserialize raster at index %d", i); - return 0; - } + if (arg->ownsdata[i]) { + POSTGIS_RT_DEBUG(4, "deserializing raster"); + arg->raster[i] = rt_raster_deserialize(arg->pgraster[i], FALSE); + if (arg->raster[i] == NULL) { + elog(ERROR, "rtpg_nmapalgebra_rastbandarg_process: Could not deserialize raster at index %d", i); + return 0; } } diff --git a/raster/rt_pg/rtpostgis.sql.in.c b/raster/rt_pg/rtpostgis.sql.in.c index aaf49a06e..28a874c17 100644 --- a/raster/rt_pg/rtpostgis.sql.in.c +++ b/raster/rt_pg/rtpostgis.sql.in.c @@ -85,6 +85,15 @@ CREATE TYPE rastbandarg AS ( nband integer ); +----------------------------------------------------------------------- +-- generic composite type of a geometry and value +----------------------------------------------------------------------- + +CREATE TYPE geomval AS ( + geom geometry, + val double precision +); + ----------------------------------------------------------------------- -- Raster Accessors ----------------------------------------------------------------------- @@ -4102,15 +4111,30 @@ CREATE OR REPLACE FUNCTION st_setvalues( $$ LANGUAGE 'plpgsql' IMMUTABLE; +-- cannot be STRICT as newvalue can be NULL +CREATE OR REPLACE FUNCTION st_setvalues( + rast raster, nband integer, + geomvalset geomval[], + keepnodata boolean DEFAULT FALSE +) + RETURNS raster + AS 'MODULE_PATHNAME', 'RASTER_setPixelValuesGeomval' + LANGUAGE 'c' IMMUTABLE; + +-- cannot be STRICT as newvalue can be NULL +CREATE OR REPLACE FUNCTION st_setvalues( + rast raster, nband integer, + geom geometry, newvalue double precision, + keepnodata boolean DEFAULT FALSE +) + RETURNS raster + AS $$ SELECT st_setvalues($1, $2, ARRAY[ROW($3, $4)]::geomval[], $5) $$ + LANGUAGE 'sql' IMMUTABLE; + ----------------------------------------------------------------------- -- Raster Processing Functions ----------------------------------------------------------------------- -CREATE TYPE geomval AS ( - geom geometry, - val double precision -); - ----------------------------------------------------------------------- -- ST_DumpAsPolygons ----------------------------------------------------------------------- diff --git a/raster/test/regress/Makefile.in b/raster/test/regress/Makefile.in index d507b46d1..ba1407bb9 100644 --- a/raster/test/regress/Makefile.in +++ b/raster/test/regress/Makefile.in @@ -103,7 +103,8 @@ TEST_MAPALGEBRA = \ rt_mapalgebra_expr \ rt_union \ rt_invdistweight4ma \ - rt_4ma + rt_4ma \ + rt_setvalues_geomval TEST_GIST = \ rt_above \ diff --git a/raster/test/regress/rt_setvalues_geomval.sql b/raster/test/regress/rt_setvalues_geomval.sql new file mode 100644 index 000000000..67e1aed31 --- /dev/null +++ b/raster/test/regress/rt_setvalues_geomval.sql @@ -0,0 +1,80 @@ +DROP TABLE IF EXISTS raster_setvalues_rast; +CREATE TABLE raster_setvalues_rast AS + SELECT 1 AS rid, ST_AddBand(ST_MakeEmptyRaster(5, 5, 0, 0, 1, -1, 0, 0, 0), 1, '8BUI', 0, 0) AS rast +; + +DROP TABLE IF EXISTS raster_setvalues_geom; +CREATE TABLE raster_setvalues_geom AS + SELECT 1 AS gid, 'SRID=0;POINT(2.5 -2.5)'::geometry geom UNION ALL + SELECT 2 AS gid, 'SRID=0;POLYGON((1 -1, 4 -1, 4 -4, 1 -4, 1 -1))'::geometry geom UNION ALL + SELECT 3 AS gid, 'SRID=0;POLYGON((0 0, 5 0, 5 -1, 1 -1, 1 -4, 0 -4, 0 0))'::geometry geom UNION ALL + SELECT 4 AS gid, 'SRID=0;MULTIPOINT(0 0, 4 4, 4 -4)'::geometry +; + +SELECT + rid, gid, ST_DumpValues(ST_SetValues(rast, 1, geom, gid)) +FROM raster_setvalues_rast t1 +CROSS JOIN raster_setvalues_geom t2 +ORDER BY rid, gid; + +SELECT + t1.rid, t2.gid, t3.gid, ST_DumpValues(ST_SetValues(rast, 1, ARRAY[ROW(t2.geom, t2.gid), ROW(t3.geom, t3.gid)]::geomval[])) +FROM raster_setvalues_rast t1 +CROSS JOIN raster_setvalues_geom t2 +CROSS JOIN raster_setvalues_geom t3 +WHERE t2.gid = 1 + AND t3.gid = 2 +ORDER BY t1.rid, t2.gid, t3.gid; + +SELECT + t1.rid, t2.gid, t3.gid, ST_DumpValues(ST_SetValues(rast, 1, ARRAY[ROW(t3.geom, t3.gid), ROW(t2.geom, t2.gid)]::geomval[])) +FROM raster_setvalues_rast t1 +CROSS JOIN raster_setvalues_geom t2 +CROSS JOIN raster_setvalues_geom t3 +WHERE t2.gid = 1 + AND t3.gid = 2 +ORDER BY t1.rid, t2.gid, t3.gid; + +SELECT + t1.rid, t2.gid, t3.gid, ST_DumpValues(ST_SetValues(rast, 1, ARRAY[ROW(t3.geom, t3.gid), ROW(t2.geom, t2.gid)]::geomval[])) +FROM raster_setvalues_rast t1 +CROSS JOIN raster_setvalues_geom t2 +CROSS JOIN raster_setvalues_geom t3 +WHERE t2.gid = 1 + AND t3.gid = 3 +ORDER BY t1.rid, t2.gid, t3.gid; + +SELECT + t1.rid, t2.gid, t3.gid, ST_DumpValues(ST_SetValues(rast, 1, ARRAY[ROW(t3.geom, t3.gid), ROW(t2.geom, t2.gid)]::geomval[])) +FROM raster_setvalues_rast t1 +CROSS JOIN raster_setvalues_geom t2 +CROSS JOIN raster_setvalues_geom t3 +WHERE t2.gid = 1 + AND t3.gid = 4; + +WITH foo AS ( + SELECT + array_agg(gid) AS gid, + ST_Union(geom) AS geom + FROM raster_setvalues_geom + WHERE gid IN (1,4) +) +SELECT + t1.rid, t2.gid, ST_DumpValues(ST_SetValues(rast, 1, ARRAY[ROW(t2.geom, 99)]::geomval[])) +FROM raster_setvalues_rast t1 +CROSS JOIN foo t2; + +WITH foo AS ( + SELECT + array_agg(gid) AS gid, + ST_Union(geom) AS geom + FROM raster_setvalues_geom + WHERE gid IN (2,3) +) +SELECT + t1.rid, t2.gid, ST_DumpValues(ST_SetValues(rast, 1, ARRAY[ROW(t2.geom, 99)]::geomval[])) +FROM raster_setvalues_rast t1 +CROSS JOIN foo t2; + +DROP TABLE IF EXISTS raster_setvalues_rast; +DROP TABLE IF EXISTS raster_setvalues_geom; diff --git a/raster/test/regress/rt_setvalues_geomval_expected b/raster/test/regress/rt_setvalues_geomval_expected new file mode 100644 index 000000000..1f95cf29e --- /dev/null +++ b/raster/test/regress/rt_setvalues_geomval_expected @@ -0,0 +1,15 @@ +NOTICE: table "raster_setvalues_rast" does not exist, skipping +NOTICE: table "raster_setvalues_geom" does not exist, skipping +NOTICE: Point is outside raster extent. Skipping +1|1|(1,"{{NULL,NULL,NULL,NULL,NULL},{NULL,NULL,NULL,NULL,NULL},{NULL,NULL,1,NULL,NULL},{NULL,NULL,NULL,NULL,NULL},{NULL,NULL,NULL,NULL,NULL}}") +1|2|(1,"{{NULL,NULL,NULL,NULL,NULL},{NULL,2,2,2,NULL},{NULL,2,2,2,NULL},{NULL,2,2,2,NULL},{NULL,NULL,NULL,NULL,NULL}}") +1|3|(1,"{{3,3,3,3,3},{3,NULL,NULL,NULL,NULL},{3,NULL,NULL,NULL,NULL},{3,NULL,NULL,NULL,NULL},{NULL,NULL,NULL,NULL,NULL}}") +1|4|(1,"{{4,NULL,NULL,NULL,NULL},{NULL,NULL,NULL,NULL,NULL},{NULL,NULL,NULL,NULL,NULL},{NULL,NULL,NULL,NULL,NULL},{NULL,NULL,NULL,NULL,4}}") +1|1|2|(1,"{{NULL,NULL,NULL,NULL,NULL},{NULL,2,2,2,NULL},{NULL,2,2,2,NULL},{NULL,2,2,2,NULL},{NULL,NULL,NULL,NULL,NULL}}") +1|1|2|(1,"{{NULL,NULL,NULL,NULL,NULL},{NULL,2,2,2,NULL},{NULL,2,1,2,NULL},{NULL,2,2,2,NULL},{NULL,NULL,NULL,NULL,NULL}}") +1|1|3|(1,"{{3,3,3,3,3},{3,NULL,NULL,NULL,NULL},{3,NULL,1,NULL,NULL},{3,NULL,NULL,NULL,NULL},{NULL,NULL,NULL,NULL,NULL}}") +NOTICE: Point is outside raster extent. Skipping +1|1|4|(1,"{{4,NULL,NULL,NULL,NULL},{NULL,NULL,NULL,NULL,NULL},{NULL,NULL,1,NULL,NULL},{NULL,NULL,NULL,NULL,NULL},{NULL,NULL,NULL,NULL,4}}") +NOTICE: Point is outside raster extent. Skipping +1|{1,4}|(1,"{{99,NULL,NULL,NULL,NULL},{NULL,NULL,NULL,NULL,NULL},{NULL,NULL,99,NULL,NULL},{NULL,NULL,NULL,NULL,NULL},{NULL,NULL,NULL,NULL,99}}") +1|{2,3}|(1,"{{99,99,99,99,99},{99,99,99,99,NULL},{99,99,99,99,NULL},{99,99,99,99,NULL},{NULL,NULL,NULL,NULL,NULL}}")