diff --git a/raster/rt_core/rt_api.c b/raster/rt_core/rt_api.c
index 13d7a4fed..0a6a220ee 100644
--- a/raster/rt_core/rt_api.c
+++ b/raster/rt_core/rt_api.c
@@ -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;
+}
diff --git a/raster/rt_core/rt_api.h b/raster/rt_core/rt_api.h
index d52d20de5..0a92d3746 100644
--- a/raster/rt_core/rt_api.h
+++ b/raster/rt_core/rt_api.h
@@ -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 -------------------------------------------------------*/
 
 /*
diff --git a/raster/rt_pg/rt_pg.c b/raster/rt_pg/rt_pg.c
index b39015d94..3399bf6fd 100644
--- a/raster/rt_pg/rt_pg.c
+++ b/raster/rt_pg/rt_pg.c
@@ -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
  */
diff --git a/raster/rt_pg/rtpostgis.sql.in.c b/raster/rt_pg/rtpostgis.sql.in.c
index 083c77b95..f16d5f2c4 100644
--- a/raster/rt_pg/rtpostgis.sql.in.c
+++ b/raster/rt_pg/rtpostgis.sql.in.c
@@ -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
diff --git a/raster/test/core/testapi.c b/raster/test/core/testapi.c
index 183de483c..448e97bfc 100644
--- a/raster/test/core/testapi.c
+++ b/raster/test/core/testapi.c
@@ -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;
diff --git a/raster/test/regress/Makefile.in b/raster/test/regress/Makefile.in
index 848112863..2a86cbe5f 100644
--- a/raster/test/regress/Makefile.in
+++ b/raster/test/regress/Makefile.in
@@ -88,6 +88,7 @@ TEST_UTILITY = \
 	rt_mapalgebra.sql \
 	rt_reclass.sql \
 	rt_resample.sql \
+	rt_asraster.sql \
 	$(NULL)
 
 TEST_GIST = \
diff --git a/raster/test/regress/rt_asraster.sql b/raster/test/regress/rt_asraster.sql
new file mode 100644
index 000000000..a74174aab
--- /dev/null
+++ b/raster/test/regress/rt_asraster.sql
@@ -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;
diff --git a/raster/test/regress/rt_asraster_expected b/raster/test/regress/rt_asraster_expected
new file mode 100644
index 000000000..36f5cf10e
--- /dev/null
+++ b/raster/test/regress/rt_asraster_expected
@@ -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
diff --git a/raster/test/regress/rt_resample.sql b/raster/test/regress/rt_resample.sql
index d755cdd69..d76cf7757 100644
--- a/raster/test/regress/rt_resample.sql
+++ b/raster/test/regress/rt_resample.sql
@@ -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();