[SCM] PostGIS branch master updated. 3.6.0rc2-495-gc9bdc50fb
git at osgeo.org
git at osgeo.org
Thu Jun 4 09:33:47 PDT 2026
This is an automated email from the git hooks/post-receive script. It was
generated because a ref change was pushed to the repository containing
the project "PostGIS".
The branch, master has been updated
via c9bdc50fbdaac675ef784463df1849e27d447872 (commit)
from 7761433a55ce7445282798f018d57819f44dc319 (commit)
Those revisions listed above that are new to this repository have
not appeared on any other notification email; so we list those
revisions in full, below.
- Log -----------------------------------------------------------------
commit c9bdc50fbdaac675ef784463df1849e27d447872
Author: Paul Ramsey <pramsey at cleverelephant.ca>
Date: Thu Jun 4 09:32:40 2026 -0700
Change target target reskew calculation from an iterative
process to a deterministic calculation. Closes #5854.
References #5918.
diff --git a/raster/rt_core/librtcore.h b/raster/rt_core/librtcore.h
index 6847319cf..618d1292d 100644
--- a/raster/rt_core/librtcore.h
+++ b/raster/rt_core/librtcore.h
@@ -1544,21 +1544,17 @@ rt_errorstate rt_raster_get_perimeter(
/*
* Compute skewed extent that covers unskewed extent.
*
- * @param envelope : unskewed extent of type rt_envelope
- * @param skew : pointer to 2-element array (x, y) of skew
- * @param scale : pointer to 2-element array (x, y) of scale
- * @param tolerance : value between 0 and 1 where the smaller the tolerance
- * results in an extent approaching the "minimum" skewed extent.
- * If value <= 0, tolerance = 0.1. If value > 1, tolerance = 1.
+ * @param extent : unskewed extent of type rt_envelope
+ * @param skew : pointer to 2-element array {skew_x, skew_y}
+ * @param scale : pointer to 2-element array {fabs(scale_x), fabs(scale_y)}
*
- * @return skewed raster who's extent covers unskewed extent, NULL on error
+ * @return skewed raster whose extent covers unskewed extent, NULL on error
*/
rt_raster
rt_raster_compute_skewed_raster(
rt_envelope extent,
double *skew,
- double *scale,
- double tolerance
+ double *scale
);
/**
diff --git a/raster/rt_core/rt_raster.c b/raster/rt_core/rt_raster.c
index 99578fe87..3ce611261 100644
--- a/raster/rt_core/rt_raster.c
+++ b/raster/rt_core/rt_raster.c
@@ -854,381 +854,96 @@ rt_raster_get_envelope(
******************************************************************************/
/*
- * Compute skewed extent that covers unskewed extent.
+ * Compute skewed raster extent that covers the given unskewed extent.
+ * Uses a direct linear-algebra solution rather than an iterative approach.
*
- * @param envelope : unskewed extent of type rt_envelope
- * @param skew : pointer to 2-element array (x, y) of skew
- * @param scale : pointer to 2-element array (x, y) of scale
- * @param tolerance : value between 0 and 1 where the smaller the tolerance
- * results in an extent approaching the "minimum" skewed extent.
- * If value <= 0, tolerance = 0.1. If value > 1, tolerance = 1.
+ * @param extent : unskewed extent of type rt_envelope
+ * @param skew : pointer to 2-element array {skew_x, skew_y}
+ * @param scale : pointer to 2-element array {fabs(scale_x), fabs(scale_y)}
*
- * @return skewed raster who's extent covers unskewed extent, NULL on error
+ * @return skewed raster whose extent covers unskewed extent, NULL on error
*/
rt_raster
rt_raster_compute_skewed_raster(
rt_envelope extent,
double *skew,
- double *scale,
- double tolerance
+ double *scale
) {
- uint32_t run = 0;
- uint32_t max_run = 1;
- double dbl_run = 0;
-
- int rtn;
- int covers = 0;
+ double sx, syg, kx, ky, det;
+ double col_nums[4], row_nums[4];
+ double col_min, col_max, row_min, row_max;
+ double K1, K2, eps;
+ double ulx, uly;
+ int width, height;
rt_raster raster;
- double _gt[6] = {0};
- double _igt[6] = {0};
- int _d[2] = {1, -1};
- int _dlast = 0;
- int _dlastpos = 0;
- double _w[2] = {0};
- double _r[2] = {0};
- double _xy[2] = {0};
int i;
- int j;
- int x;
- int y;
- LWGEOM *geom = NULL;
- GEOSGeometry *sgeom = NULL;
- GEOSGeometry *ngeom = NULL;
-
- if (
- (tolerance < 0.) ||
- FLT_EQ(tolerance, 0.)
- ) {
- tolerance = 0.1;
- }
- else if (tolerance > 1.)
- tolerance = 1;
-
- dbl_run = tolerance;
- while (dbl_run < 10) {
- dbl_run *= 10.;
- max_run *= 10;
- }
-
- /* scale must be provided */
- if (scale == NULL)
+ if (scale == NULL || skew == NULL)
return NULL;
- for (i = 0; i < 2; i++) {
- if (FLT_EQ(scale[i], 0.0))
- {
- rterror("rt_raster_compute_skewed_raster: Scale cannot be zero");
- return 0;
- }
- if (i < 1)
- _gt[1] = fabs(scale[i] * tolerance);
- else
- _gt[5] = fabs(scale[i] * tolerance);
- }
- /* conform scale-y to be negative */
- _gt[5] *= -1;
-
- /* skew not provided or skew is zero, return raster of correct dim and spatial attributes */
- if ((skew == NULL) || (FLT_EQ(skew[0], 0.0) && FLT_EQ(skew[1], 0.0)))
- {
- int _dim[2] = {
- (int) fmax((fabs(extent.MaxX - extent.MinX) + (fabs(scale[0]) / 2.)) / fabs(scale[0]), 1),
- (int) fmax((fabs(extent.MaxY - extent.MinY) + (fabs(scale[1]) / 2.)) / fabs(scale[1]), 1)
- };
-
- raster = rt_raster_new(_dim[0], _dim[1]);
- if (raster == NULL) {
- rterror("rt_raster_compute_skewed_raster: Could not create output raster");
- return NULL;
- }
-
- rt_raster_set_offsets(raster, extent.MinX, extent.MaxY);
- rt_raster_set_scale(raster, fabs(scale[0]), -1 * fabs(scale[1]));
- rt_raster_set_skews(raster, skew[0], skew[1]);
-
- return raster;
- }
-
- /* direction to shift upper-left corner */
- if (skew[0] > 0.)
- _d[0] = -1;
- if (skew[1] < 0.)
- _d[1] = 1;
-
- /* geotransform */
- _gt[0] = extent.UpperLeftX;
- _gt[2] = skew[0] * tolerance;
- _gt[3] = extent.UpperLeftY;
- _gt[4] = skew[1] * tolerance;
-
- RASTER_DEBUGF(4, "Initial geotransform: %f, %f, %f, %f, %f, %f",
- _gt[0], _gt[1], _gt[2], _gt[3], _gt[4], _gt[5]
- );
- RASTER_DEBUGF(4, "Delta: %d, %d", _d[0], _d[1]);
-
- /* simple raster */
- if ((raster = rt_raster_new(1, 1)) == NULL) {
- rterror("rt_raster_compute_skewed_raster: Out of memory allocating extent raster");
- return NULL;
- }
- rt_raster_set_geotransform_matrix(raster, _gt);
-
- /* get inverse geotransform matrix */
- if (!GDALInvGeoTransform(_gt, _igt)) {
- rterror("rt_raster_compute_skewed_raster: Could not compute inverse geotransform matrix");
- rt_raster_destroy(raster);
- return NULL;
- }
- RASTER_DEBUGF(4, "Inverse geotransform: %f, %f, %f, %f, %f, %f",
- _igt[0], _igt[1], _igt[2], _igt[3], _igt[4], _igt[5]
- );
-
- /* shift along axis */
- for (i = 0; i < 2; i++) {
- covers = 0;
- run = 0;
-
- RASTER_DEBUGF(3, "Shifting along %s axis", i < 1 ? "X" : "Y");
-
- do {
-
- /* prevent possible infinite loop */
- if (run > max_run) {
- rterror("rt_raster_compute_skewed_raster: Could not compute skewed extent due to check preventing infinite loop");
- rt_raster_destroy(raster);
- return NULL;
- }
-
- /*
- check the four corners that they are covered along the specific axis
- pixel column should be >= 0
- */
- for (j = 0; j < 4; j++) {
- switch (j) {
- /* upper-left */
- case 0:
- _xy[0] = extent.MinX;
- _xy[1] = extent.MaxY;
- break;
- /* lower-left */
- case 1:
- _xy[0] = extent.MinX;
- _xy[1] = extent.MinY;
- break;
- /* lower-right */
- case 2:
- _xy[0] = extent.MaxX;
- _xy[1] = extent.MinY;
- break;
- /* upper-right */
- case 3:
- _xy[0] = extent.MaxX;
- _xy[1] = extent.MaxY;
- break;
- }
-
- rtn = rt_raster_geopoint_to_cell(
- raster,
- _xy[0], _xy[1],
- &(_r[0]), &(_r[1]),
- _igt
- );
- if (rtn != ES_NONE) {
- rterror("rt_raster_compute_skewed_raster: Could not compute raster pixel for spatial coordinates");
- rt_raster_destroy(raster);
- return NULL;
- }
-
- RASTER_DEBUGF(4, "Point %d at cell %d x %d", j, (int) _r[0], (int) _r[1]);
-
- /* raster doesn't cover point */
- if ((int) _r[i] < 0) {
- RASTER_DEBUGF(4, "Point outside of skewed extent: %d", j);
- covers = 0;
-
- if (_dlastpos != j) {
- _dlast = (int) _r[i];
- _dlastpos = j;
- }
- else if ((int) _r[i] < _dlast) {
- RASTER_DEBUG(4, "Point going in wrong direction. Reversing direction");
- _d[i] *= -1;
- _dlastpos = -1;
- run = 0;
- }
-
- break;
- }
-
- covers++;
- }
-
- if (!covers) {
- x = 0;
- y = 0;
- if (i < 1)
- x = _d[i] * fabs(_r[i]);
- else
- y = _d[i] * fabs(_r[i]);
-
- rtn = rt_raster_cell_to_geopoint(
- raster,
- x, y,
- &(_w[0]), &(_w[1]),
- _gt
- );
- if (rtn != ES_NONE) {
- rterror("rt_raster_compute_skewed_raster: Could not compute spatial coordinates for raster pixel");
- rt_raster_destroy(raster);
- return NULL;
- }
-
- /* adjust ul */
- if (i < 1)
- _gt[0] = _w[i];
- else
- _gt[3] = _w[i];
- rt_raster_set_geotransform_matrix(raster, _gt);
- RASTER_DEBUGF(4, "Shifted geotransform: %f, %f, %f, %f, %f, %f",
- _gt[0], _gt[1], _gt[2], _gt[3], _gt[4], _gt[5]
- );
-
- /* get inverse geotransform matrix */
- if (!GDALInvGeoTransform(_gt, _igt)) {
- rterror("rt_raster_compute_skewed_raster: Could not compute inverse geotransform matrix");
- rt_raster_destroy(raster);
- return NULL;
- }
- RASTER_DEBUGF(4, "Inverse geotransform: %f, %f, %f, %f, %f, %f",
- _igt[0], _igt[1], _igt[2], _igt[3], _igt[4], _igt[5]
- );
- }
-
- run++;
- }
- while (!covers);
- }
-
- /* covers test */
- rtn = rt_raster_geopoint_to_cell(
- raster,
- extent.MaxX, extent.MinY,
- &(_r[0]), &(_r[1]),
- _igt
- );
- if (rtn != ES_NONE) {
- rterror("rt_raster_compute_skewed_raster: Could not compute raster pixel for spatial coordinates");
- rt_raster_destroy(raster);
+ if (FLT_EQ(scale[0], 0.0) || FLT_EQ(scale[1], 0.0)) {
+ rterror("rt_raster_compute_skewed_raster: Scale cannot be zero");
return NULL;
}
- RASTER_DEBUGF(4, "geopoint %f x %f at cell %d x %d", extent.MaxX, extent.MinY, (int) _r[0], (int) _r[1]);
+ sx = fabs(scale[0]);
+ syg = -fabs(scale[1]); /* gt[5]: negative for north-up */
+ kx = skew[0];
+ ky = skew[1];
+ det = sx * syg - kx * ky;
- raster->width = _r[0];
- raster->height = _r[1];
-
- /* initialize GEOS */
- initGEOS(rtinfo, lwgeom_geos_error);
-
- /* create reference LWPOLY */
- {
- LWPOLY *npoly = rt_util_envelope_to_lwpoly(extent);
- if (npoly == NULL) {
- rterror("rt_raster_compute_skewed_raster: Could not build extent's geometry for covers test");
- rt_raster_destroy(raster);
- return NULL;
- }
-
- ngeom = (GEOSGeometry *) LWGEOM2GEOS(lwpoly_as_lwgeom(npoly), 0);
- lwpoly_free(npoly);
+ if (FLT_EQ(det, 0.0)) {
+ rterror("rt_raster_compute_skewed_raster: Degenerate geotransform (determinant is zero)");
+ return NULL;
}
- do {
- covers = 0;
-
- /* construct sgeom from raster */
- if ((rt_raster_get_convex_hull(raster, &geom) != ES_NONE) || geom == NULL) {
- rterror("rt_raster_compute_skewed_raster: Could not build skewed extent's geometry for covers test");
- GEOSGeom_destroy(ngeom);
- rt_raster_destroy(raster);
- return NULL;
- }
-
- sgeom = (GEOSGeometry *) LWGEOM2GEOS(geom, 0);
- lwgeom_free(geom);
-
- covers = GEOSRelatePattern(sgeom, ngeom, "******FF*");
- GEOSGeom_destroy(sgeom);
-
- if (covers == 2) {
- rterror("rt_raster_compute_skewed_raster: Could not run covers test");
- GEOSGeom_destroy(ngeom);
- rt_raster_destroy(raster);
- return NULL;
- }
-
- if (!covers)
- {
- raster->width++;
- raster->height++;
- }
- }
- while (!covers);
-
- RASTER_DEBUGF(4, "Skewed extent does cover normal extent with dimensions %d x %d", raster->width, raster->height);
-
- raster->width = (int) ((((double) raster->width) * fabs(_gt[1]) + fabs(scale[0] / 2.)) / fabs(scale[0]));
- raster->height = (int) ((((double) raster->height) * fabs(_gt[5]) + fabs(scale[1] / 2.)) / fabs(scale[1]));
- _gt[1] = fabs(scale[0]);
- _gt[5] = -1 * fabs(scale[1]);
- _gt[2] = skew[0];
- _gt[4] = skew[1];
- rt_raster_set_geotransform_matrix(raster, _gt);
-
- /* minimize width/height */
- for (i = 0; i < 2; i++) {
- covers = 1;
- do {
- if (i < 1)
- raster->width--;
- else
- raster->height--;
-
- /* construct sgeom from raster */
- if ((rt_raster_get_convex_hull(raster, &geom) != ES_NONE) || geom == NULL) {
- rterror("rt_raster_compute_skewed_raster: Could not build skewed extent's geometry for minimizing dimensions");
- GEOSGeom_destroy(ngeom);
- rt_raster_destroy(raster);
- return NULL;
- }
-
- sgeom = (GEOSGeometry *) LWGEOM2GEOS(geom, 0);
- lwgeom_free(geom);
-
- covers = GEOSRelatePattern(sgeom, ngeom, "******FF*");
- GEOSGeom_destroy(sgeom);
-
- if (covers == 2) {
- rterror("rt_raster_compute_skewed_raster: Could not run covers test for minimizing dimensions");
- GEOSGeom_destroy(ngeom);
- rt_raster_destroy(raster);
- return NULL;
- }
- } while (covers);
-
- if (i < 1)
- raster->width++;
- else
- raster->height++;
+ /* col_num(x,y) = syg*x - kx*y, row_num(x,y) = -ky*x + sx*y */
+ col_nums[0] = syg * extent.MinX - kx * extent.MaxY;
+ col_nums[1] = syg * extent.MaxX - kx * extent.MaxY;
+ col_nums[2] = syg * extent.MaxX - kx * extent.MinY;
+ col_nums[3] = syg * extent.MinX - kx * extent.MinY;
+ row_nums[0] = -ky * extent.MinX + sx * extent.MaxY;
+ row_nums[1] = -ky * extent.MaxX + sx * extent.MaxY;
+ row_nums[2] = -ky * extent.MaxX + sx * extent.MinY;
+ row_nums[3] = -ky * extent.MinX + sx * extent.MinY;
+ col_min = col_max = col_nums[0];
+ row_min = row_max = row_nums[0];
+ for (i = 1; i < 4; i++) {
+ if (col_nums[i] < col_min) col_min = col_nums[i];
+ if (col_nums[i] > col_max) col_max = col_nums[i];
+ if (row_nums[i] < row_min) row_min = row_nums[i];
+ if (row_nums[i] > row_max) row_max = row_nums[i];
}
- GEOSGeom_destroy(ngeom);
+ K1 = (det < 0.0) ? col_max : col_min;
+ K2 = (det < 0.0) ? row_max : row_min;
+
+ /* Tiny safety margin so FP rounding can't place boundary corners just outside */
+ eps = (fabs(K1) + fabs(K2) + fabs(det)) * DBL_EPSILON * 16.0;
+ if (det < 0.0) { K1 += eps; K2 += eps; }
+ else { K1 -= eps; K2 -= eps; }
+
+ /* Solve: syg*ulx - kx*uly = K1, -ky*ulx + sx*uly = K2 */
+ ulx = (sx * K1 + kx * K2) / det;
+ uly = (ky * K1 + syg * K2) / det;
+
+ width = (int)floor((col_max - col_min) / fabs(det)) + 1;
+ height = (int)floor((row_max - row_min) / fabs(det)) + 1;
+
+ raster = rt_raster_new(width, height);
+ if (raster == NULL) {
+ rterror("rt_raster_compute_skewed_raster: Could not create output raster");
+ return NULL;
+ }
+ rt_raster_set_offsets(raster, ulx, uly);
+ rt_raster_set_scale(raster, fabs(scale[0]), -fabs(scale[1]));
+ rt_raster_set_skews(raster, skew[0], skew[1]);
return raster;
}
+
/**
* Return TRUE if the raster is empty. i.e. is NULL, width = 0 or height = 0
*
@@ -2749,8 +2464,7 @@ rt_raster_gdal_rasterize(
skewedrast = rt_raster_compute_skewed_raster(
extent,
_skew,
- _scale,
- 0.01
+ _scale
);
if (skewedrast == NULL) {
rterror("rt_raster_gdal_rasterize: Could not compute skewed raster");
diff --git a/raster/rt_core/rt_warp.c b/raster/rt_core/rt_warp.c
index 59110feec..f5854fd0f 100644
--- a/raster/rt_core/rt_warp.c
+++ b/raster/rt_core/rt_warp.c
@@ -145,34 +145,9 @@ _rti_warp_arg_destroy(_rti_warp_arg arg) {
}
/**
- * Return a warped raster using GDAL Warp API
- *
- * @param raster : raster to transform
- * @param src_srs : the raster's coordinate system in OGC WKT
- * @param dst_srs : the warped raster's coordinate system in OGC WKT
- * @param scale_x : the x size of pixels of the warped raster's pixels in
- * units of dst_srs
- * @param scale_y : the y size of pixels of the warped raster's pixels in
- * units of dst_srs
- * @param width : the number of columns of the warped raster. note that
- * width/height CANNOT be used with scale_x/scale_y
- * @param height : the number of rows of the warped raster. note that
- * width/height CANNOT be used with scale_x/scale_y
- * @param ul_xw : the X value of upper-left corner of the warped raster in
- * units of dst_srs
- * @param ul_yw : the Y value of upper-left corner of the warped raster in
- * units of dst_srs
- * @param grid_xw : the X value of point on a grid to align warped raster
- * to in units of dst_srs
- * @param grid_yw : the Y value of point on a grid to align warped raster
- * to in units of dst_srs
- * @param skew_x : the X skew of the warped raster in units of dst_srs
- * @param skew_y : the Y skew of the warped raster in units of dst_srs
- * @param resample_alg : the resampling algorithm
- * @param max_err : maximum error measured in input pixels permitted
- * (0.0 for exact calculations)
- *
- * @return the warped raster or NULL
+ * Return a warped raster using GDAL Warp API.
+ * When non-zero skew values are requested, the output extent is computed
+ * via rt_raster_compute_skewed_raster().
*/
rt_raster rt_raster_gdal_warp(
rt_raster raster,
@@ -207,103 +182,64 @@ rt_raster rt_raster_gdal_warp(
int i = 0;
int numBands = 0;
- /* flag indicating that the spatial info is being substituted */
int subspatial = 0;
- RASTER_DEBUG(3, "starting");
+ RASTER_DEBUG(3, "rt_raster_gdal_warp: starting");
assert(NULL != raster);
- /* internal variables */
arg = _rti_warp_arg_init();
if (arg == NULL) {
rterror("rt_raster_gdal_warp: Could not initialize internal variables");
return NULL;
}
- /*
- max_err must be gte zero
-
- the value 0.125 is the default used in gdalwarp.cpp on line 283
- */
if (max_err < 0.) max_err = 0.125;
- RASTER_DEBUGF(4, "max_err = %f", max_err);
/* handle srs */
if (src_srs != NULL) {
- /* reprojection taking place */
if (dst_srs != NULL && strcmp(src_srs, dst_srs) != 0) {
- RASTER_DEBUG(4, "Warp operation does include a reprojection");
+ RASTER_DEBUG(4, "rt_raster_gdal_warp: Warp includes reprojection");
arg->src.srs = rt_util_gdal_convert_sr(src_srs, 0);
arg->dst.srs = rt_util_gdal_convert_sr(dst_srs, 0);
-
if (arg->src.srs == NULL || arg->dst.srs == NULL) {
rterror("rt_raster_gdal_warp: Could not convert srs values to GDAL accepted format");
_rti_warp_arg_destroy(arg);
return NULL;
}
}
- /* no reprojection, a stub just for clarity */
- else {
- RASTER_DEBUG(4, "Warp operation does NOT include reprojection");
- }
}
else if (dst_srs != NULL) {
- /* dst_srs provided but not src_srs */
rterror("rt_raster_gdal_warp: SRS required for input raster if SRS provided for warped raster");
_rti_warp_arg_destroy(arg);
return NULL;
}
- /* load raster into a GDAL MEM dataset */
+ /* load raster into GDAL MEM dataset */
arg->src.ds = rt_raster_to_gdal_mem(raster, arg->src.srs, NULL, NULL, 0, &(arg->src.drv), &(arg->src.destroy_drv));
if (NULL == arg->src.ds) {
rterror("rt_raster_gdal_warp: Could not convert raster to GDAL MEM format");
_rti_warp_arg_destroy(arg);
return NULL;
}
- RASTER_DEBUG(3, "raster loaded into GDAL MEM dataset");
- /* special case when src_srs and dst_srs is NULL and raster's geotransform matrix is default */
+ /* special case: default geotransform raster with no SRS */
if (
src_srs == NULL && dst_srs == NULL &&
rt_raster_get_srid(raster) == SRID_UNKNOWN
) {
double gt[6];
-
-#if POSTGIS_DEBUG_LEVEL > 3
- GDALGetGeoTransform(arg->src.ds, gt);
- RASTER_DEBUGF(3, "GDAL MEM geotransform: %f, %f, %f, %f, %f, %f",
- gt[0], gt[1], gt[2], gt[3], gt[4], gt[5]);
-#endif
-
- /* default geotransform */
rt_raster_get_geotransform_matrix(raster, gt);
- RASTER_DEBUGF(3, "raster geotransform: %f, %f, %f, %f, %f, %f",
- gt[0], gt[1], gt[2], gt[3], gt[4], gt[5]);
-
- /* substitute spatial info (lack of) with a real one EPSG:32731 (WGS84/UTM zone 31s) */
if (FLT_EQ(gt[0], 0.0) && FLT_EQ(gt[3], 0.0) && FLT_EQ(gt[1], 1.0) && FLT_EQ(gt[5], -1.0) &&
FLT_EQ(gt[2], 0.0) && FLT_EQ(gt[4], 0.0))
{
double ngt[6] = {166021.4431, 0.1, 0, 10000000.0000, 0, -0.1};
-
rtwarn("Raster has default geotransform. Adjusting metadata for use of GDAL Warp API");
-
subspatial = 1;
-
GDALSetGeoTransform(arg->src.ds, ngt);
GDALFlushCache(arg->src.ds);
-
- /* EPSG:32731 */
arg->src.srs = rt_util_gdal_convert_sr("EPSG:32731", 0);
arg->dst.srs = rt_util_gdal_convert_sr("EPSG:32731", 0);
-
-#if POSTGIS_DEBUG_LEVEL > 3
- GDALGetGeoTransform(arg->src.ds, gt);
- RASTER_DEBUGF(3, "GDAL MEM geotransform: %f, %f, %f, %f, %f, %f",
- gt[0], gt[1], gt[2], gt[3], gt[4], gt[5]);
-#endif
}
}
@@ -312,7 +248,7 @@ rt_raster rt_raster_gdal_warp(
arg->transform.option.len = 2;
arg->transform.option.item = rtalloc(sizeof(char *) * (arg->transform.option.len + 1));
if (NULL == arg->transform.option.item) {
- rterror("rt_raster_gdal_warp: Could not allocation memory for transform options");
+ rterror("rt_raster_gdal_warp: Could not allocate memory for transform options");
_rti_warp_arg_destroy(arg);
return NULL;
}
@@ -322,15 +258,14 @@ rt_raster rt_raster_gdal_warp(
const char *srs = i ? arg->dst.srs : arg->src.srs;
const char *lbl = i ? "DST_SRS=" : "SRC_SRS=";
size_t sz = sizeof(char) * (strlen(lbl) + 1);
- if ( srs ) sz += strlen(srs);
+ if (srs) sz += strlen(srs);
arg->transform.option.item[i] = (char *) rtalloc(sz);
if (NULL == arg->transform.option.item[i]) {
- rterror("rt_raster_gdal_warp: Could not allocation memory for transform options");
+ rterror("rt_raster_gdal_warp: Could not allocate memory for transform options");
_rti_warp_arg_destroy(arg);
return NULL;
}
sprintf(arg->transform.option.item[i], "%s%s", lbl, srs ? srs : "");
- RASTER_DEBUGF(4, "arg->transform.option.item[%d] = %s", i, arg->transform.option.item[i]);
}
}
else
@@ -339,7 +274,7 @@ rt_raster rt_raster_gdal_warp(
/* transformation object for building dst dataset */
arg->transform.arg.transform = GDALCreateGenImgProjTransformer2(arg->src.ds, NULL, arg->transform.option.item);
if (NULL == arg->transform.arg.transform) {
- rterror("rt_raster_gdal_warp: Could not create GDAL transformation object for output dataset creation");
+ rterror("rt_raster_gdal_warp: Could not create GDAL transformation object");
_rti_warp_arg_destroy(arg);
return NULL;
}
@@ -349,70 +284,49 @@ rt_raster rt_raster_gdal_warp(
arg->src.ds, GDALGenImgProjTransform,
arg->transform.arg.transform, _gt, &(_dim[0]), &(_dim[1]), dst_extent, 0);
if (cplerr != CE_None) {
- rterror("rt_raster_gdal_warp: Could not get GDAL suggested warp output for output dataset creation");
+ rterror("rt_raster_gdal_warp: Could not get GDAL suggested warp output");
_rti_warp_arg_destroy(arg);
return NULL;
}
GDALDestroyGenImgProjTransformer(arg->transform.arg.transform);
arg->transform.arg.transform = NULL;
- /*
- don't use suggested dimensions as use of suggested scales
- on suggested extent will result in suggested dimensions
- */
_dim[0] = 0;
_dim[1] = 0;
- RASTER_DEBUGF(3, "Suggested geotransform: %f, %f, %f, %f, %f, %f",
- _gt[0], _gt[1], _gt[2], _gt[3], _gt[4], _gt[5]);
-
- /* store extent in easier-to-use object */
extent.MinX = dst_extent[0];
extent.MinY = dst_extent[1];
extent.MaxX = dst_extent[2];
extent.MaxY = dst_extent[3];
-
extent.UpperLeftX = dst_extent[0];
extent.UpperLeftY = dst_extent[3];
- RASTER_DEBUGF(3, "Suggested extent: %f, %f, %f, %f",
- extent.MinX, extent.MinY, extent.MaxX, extent.MaxY);
-
/* scale and width/height are mutually exclusive */
if (
((NULL != scale_x) || (NULL != scale_y)) &&
((NULL != width) || (NULL != height))
) {
- rterror("rt_raster_gdal_warp: Scale X/Y and width/height are mutually exclusive. Only provide one");
+ rterror("rt_raster_gdal_warp: Scale X/Y and width/height are mutually exclusive");
_rti_warp_arg_destroy(arg);
return NULL;
}
- /* user-defined width */
if (NULL != width) {
_dim[0] = abs(*width);
_scale[0] = fabs((extent.MaxX - extent.MinX) / ((double) _dim[0]));
}
- /* user-defined height */
if (NULL != height) {
_dim[1] = abs(*height);
_scale[1] = fabs((extent.MaxY - extent.MinY) / ((double) _dim[1]));
}
- /* user-defined scale */
if (
((NULL != scale_x) && (FLT_NEQ(*scale_x, 0.0))) &&
((NULL != scale_y) && (FLT_NEQ(*scale_y, 0.0)))
) {
_scale[0] = fabs(*scale_x);
_scale[1] = fabs(*scale_y);
-
- /* special override since we changed the original GT scales */
if (subspatial) {
- /*
- _scale[0] *= 10;
- _scale[1] *= 10;
- */
_scale[0] /= 10;
_scale[1] /= 10;
}
@@ -421,83 +335,59 @@ rt_raster rt_raster_gdal_warp(
((NULL != scale_x) && (NULL == scale_y)) ||
((NULL == scale_x) && (NULL != scale_y))
) {
- rterror("rt_raster_gdal_warp: Both X and Y scale values must be provided for scale");
+ rterror("rt_raster_gdal_warp: Both X and Y scale values must be provided");
_rti_warp_arg_destroy(arg);
return NULL;
}
- /* scale not defined, use suggested */
- if (FLT_EQ(_scale[0], 0.0) && FLT_EQ(_scale[1], 0.0))
- {
+ /* scale not defined: use suggested */
+ if (FLT_EQ(_scale[0], 0.0) && FLT_EQ(_scale[1], 0.0)) {
_scale[0] = fabs(_gt[1]);
_scale[1] = fabs(_gt[5]);
}
- RASTER_DEBUGF(4, "Using scale: %f x %f", _scale[0], -1 * _scale[1]);
+ RASTER_DEBUGF(4, "rt_raster_gdal_warp: Using scale: %f x %f", _scale[0], -1 * _scale[1]);
- /* user-defined skew */
+ /* skew */
if (NULL != skew_x) {
_skew[0] = *skew_x;
-
- /*
- negative scale-x affects skew
- for now, force skew to be in left-right, top-down orientation
- */
- if (
- NULL != scale_x &&
- *scale_x < 0.
- ) {
+ if (NULL != scale_x && *scale_x < 0.)
_skew[0] *= -1;
- }
}
if (NULL != skew_y) {
_skew[1] = *skew_y;
-
- /*
- positive scale-y affects skew
- for now, force skew to be in left-right, top-down orientation
- */
- if (
- NULL != scale_y &&
- *scale_y > 0.
- ) {
+ if (NULL != scale_y && *scale_y > 0)
_skew[1] *= -1;
- }
}
- RASTER_DEBUGF(4, "Using skew: %f x %f", _skew[0], _skew[1]);
+ RASTER_DEBUGF(4, "rt_raster_gdal_warp: Using skew: %f x %f", _skew[0], _skew[1]);
- /* reprocess extent if skewed */
- if (FLT_NEQ(_skew[0], 0.0) || FLT_NEQ(_skew[1], 0.0))
- {
+ /* Compute the output grid parameters when skew is requested. */
+ if (FLT_NEQ(_skew[0], 0.0) || FLT_NEQ(_skew[1], 0.0)) {
rt_raster skewedrast;
- RASTER_DEBUG(3, "Computing skewed extent's envelope");
+ RASTER_DEBUG(3, "rt_raster_gdal_warp: Computing skewed extent");
- skewedrast = rt_raster_compute_skewed_raster(
- extent,
- _skew,
- _scale,
- 0.01
- );
+ skewedrast = rt_raster_compute_skewed_raster(extent, _skew, _scale);
if (skewedrast == NULL) {
- rterror("rt_raster_gdal_warp: Could not compute skewed raster");
+ rterror("rt_raster_gdal_warp: Could not compute skewed extent");
_rti_warp_arg_destroy(arg);
return NULL;
}
- if (_dim[0] == 0)
- _dim[0] = skewedrast->width;
- if (_dim[1] == 0)
- _dim[1] = skewedrast->height;
+ if (_dim[0] == 0) _dim[0] = rt_raster_get_width(skewedrast);
+ if (_dim[1] == 0) _dim[1] = rt_raster_get_height(skewedrast);
- extent.UpperLeftX = skewedrast->ipX;
- extent.UpperLeftY = skewedrast->ipY;
+ extent.UpperLeftX = rt_raster_get_x_offset(skewedrast);
+ extent.UpperLeftY = rt_raster_get_y_offset(skewedrast);
+
+ RASTER_DEBUGF(3, "rt_raster_gdal_warp: Skewed extent: UL=(%f,%f) dim=%dx%d",
+ extent.UpperLeftX, extent.UpperLeftY, _dim[0], _dim[1]);
rt_raster_destroy(skewedrast);
}
- /* dimensions not defined, compute */
+ /* dimensions not defined, compute from extent */
if (!_dim[0])
_dim[0] = (int) fmax((fabs(extent.MaxX - extent.MinX) + (_scale[0] / 2.)) / _scale[0], 1);
if (!_dim[1])
@@ -517,21 +407,12 @@ rt_raster rt_raster_gdal_warp(
rt_raster_set_skews(rast, _skew[0], _skew[1]);
rt_raster_get_geotransform_matrix(rast, _gt);
- RASTER_DEBUGF(3, "Temp raster's geotransform: %f, %f, %f, %f, %f, %f",
+ RASTER_DEBUGF(3, "rt_raster_gdal_warp: Temp raster geotransform: %f,%f,%f,%f,%f,%f",
_gt[0], _gt[1], _gt[2], _gt[3], _gt[4], _gt[5]);
- RASTER_DEBUGF(3, "Temp raster's dimensions (width x height): %d x %d",
- _dim[0], _dim[1]);
/* user-defined upper-left corner */
- if (
- NULL != ul_xw &&
- NULL != ul_yw
- ) {
+ if (NULL != ul_xw && NULL != ul_yw) {
ul_user = 1;
-
- RASTER_DEBUGF(4, "Using user-specified upper-left corner: %f, %f", *ul_xw, *ul_yw);
-
- /* set upper-left corner */
rt_raster_set_offsets(rast, *ul_xw, *ul_yw);
extent.UpperLeftX = *ul_xw;
extent.UpperLeftY = *ul_yw;
@@ -599,7 +480,6 @@ rt_raster rt_raster_gdal_warp(
NULL
) != ES_NONE) {
rterror("rt_raster_gdal_warp: Could not compute spatial coordinates for raster pixel");
-
rt_raster_destroy(rast);
_rti_warp_arg_destroy(arg);
return NULL;
@@ -646,7 +526,6 @@ rt_raster rt_raster_gdal_warp(
NULL
) != ES_NONE) {
rterror("rt_raster_gdal_warp: Could not compute spatial coordinates for raster pixel");
-
rt_raster_destroy(rast);
_rti_warp_arg_destroy(arg);
return NULL;
@@ -662,68 +541,35 @@ rt_raster rt_raster_gdal_warp(
while (0);
}
- /*
- after this point, rt_envelope extent is no longer used
- */
-
- /* get key attributes from rast */
+ /* get key attributes */
_dim[0] = rast->width;
_dim[1] = rast->height;
rt_raster_get_geotransform_matrix(rast, _gt);
- /* scale-x is negative or scale-y is positive */
- if ((
- (NULL != scale_x) && (*scale_x < 0.)
- ) || (
- (NULL != scale_y) && (*scale_y > 0)
- )) {
+ /* handle negative scale-x or positive scale-y (same as original) */
+ if ((NULL != scale_x && *scale_x < 0.) || (NULL != scale_y && *scale_y > 0)) {
double _w[2] = {0};
-
- /* negative scale-x */
- if (
- (NULL != scale_x) &&
- (*scale_x < 0.)
- ) {
- if (rt_raster_cell_to_geopoint(
- rast,
- rast->width, 0,
- &(_w[0]), &(_w[1]),
- NULL
- ) != ES_NONE) {
+ if (NULL != scale_x && *scale_x < 0.) {
+ if (rt_raster_cell_to_geopoint(rast, rast->width, 0, &(_w[0]), &(_w[1]), NULL) != ES_NONE) {
rterror("rt_raster_gdal_warp: Could not compute spatial coordinates for raster pixel");
rt_raster_destroy(rast);
_rti_warp_arg_destroy(arg);
return NULL;
}
-
_gt[0] = _w[0];
_gt[1] = *scale_x;
-
- /* check for skew */
if (NULL != skew_x && FLT_NEQ(*skew_x, 0.0))
_gt[2] = *skew_x;
}
- /* positive scale-y */
- if (
- (NULL != scale_y) &&
- (*scale_y > 0)
- ) {
- if (rt_raster_cell_to_geopoint(
- rast,
- 0, rast->height,
- &(_w[0]), &(_w[1]),
- NULL
- ) != ES_NONE) {
+ if (NULL != scale_y && *scale_y > 0) {
+ if (rt_raster_cell_to_geopoint(rast, 0, rast->height, &(_w[0]), &(_w[1]), NULL) != ES_NONE) {
rterror("rt_raster_gdal_warp: Could not compute spatial coordinates for raster pixel");
rt_raster_destroy(rast);
_rti_warp_arg_destroy(arg);
return NULL;
}
-
_gt[3] = _w[1];
_gt[5] = *scale_y;
-
- /* check for skew */
if (NULL != skew_y && FLT_NEQ(*skew_y, 0.0))
_gt[4] = *skew_y;
}
@@ -732,12 +578,11 @@ rt_raster rt_raster_gdal_warp(
rt_raster_destroy(rast);
rast = NULL;
- RASTER_DEBUGF(3, "Applied geotransform: %f, %f, %f, %f, %f, %f",
+ RASTER_DEBUGF(3, "rt_raster_gdal_warp: Final geotransform: %f,%f,%f,%f,%f,%f",
_gt[0], _gt[1], _gt[2], _gt[3], _gt[4], _gt[5]);
- RASTER_DEBUGF(3, "Raster dimensions (width x height): %d x %d",
- _dim[0], _dim[1]);
+ RASTER_DEBUGF(3, "rt_raster_gdal_warp: Raster dimensions: %d x %d", _dim[0], _dim[1]);
- if ( _dim[0] == 0 || _dim[1] == 0 ) {
+ if (_dim[0] == 0 || _dim[1] == 0) {
rterror("rt_raster_gdal_warp: The width (%d) or height (%d) of the warped raster is zero", _dim[0], _dim[1]);
_rti_warp_arg_destroy(arg);
return NULL;
@@ -745,7 +590,6 @@ rt_raster rt_raster_gdal_warp(
/* load VRT driver */
if (!rt_util_gdal_driver_registered("VRT")) {
- RASTER_DEBUG(3, "Registering VRT driver");
GDALRegister_VRT();
arg->dst.destroy_drv = 1;
}
@@ -772,7 +616,6 @@ rt_raster rt_raster_gdal_warp(
_rti_warp_arg_destroy(arg);
return NULL;
}
- RASTER_DEBUGF(3, "Applied SRS: %s", GDALGetProjectionRef(arg->dst.ds));
}
/* set dst geotransform */
@@ -783,7 +626,7 @@ rt_raster rt_raster_gdal_warp(
return NULL;
}
- /* add bands to dst dataset */
+ /* add bands */
numBands = rt_raster_get_num_bands(raster);
for (i = 0; i < numBands; i++) {
rtband = rt_raster_get_band(raster, i);
@@ -805,8 +648,6 @@ rt_raster rt_raster_gdal_warp(
return NULL;
}
- /* get band to write data to */
- band = NULL;
band = GDALGetRasterBand(arg->dst.ds, i + 1);
if (NULL == band) {
rterror("rt_raster_gdal_warp: Could not get GDAL band for additional processing");
@@ -814,12 +655,10 @@ rt_raster rt_raster_gdal_warp(
return NULL;
}
- /* set nodata */
if (rt_band_get_hasnodata_flag(rtband) != FALSE) {
rt_band_get_nodata(rtband, &nodata);
if (GDALSetRasterNoDataValue(band, nodata) != CE_None)
rtwarn("rt_raster_gdal_warp: Could not set nodata value for band %d", i);
- RASTER_DEBUGF(3, "nodata value set to %f", GDALGetRasterNoDataValue(band, NULL));
}
}
@@ -835,7 +674,6 @@ rt_raster rt_raster_gdal_warp(
}
arg->transform.func = GDALGenImgProjTransform;
- /* use approximate transformation object */
if (max_err > 0.0) {
arg->transform.arg.transform = arg->transform.arg.approx = GDALCreateApproxTransformer(
GDALGenImgProjTransform,
@@ -846,7 +684,6 @@ rt_raster rt_raster_gdal_warp(
_rti_warp_arg_destroy(arg);
return NULL;
}
-
arg->transform.func = GDALApproxTransform;
}
@@ -858,37 +695,25 @@ rt_raster rt_raster_gdal_warp(
return NULL;
}
- /* set options */
arg->wopts->eResampleAlg = resample_alg;
arg->wopts->hSrcDS = arg->src.ds;
arg->wopts->hDstDS = arg->dst.ds;
arg->wopts->pfnTransformer = arg->transform.func;
arg->wopts->pTransformerArg = arg->transform.arg.transform;
arg->wopts->papszWarpOptions = (char **) CPLMalloc(sizeof(char *) * 2);
-
arg->wopts->papszWarpOptions[0] = (char *) CPLMalloc(sizeof(char) * (strlen("INIT_DEST=NO_DATA") + 1));
strcpy(arg->wopts->papszWarpOptions[0], "INIT_DEST=NO_DATA");
arg->wopts->papszWarpOptions[1] = NULL;
- /* set band mapping */
+ /* band mapping */
arg->wopts->nBandCount = numBands;
arg->wopts->panSrcBands = (int *) CPLMalloc(sizeof(int) * arg->wopts->nBandCount);
arg->wopts->panDstBands = (int *) CPLMalloc(sizeof(int) * arg->wopts->nBandCount);
for (i = 0; i < arg->wopts->nBandCount; i++)
arg->wopts->panDstBands[i] = arg->wopts->panSrcBands[i] = i + 1;
- /*
- * https://trac.osgeo.org/postgis/ticket/5881
- * In order to call GDALWarp with BAND_INIT=NO_DATA we need to ensure
- * that the src and dst rasters have nodata values and they are
- * matched up nicely. This block used by tested with the hasnodata
- * check on all src raster bands, but now we just do it every time
- * because that makes sense (any warped raster is likely to have
- * empty corners on output, and those corners need to be filled with
- * some kind of NODATA value).
- */
+ /* nodata mapping */
{
- RASTER_DEBUG(3, "Setting nodata mapping");
arg->wopts->padfSrcNoDataReal = (double *) CPLMalloc(numBands * sizeof(double));
arg->wopts->padfDstNoDataReal = (double *) CPLMalloc(numBands * sizeof(double));
arg->wopts->padfSrcNoDataImag = (double *) CPLMalloc(numBands * sizeof(double));
@@ -910,68 +735,43 @@ rt_raster rt_raster_gdal_warp(
_rti_warp_arg_destroy(arg);
return NULL;
}
-
- if (!rt_band_get_hasnodata_flag(band)) {
- /*
- based on line 1004 of gdalwarp.cpp
- the problem is that there is a chance that this number is a legitimate value
- */
+ if (!rt_band_get_hasnodata_flag(band))
arg->wopts->padfSrcNoDataReal[i] = -123456.789;
- }
- else {
+ else
rt_band_get_nodata(band, &(arg->wopts->padfSrcNoDataReal[i]));
- }
-
arg->wopts->padfDstNoDataReal[i] = arg->wopts->padfSrcNoDataReal[i];
arg->wopts->padfDstNoDataImag[i] = arg->wopts->padfSrcNoDataImag[i] = 0.0;
- RASTER_DEBUGF(4, "Mapped nodata value for band %d: %f (%f) => %f (%f)",
- i,
- arg->wopts->padfSrcNoDataReal[i], arg->wopts->padfSrcNoDataImag[i],
- arg->wopts->padfDstNoDataReal[i], arg->wopts->padfDstNoDataImag[i]
- );
}
}
- /* warp raster */
- RASTER_DEBUG(3, "Warping raster");
+ /* warp */
cplerr = GDALInitializeWarpedVRT(arg->dst.ds, arg->wopts);
if (cplerr != CE_None) {
rterror("rt_raster_gdal_warp: Could not warp raster");
_rti_warp_arg_destroy(arg);
return NULL;
}
- /*
- GDALSetDescription(arg->dst.ds, "/tmp/warped.vrt");
- */
GDALFlushCache(arg->dst.ds);
- RASTER_DEBUG(3, "Raster warped");
- /* convert gdal dataset to raster */
- RASTER_DEBUG(3, "Converting GDAL dataset to raster");
+ /* convert to raster */
rast = rt_raster_from_gdal_dataset(arg->dst.ds);
-
_rti_warp_arg_destroy(arg);
if (NULL == rast) {
- rterror("rt_raster_gdal_warp: Could not warp raster");
+ rterror("rt_raster_gdal_warp: Could not convert warped dataset to raster");
return NULL;
}
- /* substitute spatial, reset back to default */
+ /* reset spatial info for default-geotransform rasters */
if (subspatial) {
double gt[6] = {0, 1, 0, 0, 0, -1};
- /* See http://trac.osgeo.org/postgis/ticket/2911 */
- /* We should probably also tweak rotation here */
- /* NOTE: the times 10 is because it was divided by 10 in a section above,
- * I'm not sure the above division was needed */
gt[1] = _scale[0] * 10;
gt[5] = -1 * _scale[1] * 10;
-
rt_raster_set_geotransform_matrix(rast, gt);
rt_raster_set_srid(rast, SRID_UNKNOWN);
}
- RASTER_DEBUG(3, "done");
+ RASTER_DEBUG(3, "rt_raster_gdal_warp: done");
return rast;
}
diff --git a/raster/test/cunit/cu_raster_misc.c b/raster/test/cunit/cu_raster_misc.c
index 90888a9aa..048749193 100644
--- a/raster/test/cunit/cu_raster_misc.c
+++ b/raster/test/cunit/cu_raster_misc.c
@@ -191,17 +191,14 @@ static void test_raster_compute_skewed_raster(void) {
rast = rt_raster_compute_skewed_raster(
extent,
skew,
- scale,
- 0
+ scale
);
CU_ASSERT(rast != NULL);
- /* Check disabled: See https://trac.osgeo.org/postgis/ticket/4379
- * CU_ASSERT_EQUAL(rt_raster_get_width(rast), 2);
- */
+ CU_ASSERT_EQUAL(rt_raster_get_width(rast), 3);
CU_ASSERT_EQUAL(rt_raster_get_height(rast), 3);
- CU_ASSERT_DOUBLE_EQUAL(rt_raster_get_x_offset(rast), -0.5, DBL_EPSILON);
- CU_ASSERT_DOUBLE_EQUAL(rt_raster_get_y_offset(rast), 0, DBL_EPSILON);
+ CU_ASSERT_DOUBLE_EQUAL(rt_raster_get_x_offset(rast), -8.0/17.0, 1e-8);
+ CU_ASSERT_DOUBLE_EQUAL(rt_raster_get_y_offset(rast), -2.0/17.0, 1e-8);
cu_free_raster(rast);
}
diff --git a/raster/test/regress/rt_asraster_expected b/raster/test/regress/rt_asraster_expected
index 1e8fef17d..50110cee8 100644
--- a/raster/test/regress/rt_asraster_expected
+++ b/raster/test/regress/rt_asraster_expected
@@ -46,10 +46,10 @@ NOTICE: The rasters have different scales on the X axis
3.0||||||||||||||||
3.1|993310|100|100|1|1406.537|-869.114|0.000|0.000|-175453.086|114987.661|8BUI|0.000|t|255.000|255.000|
3.2|993310|100|100|1|1406.537|-869.114|0.000|0.000|-175453.086|114987.661|8BUI|0.000|t|255.000|255.000|
-3.3|993310|100|100|1|1406.537|-869.114|1.000|0.000|-175565.609|114987.661|8BUI|0.000|t|255.000|255.000|
-3.4|993310|100|100|1|1406.537|-869.114|0.000|1.000|-175453.086|114987.661|8BUI|0.000|t|255.000|255.000|
-3.5|993310|101|101|1|1406.537|-869.114|10.000|-5.000|-176465.793|115491.747|8BUI|0.000|t|255.000|255.000|
-3.6|993310|100|101|1|1406.537|-869.114|-5.000|10.000|-175453.086|114987.661|8BUI|0.000|t|255.000|255.000|
+3.3|993310|101|101|1|1406.537|-869.114|1.000|0.000|-175553.086|114987.661|8BUI|0.000|t|255.000|255.000|
+3.4|993310|100|101|1|1406.537|-869.114|0.000|1.000|-175453.086|114987.661|8BUI|0.000|t|255.000|255.000|
+3.5|993310|101|101|1|1406.537|-869.114|10.000|-5.000|-176458.881|115491.236|8BUI|0.000|t|255.000|255.000|
+3.6|993310|101|102|1|1406.537|-869.114|-5.000|10.000|-175453.086|114987.661|8BUI|0.000|t|255.000|255.000|
4.0||||||||||||||||
4.1|992163|150|117|1|1000.000|-1000.000|0.000|0.000|-1898000.000|-412000.000|8BUI|0.000|t|1.000|1.000|t
4.10|993310|142|88|1|1000.000|-1000.000|0.000|0.000|-176100.000|115100.000|16BUI|0.000|t|13.000|13.000|f
diff --git a/raster/test/regress/rt_gdalwarp_expected b/raster/test/regress/rt_gdalwarp_expected
index 649f8d98e..d6bc47c1f 100644
--- a/raster/test/regress/rt_gdalwarp_expected
+++ b/raster/test/regress/rt_gdalwarp_expected
@@ -14,7 +14,7 @@ NOTICE: Values must be provided for both X and Y when specifying the scale. Re
0.13|992163|11|11|1|1000.000|-1000.000|0.000|0.000|-500100.000|600950.000|t|t|t
0.14|992163|201|201|1|50.000|50.000|0.000|0.000|-500040.000|589957.000|t|t|t
0.15|992163|84|84|1|121.000|121.000|0.000|0.000|-500093.000|589875.000|t|t|t
-0.18|992163|10|10|1|1000.000|-1000.000|3.000|3.000|-500030.000|600000.000|t|t|t
+0.18|992163|11|11|1|1000.000|-1000.000|3.000|3.000|-500030.000|599999.910|t|t|t
0.25|0|5|5|1|2.000|-2.000|0.000|0.000|0.000|0.000|t|t|t
0.26|0|2|2|1|5.000|-5.000|0.000|0.000|0.000|0.000|t|t|t
0.27|0|100|100|1|0.100|-0.100|0.000|0.000|0.000|0.000|t|t|t
@@ -31,14 +31,14 @@ NOTICE: Values must be provided for both X and Y when specifying the scale. Re
1.13|992163|84|84|1|121.000|121.000|0.000|0.000|-500093.000|589875.000|t|t|t
1.14|992163|201|201|1|50.000|50.000|0.000|0.000|-500040.000|589957.000|t|t|t
1.15|992163|201|201|1|50.000|50.000|0.000|0.000|-500040.000|589957.000|t|t|t
-1.16|992163|10|10|1|1000.000|-1000.000|3.000|3.000|-500030.000|600000.000|t|t|t
-1.17|992163|10|10|1|1000.000|-1000.000|3.000|3.000|-500030.000|600000.000|t|t|t
-1.18|992163|10|10|1|1000.000|-1000.000|1.000|3.000|-500010.000|600000.000|t|t|t
-1.19|992163|20|20|1|500.000|500.000|3.000|3.000|-500065.000|590065.000|t|t|t
+1.16|992163|11|11|1|1000.000|-1000.000|3.000|3.000|-500030.000|599999.910|t|t|t
+1.17|992163|11|11|1|1000.000|-1000.000|3.000|3.000|-500030.000|599999.910|t|t|t
+1.18|992163|11|11|1|1000.000|-1000.000|1.000|3.000|-500010.000|599999.970|t|t|t
+1.19|992163|21|21|1|500.000|500.000|3.000|3.000|-500060.362|589560.362|t|t|t
1.2|992163|20|20|1|500.000|500.000|0.000|0.000|-500000.000|590000.000|t|t|t
-1.20|992163|21|21|1|500.000|500.000|0.000|6.000|-500048.000|590038.000|t|t|t
-1.21|992163|207|101|1|50.000|-100.000|3.000|0.000|-500319.000|600056.000|t|t|t
-1.22|992163|207|101|1|50.000|-100.000|3.000|0.000|-500319.000|600056.000|t|t|t
+1.20|992163|22|22|1|500.000|500.000|0.000|6.000|-500048.000|589538.000|t|t|t
+1.21|992163|208|102|1|50.000|-100.000|3.000|0.000|-500319.000|600056.000|t|t|t
+1.22|992163|208|102|1|50.000|-100.000|3.000|0.000|-500319.000|600056.000|t|t|t
1.23|992163|150|150|1|66.667|-66.667|0.000|0.000|-500000.000|600000.000|t|t|t
1.24|992163|5|5|1|2064.200|-2291.200|0.000|0.000|-500321.000|601456.000|t|t|t
1.25|||||||||||||
@@ -62,11 +62,11 @@ NOTICE: Values must be provided for both X and Y when specifying the scale. Re
3.7|992163|67|80|1|150.000|125.000|0.000|0.000|-500000.000|590000.000|t|t|t
3.8|992163|67|80|1|150.000|125.000|0.000|0.000|-500000.000|590000.000|t|t|t
4.1|992163|10|10|1|1000.000|-1000.000|0.000|0.000|-500000.000|600000.000|t|t|t
-4.2|992163|10|10|1|1000.000|-1000.000|1.000|1.000|-500010.000|600000.000|t|t|t
-4.3|992163|10|10|1|1000.000|-1000.000|0.500|0.000|-500010.000|600000.000|t|t|t
-4.4|992163|10|10|1|1000.000|-1000.000|10.000|10.000|-500100.000|600000.000|t|t|t
-4.5|992163|10|10|1|1000.000|-1000.000|10.000|10.000|-500100.000|600000.000|t|t|t
-4.6|992163|10|10|1|1000.000|-1000.000|10.000|10.000|-500100.000|600000.000|t|t|t
+4.2|992163|11|11|1|1000.000|-1000.000|1.000|1.000|-500010.000|599999.990|t|t|t
+4.3|992163|11|11|1|1000.000|-1000.000|0.500|0.000|-500005.000|600000.000|t|t|t
+4.4|992163|11|11|1|1000.000|-1000.000|10.000|10.000|-500099.990|599999.000|t|t|t
+4.5|992163|11|11|1|1000.000|-1000.000|10.000|10.000|-500099.990|599999.000|t|t|t
+4.6|992163|11|11|1|1000.000|-1000.000|10.000|10.000|-500099.990|599999.000|t|t|t
5.1|992163|10|10|1|1000.000|-1000.000|0.000|0.000|-500000.000|600000.000|t|t|t
5.10|992163|11|11|1|1000.000|-1000.000|0.000|0.000|-500005.000|600001.000|t|t|t
5.11|992163|11|11|1|1000.000|-1000.000|0.000|0.000|-500991.000|600991.000|t|t|t
diff --git a/raster/test/regress/tickets.sql b/raster/test/regress/tickets.sql
index a22717160..7a58d86e8 100644
--- a/raster/test/regress/tickets.sql
+++ b/raster/test/regress/tickets.sql
@@ -18,8 +18,8 @@ SELECT '#2532.2', NULL::geometry @ null::raster;
-- added OFFSET 0 to force PostgreSQL 12+ to materialize the cte
WITH data AS ( SELECT '#2911' l, ST_Metadata(ST_Rescale(
ST_AddBand(
- ST_MakeEmptyRaster(10, 10, 0, 0, 1, -1, 0, 0, 0),
- 1, '8BUI', 0, 0
+ ST_MakeEmptyRaster(10, 10, 0, 0, 1, -1, 0, 0, 0),
+ 1, '8BUI', 0, 0
),
2.0,
-2.0
@@ -37,38 +37,38 @@ DROP TABLE IF EXISTS test_raster_scale_big;
DROP TABLE IF EXISTS test_raster_scale_small;
CREATE TABLE test_raster_scale_regular (
- rid integer,
- rast raster
+ rid integer,
+ rast raster
);
CREATE TABLE test_raster_scale_big (
- rid integer,
- rast raster
+ rid integer,
+ rast raster
);
CREATE TABLE test_raster_scale_small (
- rid integer,
- rast raster
+ rid integer,
+ rast raster
);
CREATE OR REPLACE FUNCTION make_test_raster(
- table_suffix text,
+ table_suffix text,
rid integer,
- scale_x double precision,
- scale_y double precision DEFAULT 1.0
+ scale_x double precision,
+ scale_y double precision DEFAULT 1.0
)
RETURNS void
AS $$
DECLARE
rast raster;
- width integer := 2;
- height integer := 2;
- ul_x double precision := 0;
- ul_y double precision := 0;
- skew_x double precision := 0;
- skew_y double precision := 0;
- initvalue double precision := 1;
- nodataval double precision := 0;
+ width integer := 2;
+ height integer := 2;
+ ul_x double precision := 0;
+ ul_y double precision := 0;
+ skew_x double precision := 0;
+ skew_y double precision := 0;
+ initvalue double precision := 1;
+ nodataval double precision := 0;
BEGIN
rast := ST_MakeEmptyRaster(width, height, ul_x, ul_y, scale_x, scale_y, skew_x, skew_y, 0);
rast := ST_AddBand(rast, 1, '8BUI', initvalue, nodataval);
@@ -83,7 +83,7 @@ SELECT make_test_raster('regular', 1, 1.0000001);
SELECT make_test_raster('regular', 2, 0.9999999);
SELECT AddRasterConstraints('test_raster_scale_regular'::name, 'rast'::name, 'scale_x', 'scale_y');
SELECT r_table_name, r_raster_column, scale_x, scale_y FROM raster_columns
- WHERE r_raster_column = 'rast' AND r_table_name = 'test_raster_scale_regular';
+ WHERE r_raster_column = 'rast' AND r_table_name = 'test_raster_scale_regular';
-- Issues enforce_scaley_rast constraint violation
SELECT make_test_raster('regular', 3, 1.001, 0.9999999);
@@ -91,7 +91,7 @@ SELECT make_test_raster('regular', 3, 1.001, 0.9999999);
SELECT make_test_raster('big', 0, -1234567890123456789.0);
SELECT AddRasterConstraints('test_raster_scale_big'::name, 'rast'::name, 'scale_x', 'scale_y');
SELECT r_table_name, r_raster_column, scale_x, scale_y FROM raster_columns
- WHERE r_raster_column = 'rast' AND r_table_name = 'test_raster_scale_big';
+ WHERE r_raster_column = 'rast' AND r_table_name = 'test_raster_scale_big';
-- Issues enforce_scalex_rast constraint violation
SELECT make_test_raster('big', 1, -12345678901234567890.0);
@@ -101,7 +101,7 @@ SELECT make_test_raster('small', 1, 0.000011);
SELECT make_test_raster('small', 2, 0.00000999);
SELECT AddRasterConstraints('test_raster_scale_small'::name, 'rast'::name, 'scale_x', 'scale_y');
SELECT r_table_name, r_raster_column, scale_x, scale_y FROM raster_columns
- WHERE r_raster_column = 'rast' AND r_table_name = 'test_raster_scale_small';
+ WHERE r_raster_column = 'rast' AND r_table_name = 'test_raster_scale_small';
-- Issues enforce_scaley_rast constraint violation
SELECT make_test_raster('small', 3, 0.00001, 1.00001);
@@ -153,10 +153,12 @@ FROM (VALUES ('A0006', 300), ('A0006', 302)) t(a,b);
SELECT '#4770.b',
ST_Union(NULL::raster) OVER (PARTITION BY a ORDER BY b)
FROM (VALUES ('A0006', 300),
- ('A0006', 302)) t(a, b);
+ ('A0006', 302)) t(a, b);
SELECT '#4724.a', ST_SummaryStatsAgg(NULL::raster, NULL::int4, NULL::bool)
OVER (ORDER BY q) FROM generate_series(1,2) AS e(q);
SELECT '#4724.b', ST_SummaryStatsAgg(NULL::raster, NULL::int4, NULL::bool)
FROM generate_series(1,2) AS e(q);
+
+SELECT '#5854',Round(ST_Rotation(ST_Reskew(ST_AddBand(ST_MakeEmptyRaster(100, 430, 0, 0, 0.001, -0.001, 0, 0, 4269), '8BUI'::text, 1, 0), 0.0015))::numeric,3);
diff --git a/raster/test/regress/tickets_expected b/raster/test/regress/tickets_expected
index 46e1297be..1087ddadc 100644
--- a/raster/test/regress/tickets_expected
+++ b/raster/test/regress/tickets_expected
@@ -77,3 +77,4 @@ NOTICE: Adding maximum extent constraint
#4724.a|(0,,,,,)
#4724.a|(0,,,,,)
#4724.b|(0,,,,,)
+#5854|-0.983
-----------------------------------------------------------------------
Summary of changes:
raster/rt_core/librtcore.h | 14 +-
raster/rt_core/rt_raster.c | 404 +++++--------------------------
raster/rt_core/rt_warp.c | 312 +++++-------------------
raster/test/cunit/cu_raster_misc.c | 11 +-
raster/test/regress/rt_asraster_expected | 8 +-
raster/test/regress/rt_gdalwarp_expected | 26 +-
raster/test/regress/tickets.sql | 48 ++--
raster/test/regress/tickets_expected | 1 +
8 files changed, 167 insertions(+), 657 deletions(-)
hooks/post-receive
--
PostGIS
More information about the postgis-tickets
mailing list