PostGIS 3.0.6dev-r@@SVN_REVISION@@
Loading...
Searching...
No Matches

◆ RASTER_clip()

Datum RASTER_clip ( PG_FUNCTION_ARGS  )

Definition at line 3025 of file rtpg_mapalgebra.c.

3026{
3027 rt_pgraster *pgraster = NULL;
3028 LWGEOM *rastgeom = NULL;
3029 double gt[6] = {0};
3030 int32_t srid = SRID_UNKNOWN;
3031
3032 rt_pgraster *pgrtn = NULL;
3033 rt_raster rtn = NULL;
3034
3035 GSERIALIZED *gser = NULL;
3036 LWGEOM *geom = NULL;
3037 unsigned char *wkb = NULL;
3038 size_t wkb_len;
3039
3040 ArrayType *array;
3041 Oid etype;
3042 Datum *e;
3043 bool *nulls;
3044
3045 int16 typlen;
3046 bool typbyval;
3047 char typalign;
3048
3049 int i = 0;
3050 int j = 0;
3051 int k = 0;
3052 rtpg_clip_arg arg = NULL;
3053 LWGEOM *tmpgeom = NULL;
3054 rt_iterator itrset;
3055
3056 rt_raster _raster = NULL;
3057 rt_band band = NULL;
3058 rt_pixtype pixtype;
3059 int hasnodata;
3060 double nodataval;
3061 int noerr = 0;
3062
3063 POSTGIS_RT_DEBUG(3, "Starting...");
3064
3065 /* raster or geometry is NULL, return NULL */
3066 if (PG_ARGISNULL(0) || PG_ARGISNULL(2))
3067 PG_RETURN_NULL();
3068
3069 /* init arg */
3070 arg = rtpg_clip_arg_init();
3071 if (arg == NULL) {
3072 elog(ERROR, "RASTER_clip: Could not initialize argument structure");
3073 PG_RETURN_NULL();
3074 }
3075
3076 /* raster (0) */
3077 pgraster = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
3078
3079 /* Get raster object */
3080 arg->raster = rt_raster_deserialize(pgraster, FALSE);
3081 if (arg->raster == NULL) {
3083 PG_FREE_IF_COPY(pgraster, 0);
3084 elog(ERROR, "RASTER_clip: Could not deserialize raster");
3085 PG_RETURN_NULL();
3086 }
3087
3088 /* raster is empty, return empty raster */
3089 if (rt_raster_is_empty(arg->raster) || rt_raster_get_num_bands(arg->raster) == 0) {
3090 elog(NOTICE, "Input raster is empty or has no bands. Returning empty raster");
3091
3093 PG_FREE_IF_COPY(pgraster, 0);
3094
3095 rtn = rt_raster_new(0, 0);
3096 if (rtn == NULL) {
3097 elog(ERROR, "RASTER_clip: Could not create empty raster");
3098 PG_RETURN_NULL();
3099 }
3100
3101 pgrtn = rt_raster_serialize(rtn);
3102 rt_raster_destroy(rtn);
3103 if (NULL == pgrtn)
3104 PG_RETURN_NULL();
3105
3106 SET_VARSIZE(pgrtn, pgrtn->size);
3107 PG_RETURN_POINTER(pgrtn);
3108 }
3109
3110 /* metadata */
3112 srid = clamp_srid(rt_raster_get_srid(arg->raster));
3113
3114 /* geometry (2) */
3115 gser = PG_GETARG_GSERIALIZED_P(2);
3116 geom = lwgeom_from_gserialized(gser);
3117
3118 /* Get a 2D version of the geometry if necessary */
3119 if (lwgeom_ndims(geom) > 2) {
3120 LWGEOM *geom2d = lwgeom_force_2d(geom);
3121 lwgeom_free(geom);
3122 geom = geom2d;
3123 }
3124
3125 /* check that SRIDs match */
3126 if (srid != clamp_srid(gserialized_get_srid(gser))) {
3127 elog(NOTICE, "Geometry provided does not have the same SRID as the raster. Returning NULL");
3128
3130 PG_FREE_IF_COPY(pgraster, 0);
3131 lwgeom_free(geom);
3132 PG_FREE_IF_COPY(gser, 2);
3133
3134 PG_RETURN_NULL();
3135 }
3136
3137 /* crop (4) */
3138 if (!PG_ARGISNULL(4) && !PG_GETARG_BOOL(4))
3139 arg->extenttype = ET_FIRST;
3140
3141 /* get intersection geometry of input raster and input geometry */
3142 if (rt_raster_get_convex_hull(arg->raster, &rastgeom) != ES_NONE) {
3143
3145 PG_FREE_IF_COPY(pgraster, 0);
3146 lwgeom_free(geom);
3147 PG_FREE_IF_COPY(gser, 2);
3148
3149 elog(ERROR, "RASTER_clip: Could not get convex hull of raster");
3150 PG_RETURN_NULL();
3151 }
3152
3153 tmpgeom = lwgeom_intersection(rastgeom, geom);
3154 lwgeom_free(rastgeom);
3155 lwgeom_free(geom);
3156 PG_FREE_IF_COPY(gser, 2);
3157 geom = tmpgeom;
3158
3159 /* intersection is empty AND extent type is INTERSECTION, return empty */
3160 if (lwgeom_is_empty(geom) && arg->extenttype == ET_INTERSECTION) {
3161 elog(NOTICE, "The input raster and input geometry do not intersect. Returning empty raster");
3162
3164 PG_FREE_IF_COPY(pgraster, 0);
3165 lwgeom_free(geom);
3166
3167 rtn = rt_raster_new(0, 0);
3168 if (rtn == NULL) {
3169 elog(ERROR, "RASTER_clip: Could not create empty raster");
3170 PG_RETURN_NULL();
3171 }
3172
3173 pgrtn = rt_raster_serialize(rtn);
3174 rt_raster_destroy(rtn);
3175 if (NULL == pgrtn)
3176 PG_RETURN_NULL();
3177
3178 SET_VARSIZE(pgrtn, pgrtn->size);
3179 PG_RETURN_POINTER(pgrtn);
3180 }
3181
3182 /* nband (1) */
3183 if (!PG_ARGISNULL(1)) {
3184 array = PG_GETARG_ARRAYTYPE_P(1);
3185 etype = ARR_ELEMTYPE(array);
3186 get_typlenbyvalalign(etype, &typlen, &typbyval, &typalign);
3187
3188 switch (etype) {
3189 case INT2OID:
3190 case INT4OID:
3191 break;
3192 default:
3194 PG_FREE_IF_COPY(pgraster, 0);
3195 lwgeom_free(geom);
3196 elog(ERROR, "RASTER_clip: Invalid data type for band indexes");
3197 PG_RETURN_NULL();
3198 break;
3199 }
3200
3201 deconstruct_array(
3202 array, etype,
3203 typlen, typbyval, typalign,
3204 &e, &nulls, &(arg->numbands)
3205 );
3206
3207 arg->band = palloc(sizeof(struct rtpg_clip_band_t) * arg->numbands);
3208 if (arg->band == NULL) {
3210 PG_FREE_IF_COPY(pgraster, 0);
3211 lwgeom_free(geom);
3212 elog(ERROR, "RASTER_clip: Could not allocate memory for band arguments");
3213 PG_RETURN_NULL();
3214 }
3215
3216 for (i = 0, j = 0; i < arg->numbands; i++) {
3217 if (nulls[i]) continue;
3218
3219 switch (etype) {
3220 case INT2OID:
3221 arg->band[j].nband = DatumGetInt16(e[i]) - 1;
3222 break;
3223 case INT4OID:
3224 arg->band[j].nband = DatumGetInt32(e[i]) - 1;
3225 break;
3226 }
3227
3228 j++;
3229 }
3230
3231 if (j < arg->numbands) {
3232 arg->band = repalloc(arg->band, sizeof(struct rtpg_clip_band_t) * j);
3233 if (arg->band == NULL) {
3235 PG_FREE_IF_COPY(pgraster, 0);
3236 lwgeom_free(geom);
3237 elog(ERROR, "RASTER_clip: Could not reallocate memory for band arguments");
3238 PG_RETURN_NULL();
3239 }
3240
3241 arg->numbands = j;
3242 }
3243
3244 /* validate band */
3245 for (i = 0; i < arg->numbands; i++) {
3246 if (!rt_raster_has_band(arg->raster, arg->band[i].nband)) {
3247 elog(NOTICE, "Band at index %d not found in raster", arg->band[i].nband + 1);
3249 PG_FREE_IF_COPY(pgraster, 0);
3250 lwgeom_free(geom);
3251 PG_RETURN_NULL();
3252 }
3253
3254 arg->band[i].hasnodata = 0;
3255 arg->band[i].nodataval = 0;
3256 }
3257 }
3258 else {
3260
3261 /* raster may have no bands */
3262 if (arg->numbands) {
3263 arg->band = palloc(sizeof(struct rtpg_clip_band_t) * arg->numbands);
3264 if (arg->band == NULL) {
3265
3267 PG_FREE_IF_COPY(pgraster, 0);
3268 lwgeom_free(geom);
3269
3270 elog(ERROR, "RASTER_clip: Could not allocate memory for band arguments");
3271 PG_RETURN_NULL();
3272 }
3273
3274 for (i = 0; i < arg->numbands; i++) {
3275 arg->band[i].nband = i;
3276 arg->band[i].hasnodata = 0;
3277 arg->band[i].nodataval = 0;
3278 }
3279 }
3280 }
3281
3282 /* nodataval (3) */
3283 if (!PG_ARGISNULL(3)) {
3284 array = PG_GETARG_ARRAYTYPE_P(3);
3285 etype = ARR_ELEMTYPE(array);
3286 get_typlenbyvalalign(etype, &typlen, &typbyval, &typalign);
3287
3288 switch (etype) {
3289 case FLOAT4OID:
3290 case FLOAT8OID:
3291 break;
3292 default:
3294 PG_FREE_IF_COPY(pgraster, 0);
3295 lwgeom_free(geom);
3296 elog(ERROR, "RASTER_clip: Invalid data type for NODATA values");
3297 PG_RETURN_NULL();
3298 break;
3299 }
3300
3301 deconstruct_array(
3302 array, etype,
3303 typlen, typbyval, typalign,
3304 &e, &nulls, &k
3305 );
3306
3307 /* it doesn't matter if there are more nodataval */
3308 for (i = 0, j = 0; i < arg->numbands; i++, j++) {
3309 /* cap j to the last nodataval element */
3310 if (j >= k)
3311 j = k - 1;
3312
3313 if (nulls[j])
3314 continue;
3315
3316 arg->band[i].hasnodata = 1;
3317 switch (etype) {
3318 case FLOAT4OID:
3319 arg->band[i].nodataval = DatumGetFloat4(e[j]);
3320 break;
3321 case FLOAT8OID:
3322 arg->band[i].nodataval = DatumGetFloat8(e[j]);
3323 break;
3324 }
3325 }
3326 }
3327
3328 /* get wkb of geometry */
3329 POSTGIS_RT_DEBUG(3, "getting wkb of geometry");
3330 wkb = lwgeom_to_wkb(geom, WKB_SFSQL, &wkb_len);
3331 lwgeom_free(geom);
3332
3333 /* rasterize geometry */
3335 wkb, wkb_len,
3336 NULL,
3337 0, NULL,
3338 NULL, NULL,
3339 NULL, NULL,
3340 NULL, NULL,
3341 &(gt[1]), &(gt[5]),
3342 NULL, NULL,
3343 &(gt[0]), &(gt[3]),
3344 &(gt[2]), &(gt[4]),
3345 NULL
3346 );
3347
3348 pfree(wkb);
3349 if (arg->mask == NULL) {
3351 PG_FREE_IF_COPY(pgraster, 0);
3352 elog(ERROR, "RASTER_clip: Could not rasterize intersection geometry");
3353 PG_RETURN_NULL();
3354 }
3355
3356 /* set SRID */
3357 rt_raster_set_srid(arg->mask, srid);
3358
3359 /* run iterator */
3360
3361 /* init itrset */
3362 itrset = palloc(sizeof(struct rt_iterator_t) * 2);
3363 if (itrset == NULL) {
3365 PG_FREE_IF_COPY(pgraster, 0);
3366 elog(ERROR, "RASTER_clip: Could not allocate memory for iterator arguments");
3367 PG_RETURN_NULL();
3368 }
3369
3370 /* one band at a time */
3371 for (i = 0; i < arg->numbands; i++) {
3372 POSTGIS_RT_DEBUGF(4, "band arg %d (nband, hasnodata, nodataval) = (%d, %d, %f)",
3373 i, arg->band[i].nband, arg->band[i].hasnodata, arg->band[i].nodataval);
3374
3375 band = rt_raster_get_band(arg->raster, arg->band[i].nband);
3376
3377 /* band metadata */
3378 pixtype = rt_band_get_pixtype(band);
3379
3380 if (arg->band[i].hasnodata) {
3381 hasnodata = 1;
3382 nodataval = arg->band[i].nodataval;
3383 }
3384 else if (rt_band_get_hasnodata_flag(band)) {
3385 hasnodata = 1;
3386 rt_band_get_nodata(band, &nodataval);
3387 }
3388 else {
3389 hasnodata = 0;
3390 nodataval = rt_band_get_min_value(band);
3391 }
3392
3393 /* band is NODATA, create NODATA band and continue */
3394 if (rt_band_get_isnodata_flag(band)) {
3395 /* create raster */
3396 if (rtn == NULL) {
3397 noerr = rt_raster_from_two_rasters(arg->raster, arg->mask, arg->extenttype, &rtn, NULL);
3398 if (noerr != ES_NONE) {
3400 PG_FREE_IF_COPY(pgraster, 0);
3401 elog(ERROR, "RASTER_clip: Could not create output raster");
3402 PG_RETURN_NULL();
3403 }
3404 }
3405
3406 /* create NODATA band */
3407 if (rt_raster_generate_new_band(rtn, pixtype, nodataval, hasnodata, nodataval, i) < 0) {
3408 rt_raster_destroy(rtn);
3410 PG_FREE_IF_COPY(pgraster, 0);
3411 elog(ERROR, "RASTER_clip: Could not add NODATA band to output raster");
3412 PG_RETURN_NULL();
3413 }
3414
3415 continue;
3416 }
3417
3418 /* raster */
3419 itrset[0].raster = arg->raster;
3420 itrset[0].nband = arg->band[i].nband;
3421 itrset[0].nbnodata = 1;
3422
3423 /* mask */
3424 itrset[1].raster = arg->mask;
3425 itrset[1].nband = 0;
3426 itrset[1].nbnodata = 1;
3427
3428 /* pass to iterator */
3429 noerr = rt_raster_iterator(
3430 itrset, 2,
3431 arg->extenttype, NULL,
3432 pixtype,
3433 hasnodata, nodataval,
3434 0, 0,
3435 NULL,
3436 NULL,
3438 &_raster
3439 );
3440
3441 if (noerr != ES_NONE) {
3442 pfree(itrset);
3444 if (rtn != NULL) rt_raster_destroy(rtn);
3445 PG_FREE_IF_COPY(pgraster, 0);
3446 elog(ERROR, "RASTER_clip: Could not run raster iterator function");
3447 PG_RETURN_NULL();
3448 }
3449
3450 /* new raster */
3451 if (rtn == NULL)
3452 rtn = _raster;
3453 /* copy band */
3454 else {
3455 band = rt_raster_get_band(_raster, 0);
3456 if (band == NULL) {
3457 pfree(itrset);
3459 rt_raster_destroy(_raster);
3460 rt_raster_destroy(rtn);
3461 PG_FREE_IF_COPY(pgraster, 0);
3462 elog(NOTICE, "RASTER_clip: Could not get band from working raster");
3463 PG_RETURN_NULL();
3464 }
3465
3466 if (rt_raster_add_band(rtn, band, i) < 0) {
3467 pfree(itrset);
3469 rt_raster_destroy(_raster);
3470 rt_raster_destroy(rtn);
3471 PG_FREE_IF_COPY(pgraster, 0);
3472 elog(ERROR, "RASTER_clip: Could not add new band to output raster");
3473 PG_RETURN_NULL();
3474 }
3475
3476 rt_raster_destroy(_raster);
3477 }
3478 }
3479
3480 pfree(itrset);
3482 PG_FREE_IF_COPY(pgraster, 0);
3483
3484 pgrtn = rt_raster_serialize(rtn);
3485 rt_raster_destroy(rtn);
3486
3487 POSTGIS_RT_DEBUG(3, "Finished");
3488
3489 if (!pgrtn)
3490 PG_RETURN_NULL();
3491
3492 SET_VARSIZE(pgrtn, pgrtn->size);
3493 PG_RETURN_POINTER(pgrtn);
3494}
#define FALSE
Definition dbfopen.c:168
int32_t gserialized_get_srid(const GSERIALIZED *g)
Extract the SRID from the serialized form (it is packed into three bytes so this is a handy function)...
LWGEOM * lwgeom_from_gserialized(const GSERIALIZED *g)
Allocate a new LWGEOM from a GSERIALIZED.
int lwgeom_ndims(const LWGEOM *geom)
Return the number of dimensions (2, 3, 4) in a geometry.
Definition lwgeom.c:937
void lwgeom_free(LWGEOM *geom)
Definition lwgeom.c:1138
LWGEOM * lwgeom_force_2d(const LWGEOM *geom)
Strip out the Z/M components of an LWGEOM.
Definition lwgeom.c:775
LWGEOM * lwgeom_intersection(const LWGEOM *geom1, const LWGEOM *geom2)
uint8_t * lwgeom_to_wkb(const LWGEOM *geom, uint8_t variant, size_t *size_out)
Convert LWGEOM to a char* in WKB format.
Definition lwout_wkb.c:790
#define SRID_UNKNOWN
Unknown SRID value.
Definition liblwgeom.h:229
#define WKB_SFSQL
Definition liblwgeom.h:2122
int32_t clamp_srid(int32_t srid)
Return a valid SRID from an arbitrary integer Raises a notice if what comes out is different from wha...
Definition lwutil.c:333
int32_t rt_raster_get_srid(rt_raster raster)
Get raster's SRID.
Definition rt_raster.c:356
int rt_raster_generate_new_band(rt_raster raster, rt_pixtype pixtype, double initialvalue, uint32_t hasnodata, double nodatavalue, int index)
Generate a new inline band and add it to a raster.
Definition rt_raster.c:485
int rt_raster_add_band(rt_raster raster, rt_band band, int index)
Add band data to a raster.
Definition rt_raster.c:405
int rt_band_get_hasnodata_flag(rt_band band)
Get hasnodata flag value.
Definition rt_band.c:674
void rt_raster_destroy(rt_raster raster)
Release memory associated to a raster.
Definition rt_raster.c:82
rt_pixtype
Definition librtcore.h:185
int rt_band_get_isnodata_flag(rt_band band)
Get isnodata flag value.
Definition rt_band.c:714
rt_raster rt_raster_new(uint32_t width, uint32_t height)
Construct a raster with given dimensions.
Definition rt_raster.c:48
int rt_raster_has_band(rt_raster raster, int nband)
Return TRUE if the raster has a band of this number.
Definition rt_raster.c:1342
double rt_band_get_min_value(rt_band band)
Returns the minimal possible value for the band according to the pixel type.
Definition rt_band.c:1745
@ ES_NONE
Definition librtcore.h:180
rt_raster rt_raster_gdal_rasterize(const unsigned char *wkb, uint32_t wkb_len, const char *srs, uint32_t num_bands, rt_pixtype *pixtype, double *init, double *value, double *nodata, uint8_t *hasnodata, int *width, int *height, double *scale_x, double *scale_y, double *ul_xw, double *ul_yw, double *grid_xw, double *grid_yw, double *skew_x, double *skew_y, char **options)
Return a raster of the provided geometry.
Definition rt_raster.c:2504
uint16_t rt_raster_get_num_bands(rt_raster raster)
Definition rt_raster.c:372
rt_errorstate rt_raster_get_convex_hull(rt_raster raster, LWGEOM **hull)
Get raster's convex hull.
rt_errorstate rt_raster_iterator(rt_iterator itrset, uint16_t itrcount, rt_extenttype extenttype, rt_raster customextent, rt_pixtype pixtype, uint8_t hasnodata, double nodataval, uint16_t distancex, uint16_t distancey, rt_mask mask, void *userarg, int(*callback)(rt_iterator_arg arg, void *userarg, double *value, int *nodata), rt_raster *rtnraster)
n-raster iterator.
void rt_raster_set_srid(rt_raster raster, int32_t srid)
Set raster's SRID.
Definition rt_raster.c:363
@ ET_INTERSECTION
Definition librtcore.h:201
@ ET_FIRST
Definition librtcore.h:203
rt_errorstate rt_band_get_nodata(rt_band band, double *nodata)
Get NODATA value.
Definition rt_band.c:1730
rt_pixtype rt_band_get_pixtype(rt_band band)
Return pixeltype of this band.
Definition rt_band.c:631
void * rt_raster_serialize(rt_raster raster)
Return this raster in serialized form.
rt_errorstate rt_raster_from_two_rasters(rt_raster rast1, rt_raster rast2, rt_extenttype extenttype, rt_raster *rtnraster, double *offset)
Definition rt_raster.c:3350
void rt_raster_get_geotransform_matrix(rt_raster raster, double *gt)
Get 6-element array of raster geotransform matrix.
Definition rt_raster.c:706
rt_raster rt_raster_deserialize(void *serialized, int header_only)
Return a raster from a serialized form.
int rt_raster_is_empty(rt_raster raster)
Return TRUE if the raster is empty.
Definition rt_raster.c:1329
rt_band rt_raster_get_band(rt_raster raster, int bandNum)
Return Nth band, or NULL if unavailable.
Definition rt_raster.c:381
static int lwgeom_is_empty(const LWGEOM *geom)
Return true or false depending on whether a geometry is an "empty" geometry (no vertices members)
Definition lwinline.h:193
static rtpg_clip_arg rtpg_clip_arg_init()
static void rtpg_clip_arg_destroy(rtpg_clip_arg arg)
static int rtpg_clip_callback(rt_iterator_arg arg, void *userarg, double *value, int *nodata)
#define POSTGIS_RT_DEBUG(level, msg)
Definition rtpostgis.h:61
#define POSTGIS_RT_DEBUGF(level, msg,...)
Definition rtpostgis.h:65
rt_raster raster
Definition librtcore.h:2444
uint16_t nband
Definition librtcore.h:2445
uint8_t nbnodata
Definition librtcore.h:2446
Struct definitions.
Definition librtcore.h:2251
rtpg_clip_band band
rt_extenttype extenttype

References rtpg_clip_arg_t::band, clamp_srid(), ES_NONE, ET_FIRST, ET_INTERSECTION, rtpg_clip_arg_t::extenttype, FALSE, gserialized_get_srid(), 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: