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;
+}
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 -------------------------------------------------------*/
/*
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);
}
}
+/**
+ * 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
*/
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
-----------------------------------------------------------------------
)
RETURNS raster
AS 'MODULE_PATHNAME', 'RASTER_resample'
- LANGUAGE 'C' IMMUTABLE;
+ LANGUAGE 'C' STABLE;
CREATE OR REPLACE FUNCTION st_resample(
rast raster,
)
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,
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
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
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
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
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
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()
{
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;
rt_mapalgebra.sql \
rt_reclass.sql \
rt_resample.sql \
+ rt_asraster.sql \
$(NULL)
TEST_GIST = \
--- /dev/null
+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;
--- /dev/null
+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
rid varchar,
rast raster
);
+
CREATE OR REPLACE FUNCTION make_test_raster()
RETURNS void
AS $$
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 ');
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 (
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();
projects / postgis / commitdiff