3026{
3031
3034
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;
3055
3059 int hasnodata;
3060 double nodataval;
3061 int noerr = 0;
3062
3064
3065
3066 if (PG_ARGISNULL(0) || PG_ARGISNULL(2))
3067 PG_RETURN_NULL();
3068
3069
3071 if (arg == NULL) {
3072 elog(ERROR, "RASTER_clip: Could not initialize argument structure");
3073 PG_RETURN_NULL();
3074 }
3075
3076
3077 pgraster = (
rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
3078
3079
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
3090 elog(NOTICE, "Input raster is empty or has no bands. Returning empty raster");
3091
3093 PG_FREE_IF_COPY(pgraster, 0);
3094
3096 if (rtn == NULL) {
3097 elog(ERROR, "RASTER_clip: Could not create empty raster");
3098 PG_RETURN_NULL();
3099 }
3100
3103 if (NULL == pgrtn)
3104 PG_RETURN_NULL();
3105
3106 SET_VARSIZE(pgrtn, pgrtn->
size);
3107 PG_RETURN_POINTER(pgrtn);
3108 }
3109
3110
3113
3114
3115 gser = PG_GETARG_GSERIALIZED_P(2);
3117
3118
3122 geom = geom2d;
3123 }
3124
3125
3127 elog(NOTICE, "Geometry provided does not have the same SRID as the raster. Returning NULL");
3128
3130 PG_FREE_IF_COPY(pgraster, 0);
3132 PG_FREE_IF_COPY(gser, 2);
3133
3134 PG_RETURN_NULL();
3135 }
3136
3137
3138 if (!PG_ARGISNULL(4) && !PG_GETARG_BOOL(4))
3140
3141
3143
3145 PG_FREE_IF_COPY(pgraster, 0);
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
3156 PG_FREE_IF_COPY(gser, 2);
3157 geom = tmpgeom;
3158
3159
3161 elog(NOTICE, "The input raster and input geometry do not intersect. Returning empty raster");
3162
3164 PG_FREE_IF_COPY(pgraster, 0);
3166
3168 if (rtn == NULL) {
3169 elog(ERROR, "RASTER_clip: Could not create empty raster");
3170 PG_RETURN_NULL();
3171 }
3172
3175 if (NULL == pgrtn)
3176 PG_RETURN_NULL();
3177
3178 SET_VARSIZE(pgrtn, pgrtn->
size);
3179 PG_RETURN_POINTER(pgrtn);
3180 }
3181
3182
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);
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,
3205 );
3206
3208 if (arg->
band == NULL) {
3210 PG_FREE_IF_COPY(pgraster, 0);
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) {
3233 if (arg->
band == NULL) {
3235 PG_FREE_IF_COPY(pgraster, 0);
3237 elog(ERROR, "RASTER_clip: Could not reallocate memory for band arguments");
3238 PG_RETURN_NULL();
3239 }
3240
3242 }
3243
3244
3245 for (i = 0; i < arg->
numbands; i++) {
3247 elog(NOTICE,
"Band at index %d not found in raster", arg->
band[i].
nband + 1);
3249 PG_FREE_IF_COPY(pgraster, 0);
3251 PG_RETURN_NULL();
3252 }
3253
3256 }
3257 }
3258 else {
3260
3261
3264 if (arg->
band == NULL) {
3265
3267 PG_FREE_IF_COPY(pgraster, 0);
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++) {
3278 }
3279 }
3280 }
3281
3282
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);
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
3308 for (i = 0, j = 0; i < arg->
numbands; i++, j++) {
3309
3310 if (j >= k)
3311 j = k - 1;
3312
3313 if (nulls[j])
3314 continue;
3315
3317 switch (etype) {
3318 case FLOAT4OID:
3320 break;
3321 case FLOAT8OID:
3323 break;
3324 }
3325 }
3326 }
3327
3328
3332
3333
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
3358
3359
3360
3361
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
3371 for (i = 0; i < arg->
numbands; i++) {
3372 POSTGIS_RT_DEBUGF(4,
"band arg %d (nband, hasnodata, nodataval) = (%d, %d, %f)",
3374
3376
3377
3379
3381 hasnodata = 1;
3383 }
3385 hasnodata = 1;
3387 }
3388 else {
3389 hasnodata = 0;
3391 }
3392
3393
3395
3396 if (rtn == NULL) {
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
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
3422
3423
3425 itrset[1].
nband = 0;
3427
3428
3430 itrset, 2,
3432 pixtype,
3433 hasnodata, nodataval,
3434 0, 0,
3435 NULL,
3436 NULL,
3438 &_raster
3439 );
3440
3442 pfree(itrset);
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
3451 if (rtn == NULL)
3452 rtn = _raster;
3453
3454 else {
3456 if (band == NULL) {
3457 pfree(itrset);
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
3467 pfree(itrset);
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
3477 }
3478 }
3479
3480 pfree(itrset);
3482 PG_FREE_IF_COPY(pgraster, 0);
3483
3486
3488
3489 if (!pgrtn)
3490 PG_RETURN_NULL();
3491
3492 SET_VARSIZE(pgrtn, pgrtn->
size);
3493 PG_RETURN_POINTER(pgrtn);
3494}
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.
void lwgeom_free(LWGEOM *geom)
LWGEOM * lwgeom_force_2d(const LWGEOM *geom)
Strip out the Z/M components of an LWGEOM.
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.
#define SRID_UNKNOWN
Unknown SRID value.
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...
int32_t rt_raster_get_srid(rt_raster raster)
Get raster's SRID.
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.
int rt_raster_add_band(rt_raster raster, rt_band band, int index)
Add band data to a raster.
int rt_band_get_hasnodata_flag(rt_band band)
Get hasnodata flag value.
void rt_raster_destroy(rt_raster raster)
Release memory associated to a raster.
int rt_band_get_isnodata_flag(rt_band band)
Get isnodata flag value.
rt_raster rt_raster_new(uint32_t width, uint32_t height)
Construct a raster with given dimensions.
int rt_raster_has_band(rt_raster raster, int nband)
Return TRUE if the raster has a band of this number.
double rt_band_get_min_value(rt_band band)
Returns the minimal possible value for the band according to the pixel type.
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.
uint16_t rt_raster_get_num_bands(rt_raster raster)
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.
rt_errorstate rt_band_get_nodata(rt_band band, double *nodata)
Get NODATA value.
rt_pixtype rt_band_get_pixtype(rt_band band)
Return pixeltype of this band.
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)
void rt_raster_get_geotransform_matrix(rt_raster raster, double *gt)
Get 6-element array of raster geotransform matrix.
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.
rt_band rt_raster_get_band(rt_raster raster, int bandNum)
Return Nth band, or NULL if unavailable.
static int lwgeom_is_empty(const LWGEOM *geom)
Return true or false depending on whether a geometry is an "empty" geometry (no vertices members)
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)
#define POSTGIS_RT_DEBUGF(level, msg,...)