Addition of ST_AsRaster function to provide the ability to convert geometries into rasters.

Associated ticket is #1141.


git-svn-id: http://svn.osgeo.org/postgis/trunk@7675 b70326c6-7e19-0410-871a-916f4a2858ee
This commit is contained in:
Bborie Park 2011-07-25 15:52:36 +00:00
parent 813f391ee3
commit ed8d71db19
9 changed files with 1941 additions and 15 deletions

View file

@ -6629,3 +6629,550 @@ rt_raster rt_raster_gdal_warp(
return rast;
}
/**
* Return a raster of the provided geometry
*
* @param wkb : WKB representation of the geometry to convert
* @param wkb_len : length of the WKB representation of the geometry
* @param srs : the geometry's coordinate system in OGC WKT
* @param num_bands: number of bands in the output raster
* @param pixtype: data type of each band
* @param init: array of values to initialize each band with
* @param value: array of values for pixels of geometry
* @param nodata: array of nodata values for each band
* @param hasnodata: array flagging the presence of nodata for each band
* @param width : the number of columns of the raster
* @param height : the number of rows of the raster
* @param scale_x : the pixel width of the raster
* @param scale_y : the pixel height of the raster
* @param ul_xw : the X value of upper-left corner of the raster
* @param ul_yw : the Y value of upper-left corner of the raster
* @param grid_xw : the X value of point on grid to align raster to
* @param grid_yw : the Y value of point on grid to align raster to
* @param skew_x : the X skew of the raster
* @param skew_y : the Y skew of the raster
*
* @return the raster of the provided geometry
*/
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
) {
rt_raster rast;
int i = 0;
int noband = 0;
int banderr = 0;
int *band_list = NULL;
rt_pixtype *_pixtype = NULL;
double *_init = NULL;
double *_nodata = NULL;
uint8_t *_hasnodata = NULL;
double *_value = NULL;
double _scale_x = 0;
double _scale_y = 0;
int _width = 0;
int _height = 0;
double _skew_x = 0;
double _skew_y = 0;
OGRErr ogrerr;
OGRSpatialReferenceH src_sr = NULL;
OGRGeometryH src_geom;
OGREnvelope src_env;
int ul_user = 0;
double djunk = 0;
double grid_shift_xw = 0;
double grid_shift_yw = 0;
double grid_min_x = 0;
double grid_max_y = 0;
CPLErr cplerr;
double dst_gt[6] = {0};
GDALDriverH dst_drv = NULL;
GDALDatasetH dst_ds = NULL;
GDALRasterBandH dst_band = NULL;
RASTER_DEBUG(3, "starting");
assert(NULL != wkb);
assert(0 != wkb_len);
/* internal aliases */
_pixtype = pixtype;
_init = init;
_nodata = nodata;
_hasnodata = hasnodata;
_value = value;
/* no bands, raster is a mask */
if (num_bands < 1) {
num_bands = 1;
noband = 1;
_pixtype = (rt_pixtype *) rtalloc(sizeof(rt_pixtype));
*_pixtype = PT_8BUI;
_init = (double *) rtalloc(sizeof(double));
*_init = 0;
_nodata = (double *) rtalloc(sizeof(double));
*_nodata = 0;
_hasnodata = (uint8_t *) rtalloc(sizeof(uint8_t));
*_hasnodata = 1;
_value = (double *) rtalloc(sizeof(double));
*_value = 1;
}
assert(NULL != _pixtype);
assert(NULL != _init);
assert(NULL != _nodata);
assert(NULL != _hasnodata);
/* set OGR spatial reference */
if (NULL != srs && strlen(srs)) {
src_sr = OSRNewSpatialReference(srs);
if (NULL == src_sr) {
rterror("rt_raster_gdal_rasterize: Unable to create OSR spatial reference");
if (noband) {
rtdealloc(_pixtype);
rtdealloc(_init);
rtdealloc(_nodata);
rtdealloc(_hasnodata);
rtdealloc(_value);
}
return NULL;
}
}
/* convert WKB to OGR Geometry */
ogrerr = OGR_G_CreateFromWkb((unsigned char *) wkb, src_sr, &src_geom, wkb_len);
if (ogrerr != OGRERR_NONE) {
rterror("rt_raster_gdal_rasterize: Unable to create OGR Geometry from WKB");
if (noband) {
rtdealloc(_pixtype);
rtdealloc(_init);
rtdealloc(_nodata);
rtdealloc(_hasnodata);
rtdealloc(_value);
}
OSRDestroySpatialReference(src_sr);
OGRCleanupAll();
return NULL;
}
/* get envelope of geometry */
OGR_G_GetEnvelope(src_geom, &src_env);
RASTER_DEBUGF(3, "Suggested extent: %f, %f, %f, %f",
src_env.MinX, src_env.MaxY, src_env.MaxX, src_env.MinY);
/* user-defined scale */
if (
(NULL != scale_x) &&
(NULL != scale_y) &&
(FLT_NEQ(*scale_x, 0.0)) &&
(FLT_NEQ(*scale_y, 0.0))
) {
_scale_x = fabs(*scale_x);
_scale_y = fabs(*scale_y);
}
else if (
(NULL != width) &&
(NULL != height) &&
(FLT_NEQ(*width, 0.0)) &&
(FLT_NEQ(*height, 0.0))
) {
_width = fabs(*width);
_height = fabs(*height);
_scale_x = (src_env.MaxX - src_env.MinX) / _width;
_scale_y = (src_env.MaxY - src_env.MinY) / _height;
}
else {
rterror("rt_raster_gdal_rasterize: Values must be provided for width and height or X and Y of scale");
if (noband) {
rtdealloc(_pixtype);
rtdealloc(_init);
rtdealloc(_nodata);
rtdealloc(_hasnodata);
rtdealloc(_value);
}
OGR_G_DestroyGeometry(src_geom);
OSRDestroySpatialReference(src_sr);
OGRCleanupAll();
return NULL;
}
RASTER_DEBUGF(3, "scale (x, y) = %f, %f", _scale_x, _scale_y);
RASTER_DEBUGF(3, "dim (x, y) = %d, %d", _width, _height);
/* user-defined skew */
if (NULL != skew_x)
_skew_x = *skew_x;
if (NULL != skew_y)
_skew_y = *skew_y;
/* user-specified upper-left corner */
if (
NULL != ul_xw &&
NULL != ul_yw
) {
src_env.MinX = *ul_xw;
src_env.MaxY = *ul_yw;
ul_user = 1;
}
else if (
((NULL != ul_xw) && (NULL == ul_yw)) ||
((NULL == ul_xw) && (NULL != ul_yw))
) {
rterror("rt_raster_gdal_rasterize: Both X and Y upper-left corner values must be provided");
if (noband) {
rtdealloc(_pixtype);
rtdealloc(_init);
rtdealloc(_nodata);
rtdealloc(_hasnodata);
rtdealloc(_value);
}
OGR_G_DestroyGeometry(src_geom);
OSRDestroySpatialReference(src_sr);
OGRCleanupAll();
return NULL;
}
/* alignment only considered if upper-left corner not provided */
if (
!ul_user && (
(NULL != grid_xw) || (NULL != grid_yw)
)
) {
if (
((NULL != grid_xw) && (NULL == grid_yw)) ||
((NULL == grid_xw) && (NULL != grid_yw))
) {
rterror("rt_raster_gdal_rasterize: Both X and Y alignment values must be provided\n");
if (noband) {
rtdealloc(_pixtype);
rtdealloc(_init);
rtdealloc(_nodata);
rtdealloc(_hasnodata);
rtdealloc(_value);
}
OGR_G_DestroyGeometry(src_geom);
OSRDestroySpatialReference(src_sr);
OGRCleanupAll();
return NULL;
}
/* grid shift of upper left to match alignment grid */
grid_shift_xw = _scale_x * modf(fabs(*grid_xw - src_env.MinX) / _scale_x, &djunk);
grid_shift_yw = _scale_y * modf(fabs(*grid_yw - src_env.MaxY) / _scale_y, &djunk);
/* shift along X axis for upper left */
if (FLT_NEQ(grid_shift_xw, 0.) && FLT_NEQ(grid_shift_xw, _scale_x)) {
grid_min_x = src_env.MinX - grid_shift_xw;
grid_min_x = modf(fabs(*grid_xw - grid_min_x) / _scale_x, &djunk);
if (FLT_NEQ(grid_min_x, 0.) && FLT_NEQ(grid_min_x, 1.))
grid_shift_xw = _scale_x - grid_shift_xw;
grid_min_x = src_env.MinX - grid_shift_xw;
src_env.MinX = grid_min_x;
}
/* shift along Y axis for upper left */
if (FLT_NEQ(grid_shift_yw, 0.) && FLT_NEQ(grid_shift_yw, _scale_y)) {
grid_max_y = src_env.MaxY + grid_shift_yw;
grid_max_y = modf(fabs(*grid_yw - grid_max_y) / _scale_y, &djunk);
if (FLT_NEQ(grid_max_y, 0.))
grid_shift_yw = _scale_y - fabs(grid_shift_yw);
grid_max_y = src_env.MaxY + grid_shift_yw;
src_env.MaxY = grid_max_y;
}
RASTER_DEBUGF(3, "shift is: %f, %f", grid_shift_xw, grid_shift_yw);
RASTER_DEBUGF(3, "new ul is: %f, %f", grid_min_x, grid_max_y);
}
/* raster dimensions */
if (!_width && !_height) {
_width = (int) ((src_env.MaxX - src_env.MinX + (_scale_x / 2.)) / _scale_x);
_height = (int) ((src_env.MaxY - src_env.MinY + (_scale_y / 2.)) / _scale_y);
}
/* min and max are same, a point? */
if (
FLT_EQ(src_env.MaxX, src_env.MinX) &&
FLT_EQ(src_env.MaxY, src_env.MinY)
) {
RASTER_DEBUGF(3, "MinX = MaxX, MinY = MaxY. Explicitly setting width and height to 1",
_width, _height);
_width = 1;
_height = 1;
/* set the point to the center of the pixel */
src_env.MinX -= (_scale_x / 2.);
src_env.MaxX += (_scale_x / 2.);
src_env.MinY -= (_scale_y / 2.);
src_env.MaxY += (_scale_y / 2.);
}
/* specify geotransform */
dst_gt[0] = src_env.MinX;
dst_gt[1] = _scale_x;
dst_gt[2] = _skew_x;
dst_gt[3] = src_env.MaxY;
dst_gt[4] = _skew_y;
dst_gt[5] = -1 * _scale_y;
RASTER_DEBUGF(3, "Applied extent: %f, %f, %f, %f",
src_env.MinX, src_env.MaxY, src_env.MaxX, src_env.MinY);
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);
/* load GDAL mem */
GDALRegister_MEM();
dst_drv = GDALGetDriverByName("MEM");
if (NULL == dst_drv) {
rterror("rt_raster_gdal_rasterize: Unable to load the MEM GDAL driver");
if (noband) {
rtdealloc(_pixtype);
rtdealloc(_init);
rtdealloc(_nodata);
rtdealloc(_hasnodata);
rtdealloc(_value);
}
OGR_G_DestroyGeometry(src_geom);
OSRDestroySpatialReference(src_sr);
OGRCleanupAll();
return NULL;
}
dst_ds = GDALCreate(dst_drv, "", _width, _height, 0, GDT_Byte, NULL);
if (NULL == dst_ds) {
rterror("rt_raster_gdal_rasterize: Could not create a GDALDataset to rasterize the geometry into");
if (noband) {
rtdealloc(_pixtype);
rtdealloc(_init);
rtdealloc(_nodata);
rtdealloc(_hasnodata);
rtdealloc(_value);
}
OGR_G_DestroyGeometry(src_geom);
OSRDestroySpatialReference(src_sr);
OGRCleanupAll();
GDALDeregisterDriver(dst_drv);
GDALDestroyDriver(dst_drv);
return NULL;
}
/* set geotransform */
cplerr = GDALSetGeoTransform(dst_ds, dst_gt);
if (cplerr != CE_None) {
rterror("rt_raster_gdal_rasterize: Could not set geotransform on GDALDataset");
if (noband) {
rtdealloc(_pixtype);
rtdealloc(_init);
rtdealloc(_nodata);
rtdealloc(_hasnodata);
rtdealloc(_value);
}
OGR_G_DestroyGeometry(src_geom);
OSRDestroySpatialReference(src_sr);
OGRCleanupAll();
GDALClose(dst_ds);
GDALDeregisterDriver(dst_drv);
GDALDestroyDriver(dst_drv);
return NULL;
}
/* set SRS */
if (NULL != src_sr) {
cplerr = GDALSetProjection(dst_ds, srs);
if (cplerr != CE_None) {
rterror("rt_raster_gdal_rasterize: Could not set projection on GDALDataset");
if (noband) {
rtdealloc(_pixtype);
rtdealloc(_init);
rtdealloc(_nodata);
rtdealloc(_hasnodata);
rtdealloc(_value);
}
OGR_G_DestroyGeometry(src_geom);
OSRDestroySpatialReference(src_sr);
OGRCleanupAll();
GDALClose(dst_ds);
GDALDeregisterDriver(dst_drv);
GDALDestroyDriver(dst_drv);
return NULL;
}
}
/* create and set bands */
for (i = 0; i < num_bands; i++) {
banderr = 0;
do {
/* add band */
cplerr = GDALAddBand(dst_ds, rt_util_pixtype_to_gdal_datatype(_pixtype[i]), NULL);
if (cplerr != CE_None) {
rterror("rt_raster_gdal_rasterize: Unable to add band to GDALDataset");
banderr = 1;
break;
}
dst_band = GDALGetRasterBand(dst_ds, i + 1);
if (NULL == dst_band) {
rterror("rt_raster_gdal_rasterize: Unable to get band %d from GDALDataset", i + 1);
banderr = 1;
break;
}
/* nodata value */
if (_hasnodata[i]) {
cplerr = GDALSetRasterNoDataValue(dst_band, _nodata[i]);
if (cplerr != CE_None) {
rterror("rt_raster_gdal_rasterize: Unable to set nodata value");
banderr = 1;
break;
}
}
/* initial value */
cplerr = GDALFillRaster(dst_band, _init[i], 0);
if (cplerr != CE_None) {
rterror("rt_raster_gdal_rasterize: Unable to set initial value");
banderr = 1;
break;
}
}
while (0);
if (banderr) {
if (noband) {
rtdealloc(_pixtype);
rtdealloc(_init);
rtdealloc(_nodata);
rtdealloc(_hasnodata);
rtdealloc(_value);
}
OGR_G_DestroyGeometry(src_geom);
OSRDestroySpatialReference(src_sr);
OGRCleanupAll();
GDALClose(dst_ds);
GDALDeregisterDriver(dst_drv);
GDALDestroyDriver(dst_drv);
return NULL;
}
}
band_list = (int *) rtalloc(sizeof(double) * num_bands);
for (i = 0; i < num_bands; i++) band_list[i] = i + 1;
/* burn geometry */
cplerr = GDALRasterizeGeometries(
dst_ds,
num_bands, band_list,
1, &src_geom,
NULL, NULL,
_value,
NULL,
NULL, NULL
);
rtdealloc(band_list);
if (cplerr != CE_None) {
rterror("rt_raster_gdal_rasterize: Unable to rasterize geometry");
if (noband) {
rtdealloc(_pixtype);
rtdealloc(_init);
rtdealloc(_nodata);
rtdealloc(_hasnodata);
rtdealloc(_value);
}
OGR_G_DestroyGeometry(src_geom);
OSRDestroySpatialReference(src_sr);
OGRCleanupAll();
GDALClose(dst_ds);
GDALDeregisterDriver(dst_drv);
GDALDestroyDriver(dst_drv);
return NULL;
}
/* convert gdal dataset to raster */
RASTER_DEBUG(3, "Converting GDAL dataset to raster");
rast = rt_raster_from_gdal_dataset(dst_ds);
if (noband) {
rtdealloc(_pixtype);
rtdealloc(_init);
rtdealloc(_nodata);
rtdealloc(_hasnodata);
rtdealloc(_value);
}
OGR_G_DestroyGeometry(src_geom);
OSRDestroySpatialReference(src_sr);
OGRCleanupAll();
GDALClose(dst_ds);
GDALDeregisterDriver(dst_drv);
GDALDestroyDriver(dst_drv);
if (NULL == rast) {
rterror("rt_raster_gdal_rasterize: Unable to rasterize geometry\n");
return NULL;
}
RASTER_DEBUG(3, "done");
return rast;
}

View file

@ -954,6 +954,42 @@ rt_raster rt_raster_gdal_warp(rt_raster raster, const char *src_srs,
double *skew_x, double *skew_y,
GDALResampleAlg resample_alg, double max_err);
/**
* Return a raster of the provided geometry
*
* @param wkb : WKB representation of the geometry to convert
* @param wkb_len : length of the WKB representation of the geometry
* @param srs : the geometry's coordinate system in OGC WKT
* @param num_bands: number of bands in the output raster
* @param pixtype: data type of each band
* @param init: array of values to initialize each band with
* @param value: array of values for pixels of geometry
* @param nodata: array of nodata values for each band
* @param hasnodata: array flagging the presence of nodata for each band
* @param width : the number of columns of the raster
* @param height : the number of rows of the raster
* @param scale_x : the pixel width of the raster
* @param scale_y : the pixel height of the raster
* @param ul_xw : the X value of upper-left corner of the raster
* @param ul_yw : the Y value of upper-left corner of the raster
* @param grid_xw : the X value of point on grid to align raster to
* @param grid_yw : the Y value of point on grid to align raster to
* @param skew_x : the X skew of the raster
* @param skew_y : the Y skew of the raster
*
* @return the raster of the provided geometry
*/
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);
/*- utilities -------------------------------------------------------*/
/*

View file

@ -238,6 +238,9 @@ Datum RASTER_reclass(PG_FUNCTION_ARGS);
Datum RASTER_asGDALRaster(PG_FUNCTION_ARGS);
Datum RASTER_getGDALDrivers(PG_FUNCTION_ARGS);
/* rasterize a geometry */
Datum RASTER_asRaster(PG_FUNCTION_ARGS);
/* resample a raster */
Datum RASTER_resample(PG_FUNCTION_ARGS);
@ -4676,6 +4679,585 @@ Datum RASTER_getGDALDrivers(PG_FUNCTION_ARGS)
}
}
/**
* Rasterize a geometry
*/
PG_FUNCTION_INFO_V1(RASTER_asRaster);
Datum RASTER_asRaster(PG_FUNCTION_ARGS)
{
#ifdef GSERIALIZED_ON
GSERIALIZED *pggeom = NULL;
#else
unsigned char *pggeom = NULL;
#endif
LWGEOM *geom = NULL;
rt_raster rast = NULL;
rt_pgraster *pgrast = NULL;
unsigned char *wkb;
size_t wkb_len = 0;
unsigned char variant = WKB_SFSQL;
double scale[2] = {0};
double *scale_x = NULL;
double *scale_y = NULL;
int dim[2] = {0};
int *dim_x = NULL;
int *dim_y = NULL;
ArrayType *array;
Oid etype;
Datum *e;
bool *nulls;
int16 typlen;
bool typbyval;
char typalign;
int ndims = 1;
int *dims;
int *lbs;
int n = 0;
int i = 0;
int j = 0;
int haserr = 0;
text *pixeltypetext = NULL;
char *pixeltype = NULL;
rt_pixtype pixtype = PT_END;
rt_pixtype *pixtypes = NULL;
int pixtypes_len = 0;
double *values = NULL;
int values_len = 0;
uint8_t *hasnodatas = NULL;
double *nodatavals = NULL;
int nodatavals_len = 0;
double ulw[2] = {0};
double *ul_xw = NULL;
double *ul_yw = NULL;
double gridw[2] = {0};
double *grid_xw = NULL;
double *grid_yw = NULL;
double skew[2] = {0};
double *skew_x = NULL;
double *skew_y = NULL;
uint32_t num_bands = 0;
int srid = SRID_UNKNOWN;
char *srs = NULL;
POSTGIS_RT_DEBUG(3, "RASTER_asRaster: Starting");
/* based upon LWGEOM_asBinary function in postgis/lwgeom_ogc.c */
/* Get the geometry */
if (PG_ARGISNULL(0)) PG_RETURN_NULL();
#ifdef GSERIALIZED_ON
pggeom = (GSERIALIZED *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
geom = lwgeom_from_gserialized(pggeom);
#else
pggeom = (unsigned char *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
geom = lwgeom_deserialize(SERIALIZED_FORM(pggeom));
#endif
/* Get a 2D version of the geometry if necessary */
if (FLAGS_NDIMS(geom->flags) > 2) {
LWGEOM *geom2d = lwgeom_force_2d(geom);
lwgeom_free(geom);
geom = geom2d;
}
/* scale x */
if (!PG_ARGISNULL(1)) {
scale[0] = PG_GETARG_FLOAT8(1);
if (FLT_NEQ(scale[0], 0)) scale_x = &scale[0];
}
/* scale y */
if (!PG_ARGISNULL(2)) {
scale[1] = PG_GETARG_FLOAT8(2);
if (FLT_NEQ(scale[1], 0)) scale_y = &scale[1];
}
POSTGIS_RT_DEBUGF(3, "RASTER_asRaster: scale (x, y) = %f, %f", scale[0], scale[1]);
/* width */
if (!PG_ARGISNULL(3)) {
dim[0] = PG_GETARG_INT32(3);
if (FLT_NEQ(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];
}
POSTGIS_RT_DEBUGF(3, "RASTER_asRaster: dim (x, y) = %d, %d", dim[0], dim[1]);
/* pixeltype */
if (!PG_ARGISNULL(5)) {
array = PG_GETARG_ARRAYTYPE_P(5);
etype = ARR_ELEMTYPE(array);
get_typlenbyvalalign(etype, &typlen, &typbyval, &typalign);
switch (etype) {
case TEXTOID:
break;
default:
elog(ERROR, "RASTER_asRaster: Invalid data type for pixeltype");
lwgeom_free(geom);
PG_FREE_IF_COPY(pggeom, 0);
PG_RETURN_NULL();
break;
}
ndims = ARR_NDIM(array);
dims = ARR_DIMS(array);
lbs = ARR_LBOUND(array);
deconstruct_array(array, etype, typlen, typbyval, typalign, &e,
&nulls, &n);
if (n) {
pixtypes = (rt_pixtype *) palloc(sizeof(rt_pixtype) * n);
/* clean each pixeltype */
for (i = 0, j = 0; i < n; i++) {
if (nulls[i]) {
pixtypes[j++] = PT_64BF;
continue;
}
pixeltype = NULL;
switch (etype) {
case TEXTOID:
pixeltypetext = (text *) DatumGetPointer(e[i]);
if (NULL == pixeltypetext) break;
pixeltype = text_to_cstring(pixeltypetext);
/* trim string */
pixeltype = trim(pixeltype);
POSTGIS_RT_DEBUGF(3, "RASTER_asRaster: pixeltype is '%s'", pixeltype);
break;
}
if (strlen(pixeltype)) {
pixtype = rt_pixtype_index_from_name(pixeltype);
if (pixtype == PT_END) {
elog(ERROR, "RASTER_asRaster: Invalid pixel type provided: %s", pixeltype);
pfree(pixtypes);
lwgeom_free(geom);
PG_FREE_IF_COPY(pggeom, 0);
PG_RETURN_NULL();
}
pixtypes[j] = pixtype;
j++;
}
}
if (j > 0) {
/* trim allocation */
pixtypes = repalloc(pixtypes, j * sizeof(rt_pixtype));
pixtypes_len = j;
}
else {
pfree(pixtypes);
pixtypes = NULL;
pixtypes_len = 0;
}
}
}
#if POSTGIS_DEBUG_LEVEL > 0
for (i = 0; i < pixtypes_len; i++)
POSTGIS_RT_DEBUGF(3, "RASTER_asRaster: pixtypes[%d] = %d", i, (int) pixtypes[i]);
#endif
/* value */
if (!PG_ARGISNULL(6)) {
array = PG_GETARG_ARRAYTYPE_P(6);
etype = ARR_ELEMTYPE(array);
get_typlenbyvalalign(etype, &typlen, &typbyval, &typalign);
switch (etype) {
case FLOAT4OID:
case FLOAT8OID:
break;
default:
elog(ERROR, "RASTER_asRaster: Invalid data type for value");
if (pixtypes_len) pfree(pixtypes);
lwgeom_free(geom);
PG_FREE_IF_COPY(pggeom, 0);
PG_RETURN_NULL();
break;
}
ndims = ARR_NDIM(array);
dims = ARR_DIMS(array);
lbs = ARR_LBOUND(array);
deconstruct_array(array, etype, typlen, typbyval, typalign, &e,
&nulls, &n);
if (n) {
values = (double *) palloc(sizeof(double) * n);
for (i = 0, j = 0; i < n; i++) {
if (nulls[i]) {
values[j++] = 1;
continue;
}
switch (etype) {
case FLOAT4OID:
values[j] = (double) DatumGetFloat4(e[i]);
break;
case FLOAT8OID:
values[j] = (double) DatumGetFloat8(e[i]);
break;
}
POSTGIS_RT_DEBUGF(3, "RASTER_asRaster: values[%d] = %f", j, values[j]);
j++;
}
if (j > 0) {
/* trim allocation */
values = repalloc(values, j * sizeof(double));
values_len = j;
}
else {
pfree(values);
values = NULL;
values_len = 0;
}
}
}
#if POSTGIS_DEBUG_LEVEL > 0
for (i = 0; i < values_len; i++)
POSTGIS_RT_DEBUGF(3, "RASTER_asRaster: values[%d] = %f", i, values[i]);
#endif
/* nodataval */
if (!PG_ARGISNULL(7)) {
array = PG_GETARG_ARRAYTYPE_P(7);
etype = ARR_ELEMTYPE(array);
get_typlenbyvalalign(etype, &typlen, &typbyval, &typalign);
switch (etype) {
case FLOAT4OID:
case FLOAT8OID:
break;
default:
elog(ERROR, "RASTER_asRaster: Invalid data type for nodataval");
if (pixtypes_len) pfree(pixtypes);
if (values_len) pfree(values);
lwgeom_free(geom);
PG_FREE_IF_COPY(pggeom, 0);
PG_RETURN_NULL();
break;
}
ndims = ARR_NDIM(array);
dims = ARR_DIMS(array);
lbs = ARR_LBOUND(array);
deconstruct_array(array, etype, typlen, typbyval, typalign, &e,
&nulls, &n);
if (n) {
nodatavals = (double *) palloc(sizeof(double) * n);
hasnodatas = (uint8_t *) palloc(sizeof(uint8_t) * n);
for (i = 0, j = 0; i < n; i++) {
if (nulls[i]) {
hasnodatas[j] = 0;
nodatavals[j] = 0;
j++;
continue;
}
hasnodatas[j] = 1;
switch (etype) {
case FLOAT4OID:
nodatavals[j] = (double) DatumGetFloat4(e[i]);
break;
case FLOAT8OID:
nodatavals[j] = (double) DatumGetFloat8(e[i]);
break;
}
POSTGIS_RT_DEBUGF(3, "RASTER_asRaster: hasnodatas[%d] = %d", j, hasnodatas[j]);
POSTGIS_RT_DEBUGF(3, "RASTER_asRaster: nodatavals[%d] = %f", j, nodatavals[j]);
j++;
}
if (j > 0) {
/* trim allocation */
nodatavals = repalloc(nodatavals, j * sizeof(double));
hasnodatas = repalloc(hasnodatas, j * sizeof(uint8_t));
nodatavals_len = j;
}
else {
pfree(nodatavals);
pfree(hasnodatas);
nodatavals = NULL;
hasnodatas = NULL;
nodatavals_len = 0;
}
}
}
#if POSTGIS_DEBUG_LEVEL > 0
for (i = 0; i < nodatavals_len; i++) {
POSTGIS_RT_DEBUGF(3, "RASTER_asRaster: hasnodatas[%d] = %d", i, hasnodatas[i]);
POSTGIS_RT_DEBUGF(3, "RASTER_asRaster: nodatavals[%d] = %f", i, nodatavals[i]);
}
#endif
/* upperleftx */
if (!PG_ARGISNULL(8)) {
ulw[0] = PG_GETARG_FLOAT8(8);
ul_xw = &ulw[0];
}
/* upperlefty */
if (!PG_ARGISNULL(9)) {
ulw[1] = PG_GETARG_FLOAT8(9);
ul_yw = &ulw[1];
}
POSTGIS_RT_DEBUGF(3, "RASTER_asRaster: upperleft (x, y) = %f, %f", ulw[0], ulw[1]);
/* gridx */
if (!PG_ARGISNULL(10)) {
gridw[0] = PG_GETARG_FLOAT8(10);
grid_xw = &gridw[0];
}
/* gridy */
if (!PG_ARGISNULL(11)) {
gridw[1] = PG_GETARG_FLOAT8(11);
grid_yw = &gridw[1];
}
POSTGIS_RT_DEBUGF(3, "RASTER_asRaster: grid (x, y) = %f, %f", gridw[0], gridw[1]);
/* check dependent variables */
haserr = 0;
do {
/* only part of scale provided */
if (
(scale_x == NULL && scale_y != NULL) ||
(scale_x != NULL && scale_y == NULL)
) {
elog(NOTICE, "Values must be provided for both X and Y of scale if one is specified");
haserr = 1;
break;
}
/* only part of dimension provided */
if (
(dim_x == NULL && dim_y != NULL) ||
(dim_x != NULL && dim_y == NULL)
) {
elog(NOTICE, "Values must be provided for both width and height if one is specified");
haserr = 1;
break;
}
/* scale and dimension provided */
if (
(scale_x != NULL && scale_y != NULL) &&
(dim_x != NULL && dim_y != NULL)
) {
elog(NOTICE, "Values provided for X and Y of scale and width and height. Using the width and height");
scale_x = NULL;
scale_y = NULL;
break;
}
/* neither scale or dimension provided */
if (
(scale_x == NULL && scale_y == NULL) &&
(dim_x == NULL && dim_y == NULL)
) {
elog(NOTICE, "Values must be provided for X and Y of scale or width and height");
haserr = 1;
break;
}
/* only part of upper-left provided */
if (
(ul_xw == NULL && ul_yw != NULL) ||
(ul_xw != NULL && ul_yw == NULL)
) {
elog(NOTICE, "Values must be provided for both X and Y when specifying the upper-left corner");
haserr = 1;
break;
}
/* only part of alignment provided */
if (
(grid_xw == NULL && grid_yw != NULL) ||
(grid_xw != NULL && grid_yw == NULL)
) {
elog(NOTICE, "Values must be provided for both X and Y when specifying the alignment");
haserr = 1;
break;
}
/* upper-left and alignment provided */
if (
(ul_xw != NULL && ul_yw != NULL) &&
(grid_xw != NULL && grid_yw != NULL)
) {
elog(NOTICE, "Values provided for both X and Y of upper-left corner and alignment. Using the values of upper-left corner");
grid_xw = NULL;
grid_yw = NULL;
break;
}
}
while (0);
if (haserr) {
if (pixtypes_len) pfree(pixtypes);
if (values_len) pfree(values);
if (nodatavals_len) {
pfree(nodatavals);
pfree(hasnodatas);
}
lwgeom_free(geom);
PG_FREE_IF_COPY(pggeom, 0);
PG_RETURN_NULL();
}
/* skewx */
if (!PG_ARGISNULL(12)) {
skew[0] = PG_GETARG_FLOAT8(12);
if (FLT_NEQ(skew[0], 0)) skew_x = &skew[0];
}
/* skewy */
if (!PG_ARGISNULL(13)) {
skew[1] = PG_GETARG_FLOAT8(13);
if (FLT_NEQ(skew[1], 0)) skew_y = &skew[1];
}
POSTGIS_RT_DEBUGF(3, "RASTER_asRaster: skew (x, y) = %f, %f", skew[0], skew[1]);
/* get geometry's srid */
#ifdef GSERIALIZED_ON
srid = gserialized_get_srid(pggeom);
#else
srid = lwgeom_getsrid(pggeom);
#endif
POSTGIS_RT_DEBUGF(3, "RASTER_asRaster: srid = %d", srid);
if (srid != SRID_UNKNOWN) {
srs = getSRTextSPI(srid);
if (NULL == srs) {
elog(ERROR, "RASTER_asRaster: Could not find srtext for SRID (%d)", srid);
if (pixtypes_len) pfree(pixtypes);
if (values_len) pfree(values);
if (nodatavals_len) {
pfree(hasnodatas);
pfree(nodatavals);
}
lwgeom_free(geom);
PG_FREE_IF_COPY(pggeom, 0);
PG_RETURN_NULL();
}
POSTGIS_RT_DEBUGF(3, "RASTER_asRaster: srs is %s", srs);
}
else
srs = NULL;
/* determine number of bands */
/* MIN macro is from GDAL's cpl_port.h */
num_bands = MIN(pixtypes_len, values_len);
num_bands = MIN(num_bands, nodatavals_len);
POSTGIS_RT_DEBUGF(3, "RASTER_asRaster: pixtypes_len = %d", pixtypes_len);
POSTGIS_RT_DEBUGF(3, "RASTER_asRaster: values_len = %d", values_len);
POSTGIS_RT_DEBUGF(3, "RASTER_asRaster: nodatavals_len = %d", nodatavals_len);
POSTGIS_RT_DEBUGF(3, "RASTER_asRaster: num_bands = %d", num_bands);
/* warn of imbalanced number of band elements */
if (!(
(pixtypes_len == values_len) &&
(values_len == nodatavals_len)
)) {
elog(
NOTICE,
"Imbalanced number of values provided for pixeltype (%d), value (%d) and nodataval (%d). Using the first %d values of each parameter",
pixtypes_len,
values_len,
nodatavals_len,
num_bands
);
}
/* get wkb of geometry */
POSTGIS_RT_DEBUG(3, "RASTER_asRaster: getting wkb of geometry");
wkb = lwgeom_to_wkb(geom, variant, &wkb_len);
lwgeom_free(geom);
PG_FREE_IF_COPY(pggeom, 0);
/* rasterize geometry */
POSTGIS_RT_DEBUG(3, "RASTER_asRaster: rasterizing geometry");
/* use nodatavals for the init parameter */
rast = rt_raster_gdal_rasterize(wkb,
(uint32_t) wkb_len, srs,
num_bands, pixtypes,
nodatavals, values,
nodatavals, hasnodatas,
dim_x, dim_y,
scale_x, scale_y,
ul_xw, ul_yw,
grid_xw, grid_yw,
skew_x, skew_y
);
if (pixtypes_len) pfree(pixtypes);
if (values_len) pfree(values);
if (nodatavals_len) {
pfree(hasnodatas);
pfree(nodatavals);
}
if (!rast) {
elog(ERROR, "RASTER_asRaster: Could not create rasterize geometry");
PG_RETURN_NULL();
}
/* add target srid */
rt_raster_set_srid(rast, srid);
pgrast = rt_raster_serialize(rast);
rt_raster_destroy(rast);
if (NULL == pgrast) PG_RETURN_NULL();
POSTGIS_RT_DEBUG(3, "RASTER_asRaster: done");
SET_VARSIZE(pgrast, pgrast->size);
PG_RETURN_POINTER(pgrast);
}
/**
* Resample a raster
*/

View file

@ -1439,6 +1439,177 @@ CREATE OR REPLACE FUNCTION st_aspng(rast raster, nband int, compression int)
AS $$ SELECT st_aspng($1, ARRAY[$2], $3) $$
LANGUAGE 'SQL' IMMUTABLE STRICT;
-----------------------------------------------------------------------
-- ST_AsRaster
-----------------------------------------------------------------------
-- None of the ST_AsRaster can be strict as some parameters can be NULL
CREATE OR REPLACE FUNCTION _st_asraster(
geom geometry,
scalex double precision DEFAULT 0, scaley double precision DEFAULT 0,
width integer DEFAULT 0, height integer DEFAULT 0,
pixeltype text[] DEFAULT ARRAY['8BUI']::text[],
value double precision[] DEFAULT ARRAY[1]::double precision[],
nodataval double precision[] DEFAULT ARRAY[0]::double precision[],
upperleftx double precision DEFAULT NULL, upperlefty double precision DEFAULT NULL,
gridx double precision DEFAULT NULL, gridy double precision DEFAULT NULL,
skewx double precision DEFAULT 0, skewy double precision DEFAULT 0
)
RETURNS raster
AS 'MODULE_PATHNAME', 'RASTER_asRaster'
LANGUAGE 'C' STABLE;
CREATE OR REPLACE FUNCTION st_asraster(
geom geometry,
scalex double precision, scaley double precision,
gridx double precision DEFAULT NULL, gridy double precision DEFAULT NULL,
pixeltype text[] DEFAULT ARRAY['8BUI']::text[],
value double precision[] DEFAULT ARRAY[1]::double precision[],
nodataval double precision[] DEFAULT ARRAY[0]::double precision[],
skewx double precision DEFAULT 0, skewy double precision DEFAULT 0
)
RETURNS raster
AS $$ SELECT _st_asraster($1, $2, $3, NULL, NULL, $6, $7, $8, NULL, NULL, $4, $5, $9, $10) $$
LANGUAGE 'sql' STABLE;
CREATE OR REPLACE FUNCTION st_asraster(
geom geometry,
scalex double precision, scaley double precision,
pixeltype text[],
value double precision[] DEFAULT ARRAY[1]::double precision[],
nodataval double precision[] DEFAULT ARRAY[0]::double precision[],
upperleftx double precision DEFAULT NULL, upperlefty double precision DEFAULT NULL,
skewx double precision DEFAULT 0, skewy double precision DEFAULT 0
)
RETURNS raster
AS $$ SELECT _st_asraster($1, $2, $3, NULL, NULL, $4, $5, $6, $7, $8, NULL, NULL, $9, $10) $$
LANGUAGE 'sql' STABLE;
CREATE OR REPLACE FUNCTION st_asraster(
geom geometry,
width integer, height integer,
gridx double precision DEFAULT NULL, gridy double precision DEFAULT NULL,
pixeltype text[] DEFAULT ARRAY['8BUI']::text[],
value double precision[] DEFAULT ARRAY[1]::double precision[],
nodataval double precision[] DEFAULT ARRAY[0]::double precision[],
skewx double precision DEFAULT 0, skewy double precision DEFAULT 0
)
RETURNS raster
AS $$ SELECT _st_asraster($1, NULL, NULL, $2, $3, $6, $7, $8, NULL, NULL, $4, $5, $9, $10) $$
LANGUAGE 'sql' STABLE;
CREATE OR REPLACE FUNCTION st_asraster(
geom geometry,
width integer, height integer,
pixeltype text[],
value double precision[] DEFAULT ARRAY[1]::double precision[],
nodataval double precision[] DEFAULT ARRAY[0]::double precision[],
upperleftx double precision DEFAULT NULL, upperlefty double precision DEFAULT NULL,
skewx double precision DEFAULT 0, skewy double precision DEFAULT 0
)
RETURNS raster
AS $$ SELECT _st_asraster($1, NULL, NULL, $2, $3, $4, $5, $6, $7, $8, NULL, NULL, $9, $10) $$
LANGUAGE 'sql' STABLE;
CREATE OR REPLACE FUNCTION st_asraster(
geom geometry,
scalex double precision, scaley double precision,
gridx double precision, gridy double precision,
pixeltype text,
value double precision DEFAULT 1,
nodataval double precision DEFAULT 0,
skewx double precision DEFAULT 0, skewy double precision DEFAULT 0
)
RETURNS raster
AS $$ SELECT _st_asraster($1, $2, $3, NULL, NULL, ARRAY[$6]::text[], ARRAY[$7]::double precision[], ARRAY[$8]::double precision[], NULL, NULL, $4, $5, $9, $10) $$
LANGUAGE 'sql' STABLE;
CREATE OR REPLACE FUNCTION st_asraster(
geom geometry,
scalex double precision, scaley double precision,
pixeltype text,
value double precision DEFAULT 1,
nodataval double precision DEFAULT 0,
upperleftx double precision DEFAULT NULL, upperlefty double precision DEFAULT NULL,
skewx double precision DEFAULT 0, skewy double precision DEFAULT 0
)
RETURNS raster
AS $$ SELECT _st_asraster($1, $2, $3, NULL, NULL, ARRAY[$4]::text[], ARRAY[$5]::double precision[], ARRAY[$6]::double precision[], $7, $8, NULL, NULL, $9, $10) $$
LANGUAGE 'sql' STABLE;
CREATE OR REPLACE FUNCTION st_asraster(
geom geometry,
width integer, height integer,
gridx double precision, gridy double precision,
pixeltype text,
value double precision DEFAULT 1,
nodataval double precision DEFAULT 0,
skewx double precision DEFAULT 0, skewy double precision DEFAULT 0
)
RETURNS raster
AS $$ SELECT _st_asraster($1, NULL, NULL, $2, $3, ARRAY[$6]::text[], ARRAY[$7]::double precision[], ARRAY[$8]::double precision[], NULL, NULL, $4, $5, $9, $10) $$
LANGUAGE 'sql' STABLE;
CREATE OR REPLACE FUNCTION st_asraster(
geom geometry,
width integer, height integer,
pixeltype text,
value double precision DEFAULT 1,
nodataval double precision DEFAULT 0,
upperleftx double precision DEFAULT NULL, upperlefty double precision DEFAULT NULL,
skewx double precision DEFAULT 0, skewy double precision DEFAULT 0
)
RETURNS raster
AS $$ SELECT _st_asraster($1, NULL, NULL, $2, $3, ARRAY[$4]::text[], ARRAY[$5]::double precision[], ARRAY[$6]::double precision[], $7, $8, NULL, NULL,$9, $10) $$
LANGUAGE 'sql' STABLE;
CREATE OR REPLACE FUNCTION st_asraster(
geom geometry,
ref raster,
pixeltype text[] DEFAULT ARRAY['8BUI']::text[],
value double precision[] DEFAULT ARRAY[1]::double precision[],
nodataval double precision[] DEFAULT ARRAY[0]::double precision[]
)
RETURNS raster
AS $$
DECLARE
g geometry;
g_srid integer;
ul_x double precision;
ul_y double precision;
scale_x double precision;
scale_y double precision;
skew_x double precision;
skew_y double precision;
sr_id integer;
BEGIN
SELECT upperleftx, upperlefty, scalex, scaley, skewx, skewy, srid INTO ul_x, ul_y, scale_x, scale_y, skew_x, skew_y, sr_id FROM ST_Metadata(ref);
--RAISE NOTICE '%, %, %, %, %, %, %', ul_x, ul_y, scale_x, scale_y, skew_x, skew_y, sr_id;
-- geometry and raster has different SRID
g_srid := ST_SRID(geom);
IF g_srid != sr_id THEN
RAISE NOTICE 'The geometry''s SRID (%) is not the same as the raster''s SRID (%). The geometry will be transformed to the raster''s projection', g_srid, sr_id;
g := ST_Transform(geom, sr_id);
ELSE
g := geom;
END IF;
RETURN _st_asraster(g, scale_x, scale_y, NULL, NULL, $3, $4, $5, NULL, NULL, ul_x, ul_y, skew_x, skew_y);
END;
$$ LANGUAGE 'plpgsql' STABLE;
CREATE OR REPLACE FUNCTION st_asraster(
geom geometry,
ref raster,
pixeltype text,
value double precision DEFAULT 1,
nodataval double precision DEFAULT 0
)
RETURNS raster
AS $$ SELECT st_asraster($1, $2, ARRAY[$3]::text[], ARRAY[$4]::double precision[], ARRAY[$5]::double precision[]) $$
LANGUAGE 'sql' STABLE;
-----------------------------------------------------------------------
-- ST_Resample
-----------------------------------------------------------------------
@ -1453,7 +1624,7 @@ CREATE OR REPLACE FUNCTION _st_resample(
)
RETURNS raster
AS 'MODULE_PATHNAME', 'RASTER_resample'
LANGUAGE 'C' IMMUTABLE;
LANGUAGE 'C' STABLE;
CREATE OR REPLACE FUNCTION st_resample(
rast raster,
@ -1465,7 +1636,7 @@ CREATE OR REPLACE FUNCTION st_resample(
)
RETURNS raster
AS $$ SELECT _st_resample($1, $9, $10, $2, $3, $4, $5, $6, $7, $8) $$
LANGUAGE 'sql' IMMUTABLE;
LANGUAGE 'sql' STABLE;
CREATE OR REPLACE FUNCTION st_resample(
rast raster,
@ -1487,7 +1658,7 @@ CREATE OR REPLACE FUNCTION st_resample(
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);
END;
$$ LANGUAGE 'plpgsql' IMMUTABLE STRICT;
$$ LANGUAGE 'plpgsql' STABLE STRICT;
-----------------------------------------------------------------------
-- ST_Transform
@ -1495,17 +1666,17 @@ CREATE OR REPLACE FUNCTION st_resample(
CREATE OR REPLACE FUNCTION st_transform(rast raster, srid integer, algorithm text DEFAULT 'NearestNeighbour', maxerr double precision DEFAULT 0.125, scalex double precision DEFAULT 0, scaley double precision DEFAULT 0)
RETURNS raster
AS $$ SELECT _st_resample($1, $3, $4, $2, $5, $6) $$
LANGUAGE 'sql' IMMUTABLE STRICT;
LANGUAGE 'sql' STABLE STRICT;
CREATE OR REPLACE FUNCTION st_transform(rast raster, srid integer, scalex double precision, scaley double precision, algorithm text DEFAULT 'NearestNeighbour', maxerr double precision DEFAULT 0.125)
RETURNS raster
AS $$ SELECT _st_resample($1, $5, $6, $2, $3, $4) $$
LANGUAGE 'sql' IMMUTABLE STRICT;
LANGUAGE 'sql' STABLE STRICT;
CREATE OR REPLACE FUNCTION st_transform(rast raster, srid integer, scalexy double precision, algorithm text DEFAULT 'NearestNeighbour', maxerr double precision DEFAULT 0.125)
RETURNS raster
AS $$ SELECT _st_resample($1, $4, $5, $2, $3, $3) $$
LANGUAGE 'sql' IMMUTABLE STRICT;
LANGUAGE 'sql' STABLE STRICT;
-----------------------------------------------------------------------
-- ST_Rescale
@ -1513,12 +1684,12 @@ CREATE OR REPLACE FUNCTION st_transform(rast raster, srid integer, scalexy doubl
CREATE OR REPLACE FUNCTION st_rescale(rast raster, scalex double precision, scaley double precision, algorithm text DEFAULT 'NearestNeighbour', maxerr double precision DEFAULT 0.125)
RETURNS raster
AS $$ SELECT _st_resample($1, $4, $5, NULL, $2, $3) $$
LANGUAGE 'sql' IMMUTABLE STRICT;
LANGUAGE 'sql' STABLE STRICT;
CREATE OR REPLACE FUNCTION st_rescale(rast raster, scalexy double precision, algorithm text DEFAULT 'NearestNeighbour', maxerr double precision DEFAULT 0.125)
RETURNS raster
AS $$ SELECT _st_resample($1, $3, $4, NULL, $2, $2) $$
LANGUAGE 'sql' IMMUTABLE STRICT;
LANGUAGE 'sql' STABLE STRICT;
-----------------------------------------------------------------------
-- ST_Reskew
@ -1526,12 +1697,12 @@ CREATE OR REPLACE FUNCTION st_rescale(rast raster, scalexy double precision, alg
CREATE OR REPLACE FUNCTION st_reskew(rast raster, skewx double precision, skewy double precision, algorithm text DEFAULT 'NearestNeighbour', maxerr double precision DEFAULT 0.125)
RETURNS raster
AS $$ SELECT _st_resample($1, $4, $5, NULL, 0, 0, NULL, NULL, $2, $3) $$
LANGUAGE 'sql' IMMUTABLE STRICT;
LANGUAGE 'sql' STABLE STRICT;
CREATE OR REPLACE FUNCTION st_reskew(rast raster, skewxy double precision, algorithm text DEFAULT 'NearestNeighbour', maxerr double precision DEFAULT 0.125)
RETURNS raster
AS $$ SELECT _st_resample($1, $3, $4, NULL, 0, 0, NULL, NULL, $2, $2) $$
LANGUAGE 'sql' IMMUTABLE STRICT;
LANGUAGE 'sql' STABLE STRICT;
-----------------------------------------------------------------------
-- ST_SnapToGrid
@ -1539,17 +1710,17 @@ CREATE OR REPLACE FUNCTION st_reskew(rast raster, skewxy double precision, algor
CREATE OR REPLACE FUNCTION st_snaptogrid(rast raster, gridx double precision, gridy double precision, algorithm text DEFAULT 'NearestNeighbour', maxerr double precision DEFAULT 0.125, scalex double precision DEFAULT 0, scaley double precision DEFAULT 0)
RETURNS raster
AS $$ SELECT _st_resample($1, $4, $5, NULL, $6, $7, $2, $3) $$
LANGUAGE 'sql' IMMUTABLE STRICT;
LANGUAGE 'sql' STABLE STRICT;
CREATE OR REPLACE FUNCTION st_snaptogrid(rast raster, gridx double precision, gridy double precision, scalex double precision, scaley double precision, algorithm text DEFAULT 'NearestNeighbour', maxerr double precision DEFAULT 0.125)
RETURNS raster
AS $$ SELECT _st_resample($1, $6, $7, NULL, $4, $5, $2, $3) $$
LANGUAGE 'sql' IMMUTABLE STRICT;
LANGUAGE 'sql' STABLE STRICT;
CREATE OR REPLACE FUNCTION st_snaptogrid(rast raster, gridx double precision, gridy double precision, scalexy double precision, algorithm text DEFAULT 'NearestNeighbour', maxerr double precision DEFAULT 0.125)
RETURNS raster
AS $$ SELECT _st_resample($1, $5, $6, NULL, $4, $4, $2, $3) $$
LANGUAGE 'sql' IMMUTABLE STRICT;
LANGUAGE 'sql' STABLE STRICT;
-----------------------------------------------------------------------
-- MapAlgebra

View file

@ -1517,6 +1517,56 @@ static void testGDALWarp() {
deepRelease(raster);
}
static void testGDALRasterize() {
rt_raster raster;
char srs[] = "PROJCS[\"unnamed\",GEOGCS[\"unnamed ellipse\",DATUM[\"unknown\",SPHEROID[\"unnamed\",6370997,0]],PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433]],PROJECTION[\"Lambert_Azimuthal_Equal_Area\"],PARAMETER[\"latitude_of_center\",45],PARAMETER[\"longitude_of_center\",-100],PARAMETER[\"false_easting\",0],PARAMETER[\"false_northing\",0],UNIT[\"Meter\",1],AUTHORITY[\"EPSG\",\"2163\"]]";
const char wkb_hex[] = "010300000001000000050000000000000080841ec100000000600122410000000080841ec100000000804f22410000000040e81dc100000000804f22410000000040e81dc100000000600122410000000080841ec10000000060012241";
const char *pos = wkb_hex;
unsigned char *wkb = NULL;
int wkb_len = 0;
int i;
double scale_x = 100;
double scale_y = 100;
rt_pixtype pixtype[] = {PT_8BUI};
double init[] = {0};
double value[] = {1};
double nodata[] = {0};
uint8_t nodata_mask[] = {1};
/* hex to byte */
wkb_len = (int) ceil(((double) strlen(wkb_hex)) / 2);
wkb = (unsigned char *) malloc(sizeof(unsigned char) * wkb_len);
for (i = 0; i < wkb_len; i++) {
sscanf(pos, "%2hhx", &wkb[i]);
pos += 2;
}
raster = rt_raster_gdal_rasterize(
wkb,
wkb_len, srs,
1, pixtype,
init, value,
nodata, nodata_mask,
NULL, NULL,
&scale_x, &scale_y,
NULL, NULL,
NULL, NULL,
NULL, NULL
);
free(wkb);
CHECK(raster);
CHECK((rt_raster_get_width(raster) == 100));
CHECK((rt_raster_get_height(raster) == 100));
CHECK((rt_raster_get_num_bands(raster) != 0));
CHECK((rt_raster_get_x_offset(raster) == -500000));
CHECK((rt_raster_get_y_offset(raster) == 600000));
deepRelease(raster);
}
int
main()
{
@ -1907,6 +1957,10 @@ main()
testGDALWarp();
printf("Successfully tested rt_raster_gdal_warp\n");
printf("Testing rt_raster_gdal_rasterize\n");
testGDALRasterize();
printf("Successfully tested rt_raster_gdal_rasterize\n");
deepRelease(raster);
return EXIT_SUCCESS;

View file

@ -88,6 +88,7 @@ TEST_UTILITY = \
rt_mapalgebra.sql \
rt_reclass.sql \
rt_resample.sql \
rt_asraster.sql \
$(NULL)
TEST_GIST = \

View file

@ -0,0 +1,468 @@
DROP TABLE IF EXISTS raster_asraster_geom;
DROP TABLE IF EXISTS raster_asraster_rast;
DROP TABLE IF EXISTS raster_asraster_dst;
CREATE TABLE raster_asraster_geom (
geom geometry
);
CREATE TABLE raster_asraster_rast (
rast raster
);
CREATE TABLE raster_asraster_dst (
rid varchar,
rast raster
);
CREATE OR REPLACE FUNCTION make_test_raster()
RETURNS void
AS $$
DECLARE
width int := 10;
height int := 10;
x int;
y int;
rast raster;
BEGIN
rast := ST_MakeEmptyRaster(width, height, -500000, 600000, 1000, -1000, 0, 0, 992163);
rast := ST_AddBand(rast, 1, '64BF', 0, 0);
FOR x IN 1..width LOOP
FOR y IN 1..height LOOP
rast := ST_SetValue(rast, 1, x, y, ((x::double precision * y) + (x + y) + (x + y * x)) / (x + y + 1));
END LOOP;
END LOOP;
INSERT INTO raster_asraster_rast VALUES (rast);
RETURN;
END;
$$ LANGUAGE 'plpgsql';
SELECT make_test_raster();
DROP FUNCTION make_test_raster();
DELETE FROM "spatial_ref_sys" WHERE srid = 992163;
DELETE FROM "spatial_ref_sys" WHERE srid = 993309;
DELETE FROM "spatial_ref_sys" WHERE srid = 993310;
DELETE FROM "spatial_ref_sys" WHERE srid = 994269;
INSERT INTO "spatial_ref_sys" ("srid","auth_name","auth_srid","srtext","proj4text") VALUES (992163,'EPSG',2163,'PROJCS["unnamed",GEOGCS["unnamed ellipse",DATUM["unknown",SPHEROID["unnamed",6370997,0]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]],PROJECTION["Lambert_Azimuthal_Equal_Area"],PARAMETER["latitude_of_center",45],PARAMETER["longitude_of_center",-100],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["Meter",1],AUTHORITY["EPSG","2163"]]','+proj=laea +lat_0=45 +lon_0=-100 +x_0=0 +y_0=0 +a=6370997 +b=6370997 +units=m +no_defs ');
INSERT INTO "spatial_ref_sys" ("srid","auth_name","auth_srid","srtext","proj4text") VALUES (993309,'EPSG',3309,'PROJCS["NAD27 / California Albers",GEOGCS["NAD27",DATUM["North_American_Datum_1927",SPHEROID["Clarke 1866",6378206.4,294.9786982139006,AUTHORITY["EPSG","7008"]],AUTHORITY["EPSG","6267"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4267"]],UNIT["metre",1,AUTHORITY["EPSG","9001"]],PROJECTION["Albers_Conic_Equal_Area"],PARAMETER["standard_parallel_1",34],PARAMETER["standard_parallel_2",40.5],PARAMETER["latitude_of_center",0],PARAMETER["longitude_of_center",-120],PARAMETER["false_easting",0],PARAMETER["false_northing",-4000000],AUTHORITY["EPSG","3309"],AXIS["X",EAST],AXIS["Y",NORTH]]','+proj=aea +lat_1=34 +lat_2=40.5 +lat_0=0 +lon_0=-120 +x_0=0 +y_0=-4000000 +ellps=clrk66 +datum=NAD27 +units=m +no_defs ');
INSERT INTO "spatial_ref_sys" ("srid","auth_name","auth_srid","srtext","proj4text") VALUES (993310,'EPSG',3310,'PROJCS["NAD83 / California Albers",GEOGCS["NAD83",DATUM["North_American_Datum_1983",SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],AUTHORITY["EPSG","6269"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4269"]],UNIT["metre",1,AUTHORITY["EPSG","9001"]],PROJECTION["Albers_Conic_Equal_Area"],PARAMETER["standard_parallel_1",34],PARAMETER["standard_parallel_2",40.5],PARAMETER["latitude_of_center",0],PARAMETER["longitude_of_center",-120],PARAMETER["false_easting",0],PARAMETER["false_northing",-4000000],AUTHORITY["EPSG","3310"],AXIS["X",EAST],AXIS["Y",NORTH]]','+proj=aea +lat_1=34 +lat_2=40.5 +lat_0=0 +lon_0=-120 +x_0=0 +y_0=-4000000 +ellps=GRS80 +datum=NAD83 +units=m +no_defs ');
INSERT INTO "spatial_ref_sys" ("srid","auth_name","auth_srid","srtext","proj4text") VALUES (994269,'EPSG',4269,'GEOGCS["NAD83",DATUM["North_American_Datum_1983",SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],AUTHORITY["EPSG","6269"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4269"]]','+proj=longlat +ellps=GRS80 +datum=NAD83 +no_defs ');
INSERT INTO raster_asraster_geom VALUES (
ST_GeomFromText('MULTIPOLYGON(((-172210.499109288 114987.660953018,-175453.086381862 29201.5994536821,-151944.038528546 28257.4637483698,-151755.193144738 64618.6592845297,-129779.244489461 63766.2346307114,-132720.730482521 29365.7452160339,-110176.183408147 28076.2457866343,-113336.283431208 112064.985603184,-135659.619600536 112878.300914729,-134301.95687566 79576.8821948012,-153850.618867315 80395.4252778995,-151346.215838074 112678.410158427,-172210.499109288 114987.660953018)),((-86135.5150847774 77502.7616508612,-87105.1850870571 30678.0039829779,-69362.3449961895 29072.3373203999,-70858.5814585458 78310.0439012805,-86135.5150847774 77502.7616508612)),((-86888.5102830273 96546.8546876945,-86065.7795470885 84169.9977753228,-70801.2145468401 84976.5822106288,-72118.6159803197 97829.7405064492,-86888.5102830273 96546.8546876945)),((-50136.8809020698 111909.445130098,-48631.3614059008 44728.8885465469,-36172.0195739627 45621.806341459,-39695.018962698 109480.225649309,-50136.8809020698 111909.445130098)),((-47695.3501850868 40894.9976787795,-47761.6362577873 29399.0052930373,-34799.4262271112 30293.0638067261,-35717.8219710071 39877.2161100041,-47695.3501850868 40894.9976787795)))', 993310)
);
-- scale or width & height, pixtype, value and nodata
INSERT INTO raster_asraster_dst (rid, rast) VALUES (
1.1, (SELECT ST_AsRaster(
geom,
100, 100
) FROM raster_asraster_geom)
), (
1.2, (SELECT ST_AsRaster(
geom,
100., 100.
) FROM raster_asraster_geom)
), (
1.3, (SELECT ST_AsRaster(
geom,
500, 500
) FROM raster_asraster_geom)
), (
1.4, (SELECT ST_AsRaster(
geom,
1000., 1000.
) FROM raster_asraster_geom)
), (
1.5, (SELECT ST_AsRaster(
geom,
1000., 1000.,
'8BSI'
) FROM raster_asraster_geom)
), (
1.6, (SELECT ST_AsRaster(
geom,
1000., 1000.,
'16BUI'
) FROM raster_asraster_geom)
), (
1.7, (SELECT ST_AsRaster(
geom,
100., 100.,
'32BF'
) FROM raster_asraster_geom)
), (
1.8, (SELECT ST_AsRaster(
geom,
1000., 1000.,
ARRAY['8BSI']
) FROM raster_asraster_geom)
), (
1.9, (SELECT ST_AsRaster(
geom,
1000., 1000.,
ARRAY['16BUI']
) FROM raster_asraster_geom)
), (
1.10, (SELECT ST_AsRaster(
geom,
1000., 1000.,
ARRAY['32BF']
) FROM raster_asraster_geom)
), (
1.11, (SELECT ST_AsRaster(
geom,
100, 100,
ARRAY['8BSI']
) FROM raster_asraster_geom)
), (
1.12, (SELECT ST_AsRaster(
geom,
100, 100,
'16BUI'
) FROM raster_asraster_geom)
), (
1.13, (SELECT ST_AsRaster(
geom,
100, 100,
ARRAY['32BF'],
ARRAY[255]
) FROM raster_asraster_geom)
), (
1.14, (SELECT ST_AsRaster(
geom,
100, 100,
ARRAY['32BF'],
ARRAY[255],
ARRAY[1]
) FROM raster_asraster_geom)
), (
1.15, (SELECT ST_AsRaster(
geom,
100, 100,
ARRAY['32BF'],
ARRAY[255],
NULL
) FROM raster_asraster_geom)
), (
1.16, (SELECT ST_AsRaster(
geom,
100, 100,
ARRAY['32BF'],
ARRAY[255],
ARRAY[NULL]::double precision[]
) FROM raster_asraster_geom)
), (
1.17, (SELECT ST_AsRaster(
geom,
1000., 1000.,
ARRAY['32BF', '16BUI'],
ARRAY[255, 1],
ARRAY[NULL, 0]::double precision[]
) FROM raster_asraster_geom)
), (
1.18, (SELECT ST_AsRaster(
geom,
10, 10,
ARRAY['8BUI', '16BUI'],
ARRAY[255, 255],
ARRAY[0, NULL]::double precision[]
) FROM raster_asraster_geom)
), (
1.19, (SELECT ST_AsRaster(
geom,
1000., 1000.,
ARRAY['32BF', '16BUI', '64BF'],
ARRAY[255, 1, -1],
ARRAY[NULL, 0, NULL]::double precision[]
) FROM raster_asraster_geom)
), (
1.20, (SELECT ST_AsRaster(
geom,
1000., 1000.,
ARRAY['1BB', '2BUI'],
ARRAY[1, 1],
ARRAY[1, 0]::double precision[]
) FROM raster_asraster_geom)
);
-- upper left
INSERT INTO raster_asraster_dst (rid, rast) VALUES (
2.1, (SELECT ST_AsRaster(
geom,
1000., 1000.,
'8BUI',
255,
0,
-175453
) FROM raster_asraster_geom)
), (
2.2, (SELECT ST_AsRaster(
geom,
1000., 1000.,
'8BUI',
255,
0,
-175400, 115000
) FROM raster_asraster_geom)
), (
2.3, (SELECT ST_AsRaster(
geom,
1000., 1000.,
'8BUI',
255,
0,
-170000, 114988
) FROM raster_asraster_geom)
), (
2.4, (SELECT ST_AsRaster(
geom,
1000., 1000.,
'8BUI',
255,
0,
-170000, 110000
) FROM raster_asraster_geom)
), (
2.5, (SELECT ST_AsRaster(
geom,
1000., 1000.,
'8BUI',
255,
0,
-179000, 119000
) FROM raster_asraster_geom)
), (
2.6, (SELECT ST_AsRaster(
geom,
100, 100,
'8BUI',
255,
0,
-179000, 119000
) FROM raster_asraster_geom)
), (
2.7, (SELECT ST_AsRaster(
geom,
100, 100,
ARRAY['8BUI'],
ARRAY[255],
ARRAY[0],
-179000, 119000
) FROM raster_asraster_geom)
);
-- skew
INSERT INTO raster_asraster_dst (rid, rast) VALUES (
3.1, (SELECT ST_AsRaster(
geom,
100, 100,
'8BUI',
255,
0,
NULL, NULL,
0
) FROM raster_asraster_geom)
), (
3.2, (SELECT ST_AsRaster(
geom,
100, 100,
'8BUI',
255,
0,
NULL, NULL,
0, 0
) FROM raster_asraster_geom)
), (
3.3, (SELECT ST_AsRaster(
geom,
100, 100,
'8BUI',
255,
0,
NULL, NULL,
1, 0
) FROM raster_asraster_geom)
), (
3.4, (SELECT ST_AsRaster(
geom,
100, 100,
'8BUI',
255,
0,
NULL, NULL,
0, 1
) FROM raster_asraster_geom)
), (
3.5, (SELECT ST_AsRaster(
geom,
100, 100,
'8BUI',
255,
0,
NULL, NULL,
10, -5
) FROM raster_asraster_geom)
), (
3.6, (SELECT ST_AsRaster(
geom,
100, 100,
'8BUI',
255,
0,
NULL, NULL,
-5, 10
) FROM raster_asraster_geom)
);
-- snap to grid
INSERT INTO raster_asraster_dst (rid, rast) VALUES (
4.1, (
SELECT ST_AsRaster(
geom,
rast
)
FROM raster_asraster_geom, raster_asraster_rast
)
), (
4.2, (
SELECT ST_AsRaster(
geom,
rast,
'64BF'
)
FROM raster_asraster_geom, raster_asraster_rast
)
), (
4.3, (
SELECT ST_AsRaster(
geom,
rast,
'16BUI',
13
)
FROM raster_asraster_geom, raster_asraster_rast
)
), (
4.4, (
SELECT ST_AsRaster(
geom,
rast,
'16BUI',
13,
NULL
)
FROM raster_asraster_geom, raster_asraster_rast
)
), (
4.5, (
SELECT ST_AsRaster(
geom,
rast,
ARRAY['16BUI'],
ARRAY[13]
)
FROM raster_asraster_geom, raster_asraster_rast
)
), (
4.6, (
SELECT ST_AsRaster(
geom,
rast,
ARRAY['16BUI'],
ARRAY[13],
ARRAY[NULL]::double precision[]
)
FROM raster_asraster_geom, raster_asraster_rast
)
), (
4.7, (
SELECT ST_AsRaster(
geom,
rast,
ARRAY['16BUI'],
ARRAY[13],
ARRAY[0]
)
FROM raster_asraster_geom, raster_asraster_rast
)
), (
4.8, (SELECT ST_AsRaster(
geom,
1000., 1000.,
0, 0,
ARRAY['16BUI'],
ARRAY[13],
ARRAY[0]
)
FROM raster_asraster_geom)
), (
4.9, (SELECT ST_AsRaster(
geom,
1000., 1000.,
-175453, 114987,
ARRAY['16BUI'],
ARRAY[13],
ARRAY[0]
)
FROM raster_asraster_geom)
), (
4.10, (SELECT ST_AsRaster(
geom,
1000., 1000.,
-100, 100,
ARRAY['16BUI'],
ARRAY[13],
ARRAY[0]
)
FROM raster_asraster_geom)
), (
4.11, (SELECT ST_AsRaster(
geom,
1000., 1000.,
-100, 100,
'16BUI',
13,
0
)
FROM raster_asraster_geom)
);
SELECT
rid,
srid,
width,
height,
numbands,
round(scalex::numeric, 3) AS scalex,
round(scaley::numeric, 3) AS scaley,
round(skewx::numeric, 3) AS skewx,
round(skewy::numeric, 3) AS skewy,
round(upperleftx::numeric, 3) AS upperleftx,
round(upperlefty::numeric, 3) AS upperlefty,
pixeltype,
hasnodatavalue,
round(nodatavalue::numeric, 3) AS nodatavalue,
count > 0 AS count_check,
round(min::numeric, 3) AS min,
round(max::numeric, 3) AS max
FROM (
SELECT
rid,
(ST_MetaData(rast)).*,
(ST_SummaryStats(rast)).*,
(ST_BandMetaData(rast)).*
FROM raster_asraster_dst
ORDER BY rid
) foo;
DELETE FROM "spatial_ref_sys" WHERE srid = 992163;
DELETE FROM "spatial_ref_sys" WHERE srid = 993309;
DELETE FROM "spatial_ref_sys" WHERE srid = 993310;
DELETE FROM "spatial_ref_sys" WHERE srid = 994269;
DROP TABLE raster_asraster_geom;
DROP TABLE raster_asraster_rast;
DROP TABLE raster_asraster_dst;

View file

@ -0,0 +1,62 @@
NOTICE: table "raster_asraster_geom" does not exist, skipping
NOTICE: table "raster_asraster_rast" does not exist, skipping
NOTICE: table "raster_asraster_dst" does not exist, skipping
NOTICE: Imbalanced number of values provided for pixeltype (1), value (1) and nodataval (0). Using the first 0 values of each parameter
NOTICE: Values must be provided for both X and Y when specifying the upper-left corner
NOTICE: The geometry's SRID (993310) is not the same as the raster's SRID (992163). The geometry will be transformed to the raster's projection
NOTICE: The geometry's SRID (993310) is not the same as the raster's SRID (992163). The geometry will be transformed to the raster's projection
NOTICE: The geometry's SRID (993310) is not the same as the raster's SRID (992163). The geometry will be transformed to the raster's projection
NOTICE: The geometry's SRID (993310) is not the same as the raster's SRID (992163). The geometry will be transformed to the raster's projection
NOTICE: The geometry's SRID (993310) is not the same as the raster's SRID (992163). The geometry will be transformed to the raster's projection
NOTICE: The geometry's SRID (993310) is not the same as the raster's SRID (992163). The geometry will be transformed to the raster's projection
NOTICE: The geometry's SRID (993310) is not the same as the raster's SRID (992163). The geometry will be transformed to the raster's projection
NOTICE: Could not retrieve summary statistics of band of index 1. Returning NULL
NOTICE: Could not retrieve summary statistics of band of index 1. Returning NULL
NOTICE: Could not retrieve summary statistics of band of index 1. Returning NULL
NOTICE: Could not retrieve summary statistics of band of index 1. Returning NULL
NOTICE: Could not retrieve summary statistics of band of index 1. Returning NULL
NOTICE: Could not retrieve summary statistics of band of index 1. Returning NULL
1.1|993310|100|100|1|1406.537|-869.114|0.000|0.000|-175453.086|114987.661|8BUI|t|0.000|t|1.000|1.000
1.10|993310|141|87|1|1000.000|-1000.000|0.000|0.000|-175453.086|114987.661|32BF|t|0.000|t|1.000|1.000
1.11|993310|100|100|1|1406.537|-869.114|0.000|0.000|-175453.086|114987.661|8BUI|t|0.000|t|1.000|1.000
1.12|993310|100|100|1|1406.537|-869.114|0.000|0.000|-175453.086|114987.661|16BUI|t|0.000|t|1.000|1.000
1.13|993310|100|100|1|1406.537|-869.114|0.000|0.000|-175453.086|114987.661|32BF|t|0.000|t|255.000|255.000
1.14|993310|100|100|1|1406.537|-869.114|0.000|0.000|-175453.086|114987.661|32BF|t|0.000|t|255.000|255.000
1.15|993310|100|100|1|1406.537|-869.114|0.000|0.000|-175453.086|114987.661|8BUI|t|0.000|t|1.000|1.000
1.16|993310|100|100|1|1406.537|-869.114|0.000|0.000|-175453.086|114987.661|32BF|f|0.000|t|0.000|255.000
1.17|993310|141|87|2|1000.000|-1000.000|0.000|0.000|-175453.086|114987.661|32BF|f|0.000|t|0.000|255.000
1.18|993310|10|10|2|14065.366|-8691.142|0.000|0.000|-175453.086|114987.661|8BUI|t|0.000|t|255.000|255.000
1.19|993310|141|87|3|1000.000|-1000.000|0.000|0.000|-175453.086|114987.661|32BF|f|0.000|t|0.000|255.000
1.2|993310|1407|869|1|100.000|-100.000|0.000|0.000|-175453.086|114987.661|8BUI|t|0.000|t|1.000|1.000
1.20|993310|141|87|2|1000.000|-1000.000|0.000|0.000|-175453.086|114987.661|8BUI|t|0.000|||
1.3|993310|500|500|1|281.307|-173.823|0.000|0.000|-175453.086|114987.661|8BUI|t|0.000|t|1.000|1.000
1.4|993310|141|87|1|1000.000|-1000.000|0.000|0.000|-175453.086|114987.661|8BUI|t|0.000|t|1.000|1.000
1.5|993310|141|87|1|1000.000|-1000.000|0.000|0.000|-175453.086|114987.661|8BUI|t|0.000|t|1.000|1.000
1.6|993310|141|87|1|1000.000|-1000.000|0.000|0.000|-175453.086|114987.661|16BUI|t|0.000|t|1.000|1.000
1.7|993310|1407|869|1|100.000|-100.000|0.000|0.000|-175453.086|114987.661|32BF|t|0.000|t|1.000|1.000
1.8|993310|141|87|1|1000.000|-1000.000|0.000|0.000|-175453.086|114987.661|8BUI|t|0.000|t|1.000|1.000
1.9|993310|141|87|1|1000.000|-1000.000|0.000|0.000|-175453.086|114987.661|16BUI|t|0.000|t|1.000|1.000
2.1||||||||||||||||
2.2|993310|141|87|1|1000.000|-1000.000|0.000|0.000|-175400.000|115000.000|8BUI|t|0.000|t|255.000|255.000
2.3|993310|135|87|1|1000.000|-1000.000|0.000|0.000|-170000.000|114988.000|8BUI|t|0.000|t|255.000|255.000
2.4|993310|135|82|1|1000.000|-1000.000|0.000|0.000|-170000.000|110000.000|8BUI|t|0.000|t|255.000|255.000
2.5|993310|144|91|1|1000.000|-1000.000|0.000|0.000|-179000.000|119000.000|8BUI|t|0.000|t|255.000|255.000
2.6|993310|100|100|1|1406.537|-869.114|0.000|0.000|-179000.000|119000.000|8BUI|t|0.000|t|255.000|255.000
2.7|993310|100|100|1|1406.537|-869.114|0.000|0.000|-179000.000|119000.000|8BUI|t|0.000|t|255.000|255.000
3.1|993310|100|100|1|1406.537|-869.114|0.000|0.000|-175453.086|114987.661|8BUI|t|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|t|0.000|t|255.000|255.000
3.3|993310|100|100|1|1406.537|-869.114|1.000|0.000|-175453.086|114987.661|8BUI|t|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|t|0.000|t|255.000|255.000
3.5|993310|100|100|1|1406.537|-869.114|10.000|-5.000|-175453.086|114987.661|8BUI|t|0.000|t|255.000|255.000
3.6|993310|100|100|1|1406.537|-869.114|-5.000|10.000|-175453.086|114987.661|8BUI|t|0.000|t|255.000|255.000
4.1|992163|150|116|1|1000.000|-1000.000|0.000|0.000|-1898000.000|-412000.000|8BUI|t|0.000|t|1.000|1.000
4.10|993310|141|87|1|1000.000|-1000.000|0.000|0.000|-176100.000|115100.000|16BUI|t|0.000|t|13.000|13.000
4.11|993310|141|87|1|1000.000|-1000.000|0.000|0.000|-176100.000|115100.000|16BUI|t|0.000|t|13.000|13.000
4.2|992163|150|116|1|1000.000|-1000.000|0.000|0.000|-1898000.000|-412000.000|64BF|t|0.000|t|1.000|1.000
4.3|992163|150|116|1|1000.000|-1000.000|0.000|0.000|-1898000.000|-412000.000|16BUI|t|0.000|t|13.000|13.000
4.4|992163|150|116|1|1000.000|-1000.000|0.000|0.000|-1898000.000|-412000.000|16BUI|f|0.000|t|0.000|13.000
4.5|992163|150|116|1|1000.000|-1000.000|0.000|0.000|-1898000.000|-412000.000|16BUI|t|0.000|t|13.000|13.000
4.6|992163|150|116|1|1000.000|-1000.000|0.000|0.000|-1898000.000|-412000.000|16BUI|f|0.000|t|0.000|13.000
4.7|992163|150|116|1|1000.000|-1000.000|0.000|0.000|-1898000.000|-412000.000|16BUI|t|0.000|t|13.000|13.000
4.8|993310|141|87|1|1000.000|-1000.000|0.000|0.000|-176000.000|115000.000|16BUI|t|0.000|t|13.000|13.000
4.9|993310|142|88|1|1000.000|-1000.000|0.000|0.000|-176453.000|115987.000|16BUI|t|0.000|t|13.000|13.000

View file

@ -7,6 +7,7 @@ CREATE TABLE raster_resample_dst (
rid varchar,
rast raster
);
CREATE OR REPLACE FUNCTION make_test_raster()
RETURNS void
AS $$
@ -32,10 +33,13 @@ CREATE OR REPLACE FUNCTION make_test_raster()
END;
$$ LANGUAGE 'plpgsql';
SELECT make_test_raster();
DROP FUNCTION make_test_raster();
DELETE FROM "spatial_ref_sys" WHERE srid = 992163;
DELETE FROM "spatial_ref_sys" WHERE srid = 993309;
DELETE FROM "spatial_ref_sys" WHERE srid = 993310;
DELETE FROM "spatial_ref_sys" WHERE srid = 994269;
INSERT INTO "spatial_ref_sys" ("srid","auth_name","auth_srid","srtext","proj4text") VALUES (992163,'EPSG',2163,'PROJCS["unnamed",GEOGCS["unnamed ellipse",DATUM["unknown",SPHEROID["unnamed",6370997,0]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]],PROJECTION["Lambert_Azimuthal_Equal_Area"],PARAMETER["latitude_of_center",45],PARAMETER["longitude_of_center",-100],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["Meter",1],AUTHORITY["EPSG","2163"]]','+proj=laea +lat_0=45 +lon_0=-100 +x_0=0 +y_0=0 +a=6370997 +b=6370997 +units=m +no_defs ');
INSERT INTO "spatial_ref_sys" ("srid","auth_name","auth_srid","srtext","proj4text") VALUES (993309,'EPSG',3309,'PROJCS["NAD27 / California Albers",GEOGCS["NAD27",DATUM["North_American_Datum_1927",SPHEROID["Clarke 1866",6378206.4,294.9786982139006,AUTHORITY["EPSG","7008"]],AUTHORITY["EPSG","6267"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4267"]],UNIT["metre",1,AUTHORITY["EPSG","9001"]],PROJECTION["Albers_Conic_Equal_Area"],PARAMETER["standard_parallel_1",34],PARAMETER["standard_parallel_2",40.5],PARAMETER["latitude_of_center",0],PARAMETER["longitude_of_center",-120],PARAMETER["false_easting",0],PARAMETER["false_northing",-4000000],AUTHORITY["EPSG","3309"],AXIS["X",EAST],AXIS["Y",NORTH]]','+proj=aea +lat_1=34 +lat_2=40.5 +lat_0=0 +lon_0=-120 +x_0=0 +y_0=-4000000 +ellps=clrk66 +datum=NAD27 +units=m +no_defs ');
INSERT INTO "spatial_ref_sys" ("srid","auth_name","auth_srid","srtext","proj4text") VALUES (993310,'EPSG',3310,'PROJCS["NAD83 / California Albers",GEOGCS["NAD83",DATUM["North_American_Datum_1983",SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],AUTHORITY["EPSG","6269"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4269"]],UNIT["metre",1,AUTHORITY["EPSG","9001"]],PROJECTION["Albers_Conic_Equal_Area"],PARAMETER["standard_parallel_1",34],PARAMETER["standard_parallel_2",40.5],PARAMETER["latitude_of_center",0],PARAMETER["longitude_of_center",-120],PARAMETER["false_easting",0],PARAMETER["false_northing",-4000000],AUTHORITY["EPSG","3310"],AXIS["X",EAST],AXIS["Y",NORTH]]','+proj=aea +lat_1=34 +lat_2=40.5 +lat_0=0 +lon_0=-120 +x_0=0 +y_0=-4000000 +ellps=GRS80 +datum=NAD83 +units=m +no_defs ');
@ -526,7 +530,7 @@ SELECT
round(skewy::numeric, 3) AS skewy,
round(upperleftx::numeric, 3) AS upperleftx,
round(upperlefty::numeric, 3) AS upperlefty,
count > 0,
count > 0 AS count_check,
round(min::numeric, 3) >= 1.667 AS min_check,
round(max::numeric, 3) <= 100.995 AS max_check
FROM (
@ -537,10 +541,11 @@ FROM (
FROM raster_resample_dst
ORDER BY rid
) foo;
DELETE FROM "spatial_ref_sys" WHERE srid = 992163;
DELETE FROM "spatial_ref_sys" WHERE srid = 993309;
DELETE FROM "spatial_ref_sys" WHERE srid = 993310;
DELETE FROM "spatial_ref_sys" WHERE srid = 994269;
DROP TABLE raster_resample_src;
DROP TABLE raster_resample_dst;
DROP FUNCTION make_test_raster();