* $Id$
*
* WKTRaster - Raster Types for PostGIS
- * http://www.postgis.org/support/wiki/index.php?WKTRasterHomePage
+ * http://trac.osgeo.org/postgis/wiki/WKTRaster
*
* Copyright (C) 2011 Regents of the University of California
* <bkpark@ucdavis.edu>
*
* @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
- * @param scale_x : the pixel width of the warped raster
- * @param scale_y : the pixel height of the warped raster
- * @param ul_xw : the X value of upper-left corner of the warped raster
- * @param ul_yw : the Y value of upper-left corner of the warped raster
- * @param grid_xw : the X value of point on a grid to align warped raster to
- * @param grid_yw : the Y value of point on a grid to align warped raster to
- * @param skew_x : the X skew of the warped raster
- * @param skew_y : the Y skew of the warped raster
+ * @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)
rt_raster raster, const char *src_srs,
const char *dst_srs,
double *scale_x, double *scale_y,
+ int *width, int *height,
double *ul_xw, double *ul_yw,
double *grid_xw, double *grid_yw,
double *skew_x, double *skew_y,
double dst_gt[6] = {0};
double dst_extent[4];
- int width = 0;
- int height = 0;
+ int _width = 0;
+ int _height = 0;
int ul_user = 0;
double min_x = 0;
double min_y = 0;
/* get approximate output georeferenced bounds and resolution */
cplerr = GDALSuggestedWarpOutput2(src_ds, GDALGenImgProjTransform,
- transform_arg, dst_gt, &width, &height, dst_extent, 0);
+ transform_arg, dst_gt, &_width, &_height, dst_extent, 0);
GDALDestroyGenImgProjTransformer(transform_arg);
if (cplerr != CE_None) {
rterror("rt_raster_gdal_warp: Unable to get GDAL suggested warp output for output dataset creation\n");
if (NULL != skew_y)
dst_gt[4] = *skew_y;
+ /* 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.\n");
+
+ GDALClose(src_ds);
+
+ for (i = 0; i < transform_opts_len; i++) rtdealloc(transform_opts[i]);
+ rtdealloc(transform_opts);
+
+ GDALDeregisterDriver(src_drv);
+ GDALDestroyDriver(src_drv);
+
+ return NULL;
+ }
+
+ /* user-defined width/height */
+ if ((NULL != width) && (*width > 0.)) {
+ _width = *width;
+ dst_gt[1] = (dst_extent[2] - dst_extent[0]) / ((double) _width);
+ pix_x = 0;
+ }
+ if ((NULL != height) && (*height > 0.)) {
+ _height = *height;
+ dst_gt[5] = -1 * fabs((dst_extent[3] - dst_extent[1]) / ((double) _height));
+ pix_y = 0;
+ }
+
/* user-defined scale */
if (
- (NULL != scale_x) &&
- (FLT_NEQ(*scale_x, 0.0)) &&
- (NULL != scale_y) &&
- (FLT_NEQ(*scale_y, 0.0))
+ ((NULL != scale_x) && (FLT_NEQ(*scale_x, 0.0))) &&
+ ((NULL != scale_y) && (FLT_NEQ(*scale_y, 0.0)))
) {
pix_x = fabs(*scale_x);
pix_y = fabs(*scale_y);
/* adjust width and height to account new upper left */
if (ul_user) {
/* use suggested lower right corner */
- max_x = dst_gt[0] + dst_gt[1] * width;
- min_y = dst_gt[3] + dst_gt[5] * height;
+ max_x = dst_gt[0] + dst_gt[1] * _width;
+ min_y = dst_gt[3] + dst_gt[5] * _height;
+
+ /* user defined width */
+ if ((NULL != width) && (*width > 0.))
+ grid_pix_x = fabs((max_x - min_x) / ((double) _width));
+ else
+ _width = (int) ((max_x - min_x + (grid_pix_x / 2.)) / grid_pix_x);
+
+ /* user defined height */
+ if ((NULL != height) && (*height > 0.))
+ grid_pix_y = fabs((max_y - min_y) / ((double) _height));
+ else
+ _height = (int) ((max_y - min_y + (grid_pix_y / 2.)) / grid_pix_y);
- width = (int) ((max_x - min_x + (grid_pix_x / 2.)) / grid_pix_x);
- height = (int) ((max_y - min_y + (grid_pix_y / 2.)) / grid_pix_y);
dst_gt[1] = grid_pix_x;
dst_gt[5] = -1 * grid_pix_y;
- RASTER_DEBUGF(3, "new dimensions: %d x %d", width, height);
+ RASTER_DEBUGF(3, "new dimensions: %d x %d", _width, _height);
}
RASTER_DEBUGF(3, "shift is: %f, %f", grid_shift_xw, grid_shift_yw);
}
/* lower-right corner */
- max_x = min_x + dst_gt[1] * width;
- min_y = max_y + dst_gt[5] * height;
+ max_x = min_x + dst_gt[1] * _width;
+ min_y = max_y + dst_gt[5] * _height;
- width = (int) ((max_x - min_x + (pix_x / 2.)) / pix_x);
- height = (int) ((max_y - min_y + (pix_y / 2.)) / pix_y);
+ _width = (int) ((max_x - min_x + (pix_x / 2.)) / pix_x);
+ _height = (int) ((max_y - min_y + (pix_y / 2.)) / pix_y);
dst_gt[0] = min_x;
dst_gt[3] = max_y;
dst_gt[1] = pix_x;
dst_gt[5] = -1 * pix_y;
+
+ RASTER_DEBUGF(3, "new dimensions: %d x %d", _width, _height);
}
/* user-defined upper-left corner */
else if (ul_user) {
RASTER_DEBUGF(3, "Applied geotransform: %f, %f, %f, %f, %f, %f",
dst_gt[0], dst_gt[1], dst_gt[2], dst_gt[3], dst_gt[4], dst_gt[5]);
RASTER_DEBUGF(3, "Raster dimensions (width x height): %d x %d",
- width, height);
+ _width, _height);
- if (FLT_EQ(width, 0.0) || FLT_EQ(height, 0.0)) {
- rterror("rt_raster_gdal_warp: The width (%f) or height (%f) of the warped raster is zero\n", width, height);
+ if (FLT_EQ(_width, 0.0) || FLT_EQ(_height, 0.0)) {
+ rterror("rt_raster_gdal_warp: The width (%f) or height (%f) of the warped raster is zero\n", _width, _height);
GDALClose(src_ds);
}
/* create dst dataset */
- dst_ds = GDALCreate(dst_drv, "", width, height, 0, GDT_Byte, dst_options);
+ dst_ds = GDALCreate(dst_drv, "", _width, _height, 0, GDT_Byte, dst_options);
if (NULL == dst_ds) {
rterror("rt_raster_gdal_warp: Unable to create GDAL VRT dataset\n");
*
* @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
- * @param scale_x : the pixel width of the warped raster
- * @param scale_y : the pixel height of the warped raster
- * @param ul_xw : the X value of upper-left corner of the warped raster
- * @param ul_yw : the Y value of upper-left corner of the warped raster
- * @param grid_xw : the X value of point on a grid to align warped raster to
- * @param grid_yw : the Y value of point on a grid to align warped raster to
- * @param skew_x : the X skew of the warped raster
- * @param skew_y : the Y skew of the warped raster
+ * @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)
rt_raster rt_raster_gdal_warp(rt_raster raster, const char *src_srs,
const char *dst_srs,
double *scale_x, double *scale_y,
+ int *width, int *height,
double *ul_xw, double *ul_yw,
double *grid_xw, double *grid_yw,
double *skew_x, double *skew_y,
/* width */
if (!PG_ARGISNULL(3)) {
dim[0] = PG_GETARG_INT32(3);
- if (FLT_NEQ(dim[0], 0)) dim_x = &dim[0];
+ if (dim[0] < 0) dim[0] = 0;
+ if (dim[0] != 0) dim_x = &dim[0];
}
/* height */
if (!PG_ARGISNULL(4)) {
dim[1] = PG_GETARG_INT32(4);
- if (FLT_NEQ(dim[1], 0)) dim_y = &dim[1];
+ if (dim[1] < 0) dim[1] = 0;
+ if (dim[1] != 0) dim_y = &dim[1];
}
POSTGIS_RT_DEBUGF(3, "RASTER_asRaster: dim (x, y) = %d, %d", dim[0], dim[1]);
double *skew_x = NULL;
double *skew_y = NULL;
+ int dim[2] = {0};
+ int *dim_x = NULL;
+ int *dim_y = NULL;
+
POSTGIS_RT_DEBUG(3, "RASTER_resample: Starting");
/* pgraster is null, return null */
if (FLT_NEQ(skew[1], 0)) skew_y = &skew[1];
}
+ /* width */
+ if (!PG_ARGISNULL(10)) {
+ dim[0] = PG_GETARG_INT32(10);
+ if (dim[0] < 0) dim[0] = 0;
+ if (dim[0] > 0) dim_x = &dim[0];
+ }
+
+ /* height */
+ if (!PG_ARGISNULL(11)) {
+ dim[1] = PG_GETARG_INT32(11);
+ if (dim[1] < 0) dim[1] = 0;
+ if (dim[1] > 0) dim_y = &dim[1];
+ }
+
/* check that at least something is to be done */
if (
(clamp_srid(dst_srid) == SRID_UNKNOWN) &&
(grid_xw == NULL) &&
(grid_yw == NULL) &&
(skew_x == NULL) &&
- (skew_y == NULL)
+ (skew_y == NULL) &&
+ (dim_x == NULL) &&
+ (dim_y == NULL)
) {
elog(NOTICE, "No resampling parameters provided. Returning original raster");
rt_raster_destroy(raster);
rt_raster_destroy(raster);
PG_RETURN_POINTER(pgraster);
}
+ /* scale and width/height provided */
+ else if (
+ (scale_x != NULL || scale_y != NULL) &&
+ (dim_x != NULL || dim_y != NULL)
+ ) {
+ elog(NOTICE, "Scale X/Y and width/height are mutually exclusive. Only provide one. Returning original raster");
+ rt_raster_destroy(raster);
+ PG_RETURN_POINTER(pgraster);
+ }
/* get srses from srids */
/* source srs */
rast = rt_raster_gdal_warp(raster, src_srs,
dst_srs,
scale_x, scale_y,
+ dim_x, dim_y,
NULL, NULL,
grid_xw, grid_yw,
skew_x, skew_y,
srid integer DEFAULT NULL,
scalex double precision DEFAULT 0, scaley double precision DEFAULT 0,
gridx double precision DEFAULT NULL, gridy double precision DEFAULT NULL,
- skewx double precision DEFAULT 0, skewy double precision DEFAULT 0
+ skewx double precision DEFAULT 0, skewy double precision DEFAULT 0,
+ width integer DEFAULT NULL, height integer DEFAULT NULL
)
RETURNS raster
AS 'MODULE_PATHNAME', 'RASTER_resample'
AS $$ SELECT _st_resample($1, $9, $10, $2, $3, $4, $5, $6, $7, $8) $$
LANGUAGE 'sql' STABLE;
+CREATE OR REPLACE FUNCTION st_resample(
+ rast raster,
+ width integer, height integer,
+ srid integer DEFAULT NULL,
+ gridx double precision DEFAULT NULL, gridy double precision DEFAULT NULL,
+ skewx double precision DEFAULT 0, skewy double precision DEFAULT 0,
+ algorithm text DEFAULT 'NearestNeighbour', maxerr double precision DEFAULT 0.125
+)
+ RETURNS raster
+ AS $$ SELECT _st_resample($1, $9, $10, $4, NULL, NULL, $5, $6, $7, $8, $2, $3) $$
+ LANGUAGE 'sql' STABLE;
+
CREATE OR REPLACE FUNCTION st_resample(
rast raster,
ref raster,
algorithm text DEFAULT 'NearestNeighbour',
- maxerr double precision DEFAULT 0.125
+ maxerr double precision DEFAULT 0.125,
+ usescale boolean DEFAULT TRUE
)
RETURNS raster
AS $$
DECLARE
sr_id int;
+ dim_x int;
+ dim_y int;
scale_x double precision;
scale_y double precision;
grid_x double precision;
skew_x double precision;
skew_y double precision;
BEGIN
- SELECT srid, scalex, scaley, upperleftx, upperlefty, skewx, skewy INTO sr_id, scale_x, scale_y, grid_x, grid_y, skew_x, skew_y FROM st_metadata($2);
- RETURN _st_resample($1, $3, $4, sr_id, scale_x, scale_y, grid_x, grid_y, skew_x, skew_y);
+ SELECT srid, width, height, scalex, scaley, upperleftx, upperlefty, skewx, skewy INTO sr_id, dim_x, dim_y, scale_x, scale_y, grid_x, grid_y, skew_x, skew_y FROM st_metadata($2);
+
+ IF usescale IS TRUE THEN
+ dim_x := NULL;
+ dim_y := NULL;
+ ELSE
+ scale_x := NULL;
+ scale_y := NULL;
+ END IF;
+
+ RETURN _st_resample($1, $3, $4, sr_id, scale_x, scale_y, grid_x, grid_y, skew_x, skew_y, dim_x, dim_y);
END;
$$ LANGUAGE 'plpgsql' STABLE STRICT;
+CREATE OR REPLACE FUNCTION st_resample(
+ rast raster,
+ ref raster,
+ usescale boolean,
+ algorithm text DEFAULT 'NearestNeighbour',
+ maxerr double precision DEFAULT 0.125
+)
+ RETURNS raster
+ AS $$ SELECT st_resample($1, $2, $4, $5, $3) $$
+ LANGUAGE 'sql' STABLE STRICT;
+
-----------------------------------------------------------------------
-- ST_Transform
-----------------------------------------------------------------------
-- arg names changed
DROP FUNCTION IF EXISTS _ST_Resample(raster, text, double precision, integer, double precision, double precision, double precision, double precision, double precision, double precision);
+-- signature changed
+DROP FUNCTION IF EXISTS ST_Resample(raster, raster, text, double precision);
+
-- default parameters added
DROP FUNCTION IF EXISTS ST_HasNoBand(raster);
NULL, NULL,
NULL, NULL,
NULL, NULL,
+ NULL, NULL,
GRA_NearestNeighbour, -1
);
CHECK(rast);
1.5, (SELECT ST_Resample(
rast,
993310,
- 500, 500,
+ 500., 500.,
NULL, NULL,
0, 0,
'NearestNeighbor', 0.125
1.6, (SELECT ST_Resample(
rast,
NULL,
- 100, NULL
+ 100., NULL
) FROM raster_resample_src)
), (
1.7, (SELECT ST_Resample(
rast,
NULL,
- NULL::double precision, 100
+ NULL::double precision, 100.
) FROM raster_resample_src)
), (
1.8, (SELECT ST_Resample(
rast,
NULL,
- 500, 500
+ 500., 500.
) FROM raster_resample_src)
), (
1.9, (SELECT ST_Resample(
rast,
NULL,
- 250, 250,
+ 250., 250.,
NULL, NULL,
NULL, NULL
) FROM raster_resample_src)
1.10, (SELECT ST_Resample(
rast,
NULL,
- 250, 250,
+ 250., 250.,
NULL, NULL,
NULL, NULL,
'Bilinear', 0
1.15, (SELECT ST_Resample(
rast,
NULL,
- 50, 50,
+ 50., 50.,
-290, 7
) FROM raster_resample_src)
), (
1.16, (SELECT ST_Resample(
rast,
NULL,
- 121, 121,
+ 121., 121.,
0, 0
) FROM raster_resample_src)
), (
1.17, (SELECT ST_Resample(
rast,
993310,
- 50, 50,
+ 50., 50.,
-290, 7
) FROM raster_resample_src)
), (
1.18, (SELECT ST_Resample(
rast,
993309,
- 50, 50,
+ 50., 50.,
-290, 7
) FROM raster_resample_src)
), (
1.22, (SELECT ST_Resample(
rast,
993310,
- 500, 500,
+ 500., 500.,
NULL, NULL,
3, 3,
'Cubic', 0
1.23, (SELECT ST_Resample(
rast,
993310,
- 500, 500,
+ 500., 500.,
-12048, 14682,
0, 6,
'CubicSpline'
rast,
ST_MakeEmptyRaster(5, 5, -654321, 123456, 50, -100, 3, 0, 992163)
) FROM raster_resample_src)
+), (
+ 1.25, (SELECT ST_Resample(
+ rast,
+ ST_MakeEmptyRaster(5, 5, -654321, 123456, 50, -100, 3, 0, 992163),
+ TRUE
+ ) FROM raster_resample_src)
+), (
+ 1.26, (SELECT ST_Resample(
+ rast,
+ 150, 150
+ ) FROM raster_resample_src)
+), (
+ 1.27, (SELECT ST_Resample(
+ rast,
+ 150, 150,
+ 993310
+ ) FROM raster_resample_src)
+), (
+ 1.28, (SELECT ST_Resample(
+ rast,
+ ST_MakeEmptyRaster(5, 5, -654321, 123456, 100, 100, 0, 0, 992163),
+ FALSE
+ ) FROM raster_resample_src)
);
-- ST_Transform
1.22|993310|24|24|1|500.000|-500.000|3.000|3.000|950732.188|1409281.783|t|t|t
1.23|993310|25|25|1|500.000|-500.000|0.000|6.000|950452.000|1409682.000|t|t|t
1.24|992163|200|101|1|50.000|-100.000|3.000|0.000|-500021.000|600056.000|t|t|t
+1.25|992163|200|101|1|50.000|-100.000|3.000|0.000|-500021.000|600056.000|t|t|t
+1.26|992163|150|150|1|66.667|-66.667|0.000|0.000|-500000.000|600000.000|t|t|t
+1.27|993310|150|150|1|80.792|-80.792|0.000|0.000|950732.188|1409281.783|t|t|t
+1.28|992163|5|5|1|2064.200|-2291.200|0.000|0.000|-500321.000|601456.000|t|t|t
1.3|993309|12|12|1|1009.916|-1009.916|0.000|0.000|950762.305|1409088.896|t|t|t
1.4|994269|12|8|1|0.012|-0.012|0.000|0.000|-107.029|50.206|t|t|t
1.5|993310|24|24|1|500.000|-500.000|0.000|0.000|950732.188|1409281.783|t|t|t
projects / postgis / commitdiff