From: Bborie Park Date: Mon, 7 Nov 2011 21:07:56 +0000 (+0000) Subject: Addition of 2-raster map algebra function ST_MapAlgebra2Expr. Next is ST_MapAlgebra2Fct. X-Git-Tag: 2.0.0alpha1~752 X-Git-Url: https://granicus.if.org/sourcecode?a=commitdiff_plain;h=c395e3cfe2cfa9aa362dfa828fca570d47fe0943;p=postgis Addition of 2-raster map algebra function ST_MapAlgebra2Expr. Next is ST_MapAlgebra2Fct. Ticket is #1268. git-svn-id: http://svn.osgeo.org/postgis/trunk@8108 b70326c6-7e19-0410-871a-916f4a2858ee --- diff --git a/raster/rt_pg/rt_pg.c b/raster/rt_pg/rt_pg.c index f3977cc8b..07143e16a 100644 --- a/raster/rt_pg/rt_pg.c +++ b/raster/rt_pg/rt_pg.c @@ -243,6 +243,12 @@ Datum RASTER_intersects(PG_FUNCTION_ARGS); /* determine if two rasters are aligned */ Datum RASTER_sameAlignment(PG_FUNCTION_ARGS); +/* two-raster MapAlgebra */ +Datum RASTER_mapAlgebra2Expr(PG_FUNCTION_ARGS); +/* +Datum RASTER_mapAlgebra2Fct(PG_FUNCTION_ARGS); +*/ + /* Replace function taken from * http://ubuntuforums.org/showthread.php?s=aa6f015109fd7e4c7e30d2fd8b717497&t=141670&page=3 */ @@ -7957,6 +7963,894 @@ Datum RASTER_sameAlignment(PG_FUNCTION_ARGS) PG_RETURN_BOOL(aligned); } +/** + * Two raster MapAlgebra + */ +PG_FUNCTION_INFO_V1(RASTER_mapAlgebra2Expr); +Datum RASTER_mapAlgebra2Expr(PG_FUNCTION_ARGS) +{ + const int set_count = 2; + rt_pgraster *pgrast; + rt_raster rast[2] = {NULL}; + int _isempty[2] = {0}; + uint32_t bandindex[2] = {0}; + rt_raster _rast[2] = {NULL}; + rt_band _band[2] = {NULL}; + + char *pixtypename = NULL; + rt_pixtype pixtype = PT_END; + char *extenttypename = NULL; + rt_extenttype extenttype = ET_INTERSECTION; + double nodataval = 0; + int _hasnodata[2] = {0}; + double _nodataval[2] = {0}; + + int hasnodatanodataval = 0; + double nodatanodataval = 0; + + int spicount = 3; + uint16_t exprpos[3] = {4, 7, 8}; + char *expr = NULL; + char *sql = NULL; + SPIPlanPtr spiplan[3] = {NULL}; + uint16_t spiempty = 0; + Oid *argtype = NULL; + uint32_t argcount[3] = {0}; + int argexists[3][2] = {{0}}; + Datum *values = NULL; + bool *nulls = NULL; + TupleDesc tupdesc; + SPITupleTable *tuptable = NULL; + HeapTuple tuple; + Datum pixeldatum; + bool isnull = FALSE; + + double _rastoffset[2][4] = {{0.}}; + int _haspixel[2] = {0}; + double _pixel[2] = {0}; + uint16_t _dim[2][2] = {{0}}; + + rt_raster raster = NULL; + rt_band band = NULL; + uint16_t dim[2] = {0}; + int haspixel = 0; + double pixel = 0.; + + uint32_t i = 0; + uint32_t j = 0; + uint32_t k = 0; + uint32_t x = 0; + uint32_t y = 0; + int _x = 0; + int _y = 0; + int err; + double gt[6] = {0.}; + int aligned; + int len; + + POSTGIS_RT_DEBUG(3, "Starting RASTER_mapAlgebra2Expr"); + + for (i = 0, j = 0; i < set_count; i++) { + if (!PG_ARGISNULL(j)) { + pgrast = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(j)); + j++; + + /* raster */ + rast[i] = rt_raster_deserialize(pgrast, FALSE); + if (!rast[i]) { + elog(ERROR, "RASTER_mapAlgebra2Expr: Could not deserialize the %s raster", i < 1 ? "first" : "second"); + for (k = 0; k < i; k++) { + if (rast[k] != NULL) rt_raster_destroy(rast[k]); + } + PG_RETURN_NULL(); + } + + /* empty */ + _isempty[i] = rt_raster_is_empty(rast[i]); + + /* band index */ + if (!PG_ARGISNULL(j)) { + bandindex[i] = PG_GETARG_INT32(j); + } + j++; + } + else { + _isempty[i] = 1; + j += 2; + } + + POSTGIS_RT_DEBUGF(3, "_isempty[%d] = %d", i, _isempty[i]); + } + + /* both rasters are NULL */ + if (rast[0] == NULL && rast[1] == NULL) { + elog(NOTICE, "The two rasters provided are NULL. Returning NULL"); + PG_RETURN_NULL(); + } + + /* both rasters are empty */ + if (_isempty[0] && _isempty[1]) { + elog(NOTICE, "The two rasters provided are empty. Returning empty raster"); + + raster = rt_raster_new(0, 0); + if (raster == NULL) { + elog(ERROR, "RASTER_mapAlgebra2Expr: Unable to create empty raster"); + for (k = 0; k < i; k++) { + if (rast[k] != NULL) rt_raster_destroy(rast[k]); + } + PG_RETURN_NULL(); + } + rt_raster_set_scale(raster, 0, 0); + + pgrast = rt_raster_serialize(raster); + rt_raster_destroy(raster); + if (!pgrast) PG_RETURN_NULL(); + + SET_VARSIZE(pgrast, pgrast->size); + PG_RETURN_POINTER(pgrast); + } + + /* replace the empty or NULL raster with one matching the other */ + if ( + (rast[0] == NULL || _isempty[0]) || + (rast[1] == NULL || _isempty[1]) + ) { + if (rast[0] == NULL || _isempty[0]) { + i = 0; + j = 1; + } + else { + i = 1; + j = 0; + } + + _rast[j] = rast[j]; + + /* raster is empty, destroy it */ + if (_rast[i] != NULL) + rt_raster_destroy(_rast[i]); + + _dim[i][0] = rt_raster_get_width(_rast[j]); + _dim[i][1] = rt_raster_get_height(_rast[j]); + _dim[j][0] = rt_raster_get_width(_rast[j]); + _dim[j][1] = rt_raster_get_height(_rast[j]); + + _rast[i] = rt_raster_new( + _dim[j][0], + _dim[j][1] + ); + if (_rast[i] == NULL) { + elog(ERROR, "RASTER_mapAlgebra2Expr: Unable to create nodata raster"); + rt_raster_destroy(_rast[j]); + PG_RETURN_NULL(); + } + rt_raster_set_srid(_rast[i], rt_raster_get_srid(_rast[j])); + + rt_raster_get_geotransform_matrix(_rast[j], gt); + rt_raster_set_geotransform_matrix(_rast[i], gt); + + } + else { + _rast[0] = rast[0]; + _dim[0][0] = rt_raster_get_width(_rast[0]); + _dim[0][1] = rt_raster_get_height(_rast[0]); + + _rast[1] = rast[1]; + _dim[1][0] = rt_raster_get_width(_rast[1]); + _dim[1][1] = rt_raster_get_height(_rast[1]); + } + + /* SRID must match */ + if (rt_raster_get_srid(_rast[0]) != rt_raster_get_srid(_rast[1])) { + elog(NOTICE, "The two rasters provided have different SRIDs. Returning NULL"); + for (k = 0; k < set_count; k++) rt_raster_destroy(_rast[k]); + PG_RETURN_NULL(); + } + + /* same alignment */ + err = rt_raster_same_alignment(_rast[0], _rast[1], &aligned); + if (!err) { + elog(ERROR, "RASTER_mapAlgebra2Expr: Unable to test for alignment on the two rasters"); + for (k = 0; k < set_count; k++) rt_raster_destroy(_rast[k]); + PG_RETURN_NULL(); + } + if (!aligned) { + elog(NOTICE, "The two rasters provided do not have the same alignment. Returning NULL"); + for (k = 0; k < set_count; k++) rt_raster_destroy(_rast[k]); + PG_RETURN_NULL(); + } + + /* pixel type */ + if (!PG_ARGISNULL(5)) { + pixtypename = text_to_cstring(PG_GETARG_TEXT_P(5)); + /* Get the pixel type index */ + pixtype = rt_pixtype_index_from_name(pixtypename); + if (pixtype == PT_END ) { + elog(ERROR, "RASTER_mapAlgebra2Expr: Invalid pixel type: %s", pixtypename); + for (k = 0; k < set_count; k++) rt_raster_destroy(_rast[k]); + PG_RETURN_NULL(); + } + } + + /* extent type */ + if (!PG_ARGISNULL(6)) { + extenttypename = strtoupper(trim(text_to_cstring(PG_GETARG_TEXT_P(6)))); + extenttype = rt_util_extent_type(extenttypename); + } + + /* get raster offsets */ + if (!rt_raster_geopoint_to_cell( + _rast[1], + rt_raster_get_x_offset(_rast[0]), rt_raster_get_y_offset(_rast[0]), + &(_rastoffset[1][0]), &(_rastoffset[1][1]), + NULL + )) { + elog(ERROR, "RASTER_mapAlgebra2Expr: Unable to compute offsets of the second raster relative to the first raster"); + for (k = 0; k < set_count; k++) rt_raster_destroy(_rast[k]); + PG_RETURN_NULL(); + } + _rastoffset[1][0] = -1 * _rastoffset[1][0]; + _rastoffset[1][1] = -1 * _rastoffset[1][1]; + _rastoffset[1][2] = _rastoffset[1][0] + _dim[1][0] - 1; + _rastoffset[1][3] = _rastoffset[1][1] + _dim[1][1] - 1; + + i = 2; + POSTGIS_RT_DEBUGF(3, "extenttype: %d %s", extenttype, extenttypename); + switch (extenttype) { + case ET_FIRST: + i = 0; + _rastoffset[0][0] = 0.; + _rastoffset[0][1] = 0.; + case ET_SECOND: + if (i > 1) { + i = 1; + _rastoffset[0][0] = -1 * _rastoffset[1][0]; + _rastoffset[0][1] = -1 * _rastoffset[1][1]; + _rastoffset[1][0] = 0.; + _rastoffset[1][1] = 0.; + } + + if ( + _isempty[i] && ( + (extenttype == ET_FIRST && i == 0) || + (extenttype == ET_SECOND && i == 1) + ) + ) { + elog(NOTICE, "The %s raster is NULL. Returning NULL", (i != 1 ? "FIRST" : "SECOND")); + for (k = 0; k < set_count; k++) rt_raster_destroy(_rast[k]); + PG_RETURN_NULL(); + } + + dim[0] = _dim[i][0]; + dim[1] = _dim[i][1]; + raster = rt_raster_new( + dim[0], + dim[1] + ); + if (raster == NULL) { + elog(ERROR, "RASTER_mapAlgebra2Expr: Unable to create output raster"); + for (k = 0; k < set_count; k++) rt_raster_destroy(_rast[k]); + PG_RETURN_NULL(); + } + rt_raster_set_srid(raster, rt_raster_get_srid(_rast[i])); + rt_raster_get_geotransform_matrix(_rast[i], gt); + rt_raster_set_geotransform_matrix(raster, gt); + + /* specified band not found */ + if (rt_raster_has_no_band(_rast[i], bandindex[i] - 1)) { + elog(NOTICE, "The %s raster does not have the band at index %d. Returning no band raster of correct extent", + (i != 1 ? "FIRST" : "SECOND"), bandindex[i] + ); + + for (k = 0; k < set_count; k++) rt_raster_destroy(_rast[k]); + + pgrast = rt_raster_serialize(raster); + rt_raster_destroy(raster); + if (!pgrast) PG_RETURN_NULL(); + + SET_VARSIZE(pgrast, pgrast->size); + PG_RETURN_POINTER(pgrast); + } + break; + case ET_UNION: { + double ip[2] = {0}; + double offset[4] = {0}; + + rt_raster_get_geotransform_matrix(_rast[0], gt); + + /* new upper-left */ + ip[0] = rt_raster_get_x_offset(_rast[1]); + ip[1] = rt_raster_get_y_offset(_rast[1]); + if (ip[0] < gt[0]) + gt[0] = ip[0]; + if (ip[1] < gt[3]) + gt[3] = ip[1]; + + /* new width and height */ + offset[0] = 0; + if (_rastoffset[1][0] < 0) + offset[0] = _rastoffset[1][0]; + offset[1] = 0; + if (_rastoffset[1][1] < 0) + offset[1] = _rastoffset[1][1]; + + offset[2] = _dim[0][0] - 1; + if ((int) _rastoffset[1][2] >= _dim[0][0]) + offset[2] = _rastoffset[1][2]; + offset[3] = _dim[0][1] - 1; + if ((int) _rastoffset[1][3] >= _dim[0][1]) + offset[3] = _rastoffset[1][3]; + + dim[0] = offset[2] - offset[0] + 1; + dim[1] = offset[3] - offset[1] + 1; + + raster = rt_raster_new( + dim[0], + dim[1] + ); + if (raster == NULL) { + elog(ERROR, "RASTER_mapAlgebra2Expr: Unable to create output raster"); + for (k = 0; k < set_count; k++) rt_raster_destroy(_rast[k]); + PG_RETURN_NULL(); + } + rt_raster_set_srid(raster, rt_raster_get_srid(_rast[0])); + rt_raster_set_geotransform_matrix(raster, gt); + + if (!rt_raster_geopoint_to_cell( + _rast[0], + gt[0], gt[3], + &(_rastoffset[0][0]), &(_rastoffset[0][1]), + NULL + )) { + elog(ERROR, "RASTER_mapAlgebra2Expr: Unable to get offsets of the FIRST raster relative to the output raster"); + for (k = 0; k < set_count; k++) rt_raster_destroy(_rast[k]); + rt_raster_destroy(raster); + PG_RETURN_NULL(); + } + _rastoffset[0][0] *= -1; + _rastoffset[0][1] *= -1; + + if (!rt_raster_geopoint_to_cell( + _rast[1], + gt[0], gt[3], + &(_rastoffset[1][0]), &(_rastoffset[1][1]), + NULL + )) { + elog(ERROR, "RASTER_mapAlgebra2Expr: Unable to get offsets of the SECOND raster relative to the output raster"); + for (k = 0; k < set_count; k++) rt_raster_destroy(_rast[k]); + rt_raster_destroy(raster); + PG_RETURN_NULL(); + } + _rastoffset[1][0] *= -1; + _rastoffset[1][1] *= -1; + + } break; + case ET_INTERSECTION: { + double offset[4] = {0}; + double ip[2] = {0}; + + /* no intersection */ + if ( + _isempty[0] || + _isempty[1] || + (_rastoffset[1][2] < 0 || _rastoffset[1][0] > (_dim[0][0] - 1)) || + (_rastoffset[1][3] < 0 || _rastoffset[1][1] > (_dim[0][1] - 1)) + ) { + elog(NOTICE, "The two rasters provided have no intersection. Returning no band raster"); + + for (k = 0; k < set_count; k++) rt_raster_destroy(_rast[k]); + + raster = rt_raster_new(0, 0); + rt_raster_set_scale(raster, 0, 0); + rt_raster_set_srid(raster, rt_raster_get_srid(_rast[0])); + + pgrast = rt_raster_serialize(raster); + rt_raster_destroy(raster); + if (!pgrast) PG_RETURN_NULL(); + + SET_VARSIZE(pgrast, pgrast->size); + PG_RETURN_POINTER(pgrast); + } + + if (_rastoffset[1][0] > 0) + offset[0] = _rastoffset[1][0]; + if (_rastoffset[1][1] > 0) + offset[1] = _rastoffset[1][1]; + + offset[2] = _dim[0][0] - 1; + if (_rastoffset[1][2] < _dim[0][0]) + offset[2] = _rastoffset[1][2]; + offset[3] = _dim[0][1] - 1; + if (_rastoffset[1][3] < _dim[0][1]) + offset[3] = _rastoffset[1][3]; + + dim[0] = offset[2] - offset[0] + 1; + dim[1] = offset[3] - offset[1] + 1; + raster = rt_raster_new( + dim[0], + dim[1] + ); + if (raster == NULL) { + elog(ERROR, "RASTER_mapAlgebra2Expr: Unable to create output raster"); + for (k = 0; k < set_count; k++) rt_raster_destroy(_rast[k]); + PG_RETURN_NULL(); + } + rt_raster_set_srid(raster, rt_raster_get_srid(_rast[0])); + + rt_raster_get_geotransform_matrix(_rast[0], gt); + if (!rt_raster_cell_to_geopoint( + _rast[0], + offset[0], offset[1], + &(ip[0]), &(ip[1]), + gt + )) { + elog(ERROR, "RASTER_mapAlgebra2Expr: Unable to get spatial coordinates of upper-left pixel of output raster"); + for (k = 0; k < set_count; k++) rt_raster_destroy(_rast[k]); + rt_raster_destroy(raster); + PG_RETURN_NULL(); + } + + gt[0] = ip[0]; + gt[3] = ip[1]; + rt_raster_set_geotransform_matrix(raster, gt); + + if (!rt_raster_geopoint_to_cell( + _rast[0], + gt[0], gt[3], + &(_rastoffset[0][0]), &(_rastoffset[0][1]), + NULL + )) { + elog(ERROR, "RASTER_mapAlgebra2Expr: Unable to get pixel coordinates to compute the offsets of the FIRST raster relative to the output raster"); + for (k = 0; k < set_count; k++) rt_raster_destroy(_rast[k]); + rt_raster_destroy(raster); + PG_RETURN_NULL(); + } + _rastoffset[0][0] *= -1; + _rastoffset[0][1] *= -1; + + if (!rt_raster_geopoint_to_cell( + _rast[1], + gt[0], gt[3], + &(_rastoffset[1][0]), &(_rastoffset[1][1]), + NULL + )) { + elog(ERROR, "RASTER_mapAlgebra2Expr: Unable to get pixel coordinates to compute the offsets of the SECOND raster relative to the output raster"); + for (k = 0; k < set_count; k++) rt_raster_destroy(_rast[k]); + rt_raster_destroy(raster); + PG_RETURN_NULL(); + } + _rastoffset[1][0] *= -1; + _rastoffset[1][1] *= -1; + + } break; + } + + /* both rasters do not have specified bands */ + if ( + (!_isempty[0] && rt_raster_has_no_band(_rast[0], bandindex[0] - 1)) && + (!_isempty[1] && rt_raster_has_no_band(_rast[1], bandindex[1] - 1)) + ) { + elog(NOTICE, "The two rasters provided do not have the respectively specified band indices. Returning no band raster of correct extent"); + + for (k = 0; k < set_count; k++) rt_raster_destroy(_rast[k]); + + pgrast = rt_raster_serialize(raster); + rt_raster_destroy(raster); + if (!pgrast) PG_RETURN_NULL(); + + SET_VARSIZE(pgrast, pgrast->size); + PG_RETURN_POINTER(pgrast); + } + + /* get bands */ + for (i = 0; i < set_count; i++) { + if (_isempty[i] || rt_raster_has_no_band(_rast[i], bandindex[i] - 1)) { + _hasnodata[i] = 1; + _nodataval[i] = 0; + + continue; + } + + _band[i] = rt_raster_get_band(_rast[i], bandindex[i] - 1); + if (_band[i] == NULL) { + elog(ERROR, "RASTER_mapAlgebra2Expr: Unable to get band %d of the %s raster", + bandindex[i], + (i < 1 ? "FIRST" : "SECOND") + ); + for (k = 0; k < set_count; k++) rt_raster_destroy(_rast[k]); + rt_raster_destroy(raster); + PG_RETURN_NULL(); + } + + _hasnodata[i] = rt_band_get_hasnodata_flag(_band[i]); + if (_hasnodata[i]) + _nodataval[i] = rt_band_get_nodata(_band[i]); + } + + /* pixtype is PT_END, get pixtype based upon extent */ + if (pixtype == PT_END) { + if ((extenttype == ET_SECOND && !_isempty[1]) || _isempty[0]) + pixtype = rt_band_get_pixtype(_band[1]); + else + pixtype = rt_band_get_pixtype(_band[0]); + } + + /* nodata value for new band */ + if (extenttype == ET_SECOND && !_isempty[1] && _hasnodata[1]) { + nodataval = _nodataval[1]; + } + else if (!_isempty[0] && _hasnodata[0]) { + nodataval = _nodataval[0]; + } + else if (!_isempty[1] && _hasnodata[1]) { + nodataval = _nodataval[1]; + } + else { + elog(NOTICE, "Neither raster provided has a NODATA value for the specified band indices. NODATA value set to minimum possible for %s", rt_pixtype_name(pixtype)); + nodataval = rt_pixtype_get_min_value(pixtype); + } + + /* add band to output raster */ + if (rt_raster_generate_new_band( + raster, + pixtype, nodataval, + 1, nodataval, + 0 + ) < 0) { + elog(ERROR, "RASTER_mapAlgebra2Expr: Unable to add new band to output raster"); + for (k = 0; k < set_count; k++) rt_raster_destroy(_rast[k]); + rt_raster_destroy(raster); + PG_RETURN_NULL(); + } + + /* get output band */ + band = rt_raster_get_band(raster, 0); + if (band == NULL) { + elog(ERROR, "RASTER_mapAlgebra2Expr: Unable to get newly added band of output raster"); + for (k = 0; k < set_count; k++) rt_raster_destroy(_rast[k]); + rt_raster_destroy(raster); + PG_RETURN_NULL(); + } + + POSTGIS_RT_DEBUGF(4, "offsets = (%d, %d, %d, %d)", + (int) _rastoffset[0][0], + (int) _rastoffset[0][1], + (int) _rastoffset[1][0], + (int) _rastoffset[1][1] + ); + + POSTGIS_RT_DEBUGF(4, "metadata = (%f, %f, %d, %d, %f, %f, %f, %f, %d)", + rt_raster_get_x_offset(raster), + rt_raster_get_y_offset(raster), + rt_raster_get_width(raster), + rt_raster_get_height(raster), + rt_raster_get_x_scale(raster), + rt_raster_get_y_scale(raster), + rt_raster_get_x_skew(raster), + rt_raster_get_y_skew(raster), + rt_raster_get_srid(raster) + ); + + POSTGIS_RT_DEBUGF(4, "bandmetadata = (%s, %d, %f)", + rt_pixtype_name(rt_band_get_pixtype(band)), + rt_band_get_hasnodata_flag(band), + rt_band_get_nodata(band) + ); + + /* connect SPI */ + if (SPI_connect() != SPI_OK_CONNECT) { + elog(ERROR, "RASTER_mapAlgebra2Expr: Unable to connect to the SPI manager"); + for (k = 0; k < set_count; k++) rt_raster_destroy(_rast[k]); + rt_raster_destroy(raster); + PG_RETURN_NULL(); + } + + /* + process expressions + + exprpos elements are: + 4 - expression => spiplan[0] + 7 - nodata1expr => spiplan[1] + 8 - nodata2expr => spiplan[2] + */ + for (i = 0; i < spicount; i++) { + if (!PG_ARGISNULL(exprpos[i])) { + expr = strtoupper(text_to_cstring(PG_GETARG_TEXT_P(exprpos[i]))); + POSTGIS_RT_DEBUGF(3, "raw expr #%d: %s", i, expr); + len = strlen(expr); + expr = replace(expr, "RAST1", "$1", NULL); + if (strlen(expr) != len) { + argcount[i]++; + argexists[i][0] = 1; + } + len = strlen(expr); + expr = replace(expr, "RAST2", (argexists[i][0] ? "$2" : "$1"), NULL); + if (strlen(expr) != len) { + argcount[i]++; + argexists[i][1] = 1; + } + + len = strlen("SELECT (") + strlen(expr) + strlen(")::double precision"); + sql = (char *) palloc(len + 1); + if (sql == NULL) { + elog(ERROR, "RASTER_mapAlgebra2Expr: Unable to allocate memory for expression parameter %d", exprpos[i]); + for (k = 0; k < set_count; k++) rt_raster_destroy(_rast[k]); + rt_raster_destroy(raster); + for (k = 0; k < spicount; k++) SPI_freeplan(spiplan[k]); + SPI_finish(); + PG_RETURN_NULL(); + } + + strncpy(sql, "SELECT (", strlen("SELECT (")); + strncpy(sql + strlen("SELECT ("), expr, strlen(expr)); + strncpy(sql + strlen("SELECT (") + strlen(expr), ")::double precision", strlen(")::double precision")); + sql[len] = '\0'; + + POSTGIS_RT_DEBUGF(3, "sql #%d: %s", i, sql); + + /* create prepared plan */ + argtype = (Oid *) palloc(argcount[i] * sizeof(Oid)); + if (argtype == NULL) { + elog(ERROR, "RASTER_mapAlgebra2Expr: Unable to allocate memory for prepared plan argtypes of expression parameter %d", exprpos[i]); + + for (k = 0; k < set_count; k++) rt_raster_destroy(_rast[k]); + rt_raster_destroy(raster); + + for (k = 0; k < spicount; k++) SPI_freeplan(spiplan[k]); + SPI_finish(); + + pfree(sql); + + PG_RETURN_NULL(); + } + for (j = 0; j < argcount[i]; j++) argtype[j] = FLOAT8OID; + + spiplan[i] = SPI_prepare(sql, argcount[i], argtype); + if (spiplan[i] == NULL) { + elog(ERROR, "RASTER_mapAlgebra2Expr: Unable to create prepared plan of expression parameter %d", exprpos[i]); + + for (k = 0; k < set_count; k++) rt_raster_destroy(_rast[k]); + rt_raster_destroy(raster); + + for (k = 0; k < spicount; k++) SPI_freeplan(spiplan[k]); + SPI_finish(); + + pfree(sql); + pfree(argtype); + + PG_RETURN_NULL(); + } + + pfree(sql); + pfree(argtype); + } + else + spiempty++; + } + + /* nodatanodataval */ + if (!PG_ARGISNULL(9)) { + hasnodatanodataval = 1; + nodatanodataval = PG_GETARG_FLOAT8(9); + } + else + hasnodatanodataval = 0; + + /* loop over pixels */ + /* if spiempty == 4, no need to do this */ + if (spiempty != 4) { + for (x = 0; x < dim[0]; x++) { + for (y = 0; y < dim[1]; y++) { + + /* get pixel from each raster */ + for (i = 0; i < set_count; i++) { + _haspixel[i] = 0; + _pixel[i] = 0; + + /* row/column */ + _x = x - (int) _rastoffset[i][0]; + _y = y - (int) _rastoffset[i][1]; + + /* get pixel value */ + if ( + !_isempty[i] && + (_x >= 0 && _x < _dim[i][0]) && + (_y >= 0 && _y < _dim[i][1]) + ) { + err = rt_band_get_pixel(_band[i], _x, _y, &(_pixel[i])); + if (err < 0) { + elog(ERROR, "RASTER_mapAlgebra2Expr: Unable to get pixel of %s raster", (i < 1 ? "FIRST" : "SECOND")); + + for (k = 0; k < set_count; k++) rt_raster_destroy(_rast[k]); + rt_raster_destroy(raster); + + for (k = 0; k < spicount; k++) SPI_freeplan(spiplan[k]); + SPI_finish(); + + PG_RETURN_NULL(); + } + if (!_hasnodata[i] || FLT_NEQ(_nodataval[i], _pixel[i])) + _haspixel[i] = 1; + } + + POSTGIS_RT_DEBUGF(5, "r%d(%d, %d) = %d, %f", + i, + _x, + _y, + _haspixel[i], + _pixel[i] + ); + } + + haspixel = 0; + + /* which prepared plan to use? */ + /* !pixel0 && !pixel1 */ + /* use nodatanodataval */ + if (!_haspixel[0] && !_haspixel[1]) + i = 3; + /* pixel0 && !pixel1 */ + /* run spiplan[2] (nodata2expr) */ + else if (_haspixel[0] && !_haspixel[1]) + i = 2; + /* !pixel0 && pixel1 */ + /* run spiplan[1] (nodata1expr) */ + else if (!_haspixel[0] && _haspixel[1]) + i = 1; + /* pixel0 && pixel1 */ + /* run spiplan[0] (expression) */ + else + i = 0; + + /* process values */ + if (i == 3) { + if (hasnodatanodataval) { + haspixel = 1; + pixel = nodatanodataval; + } + } + else if (spiplan[i] != NULL) { + POSTGIS_RT_DEBUGF(5, "Using prepared plan: %d", i); + + /* expression has argument(s) */ + if (argcount[i]) { + values = (Datum *) palloc(sizeof(Datum) * argcount[i]); + if (values == NULL) { + elog(ERROR, "RASTER_mapAlgebra2Expr: Unable to allocate memory for value parameters of prepared statement %d", i); + + for (k = 0; k < set_count; k++) rt_raster_destroy(_rast[k]); + rt_raster_destroy(raster); + + for (k = 0; k < spicount; k++) SPI_freeplan(spiplan[k]); + SPI_finish(); + + PG_RETURN_NULL(); + } + + nulls = (bool *) palloc(sizeof(bool) * argcount[i]); + if (nulls == NULL) { + elog(ERROR, "RASTER_mapAlgebra2Expr: Unable to allocate memory for NULL parameters of prepared statement %d", i); + + pfree(values); + + for (k = 0; k < set_count; k++) rt_raster_destroy(_rast[k]); + rt_raster_destroy(raster); + + for (k = 0; k < spicount; k++) SPI_freeplan(spiplan[k]); + SPI_finish(); + + PG_RETURN_NULL(); + } + memset(nulls, FALSE, argcount[i]); + + /* two arguments */ + if (argcount[i] > 1) { + for (j = 0; j < argcount[i]; j++) { + if (_isempty[j] || !_haspixel[j]) + nulls[j] = TRUE; + else + values[j] = Float8GetDatum(_pixel[j]); + } + } + /* only one argument */ + else { + if (argexists[i][0]) + j = 0; + else + j = 1; + + if (_isempty[j] || !_haspixel[j]) { + POSTGIS_RT_DEBUG(5, "using null"); + nulls[0] = TRUE; + } + else { + POSTGIS_RT_DEBUGF(5, "using value %f", _pixel[j]); + values[0] = Float8GetDatum(_pixel[j]); + } + } + } + + /* run prepared plan */ + err = SPI_execute_plan(spiplan[i], values, nulls, TRUE, 1); + if (values != NULL) pfree(values); + if (nulls != NULL) pfree(nulls); + if (err != SPI_OK_SELECT || SPI_tuptable == NULL || SPI_processed != 1) { + elog(ERROR, "RASTER_mapAlgebra2Expr: Unexpected error when running prepared statement %d", i); + + for (k = 0; k < set_count; k++) rt_raster_destroy(_rast[k]); + rt_raster_destroy(raster); + + for (k = 0; k < spicount; k++) SPI_freeplan(spiplan[k]); + SPI_finish(); + + PG_RETURN_NULL(); + } + + /* get output of prepared plan */ + tupdesc = SPI_tuptable->tupdesc; + tuptable = SPI_tuptable; + tuple = tuptable->vals[0]; + + pixeldatum = SPI_getbinval(tuple, tupdesc, 1, &isnull); + if (SPI_result == SPI_ERROR_NOATTRIBUTE) { + elog(ERROR, "RASTER_mapAlgebra2Expr: Unable to get result of prepared statement %d", i); + + for (k = 0; k < set_count; k++) rt_raster_destroy(_rast[k]); + rt_raster_destroy(raster); + + if (SPI_tuptable) SPI_freetuptable(tuptable); + for (k = 0; k < spicount; k++) SPI_freeplan(spiplan[k]); + SPI_finish(); + + PG_RETURN_NULL(); + } + + if (!isnull) { + haspixel = 1; + pixel = DatumGetFloat8(pixeldatum); + } + + if (SPI_tuptable) SPI_freetuptable(tuptable); + } + + /* burn pixel if haspixel != 0 */ + if (haspixel) { + if (rt_band_set_pixel(band, x, y, pixel) < 0) { + elog(ERROR, "RASTER_mapAlgebra2Expr: Unable to set pixel value of output raster"); + + for (k = 0; k < set_count; k++) rt_raster_destroy(_rast[k]); + rt_raster_destroy(raster); + + for (k = 0; k < spicount; k++) SPI_freeplan(spiplan[k]); + SPI_finish(); + + PG_RETURN_NULL(); + } + } + + POSTGIS_RT_DEBUGF(5, "(x, y, val) = (%d, %d, %f)", x, y, haspixel ? pixel : nodataval); + } + } + } + + /* CLEANUP */ + for (i = 0; i < spicount; i++) { + if (spiplan[i] != NULL) SPI_freeplan(spiplan[i]); + } + SPI_finish(); + + for (k = 0; k < set_count; k++) { + if (_rast[k] != NULL) rt_raster_destroy(_rast[k]); + } + + pgrast = rt_raster_serialize(raster); + rt_raster_destroy(raster); + if (!pgrast) PG_RETURN_NULL(); + + SET_VARSIZE(pgrast, pgrast->size); + PG_RETURN_POINTER(pgrast); +} + /* ---------------------------------------------------------------- */ /* Memory allocation / error reporting hooks */ /* ---------------------------------------------------------------- */ diff --git a/raster/rt_pg/rt_pg.h b/raster/rt_pg/rt_pg.h index b35044895..608647b13 100644 --- a/raster/rt_pg/rt_pg.h +++ b/raster/rt_pg/rt_pg.h @@ -97,37 +97,6 @@ typedef struct rt_pgband_t { /* Header of PostgreSQL-stored RASTER value, * and binary representation of it */ -typedef struct rt_pgraster_header_t { - - /*---[ 8 byte boundary ]---{ */ - uint32_t size; /* required by postgresql: 4 bytes */ - uint16_t version; /* format version (this is version 0): 2 bytes */ - uint16_t numBands; /* Number of bands: 2 bytes */ - - /* }---[ 8 byte boundary ]---{ */ - double scaleX; /* pixel width: 8 bytes */ - - /* }---[ 8 byte boundary ]---{ */ - double scaleY; /* pixel height: 8 bytes */ - - /* }---[ 8 byte boundary ]---{ */ - double ipX; /* insertion point X: 8 bytes */ - - /* }---[ 8 byte boundary ]---{ */ - double ipY; /* insertion point Y: 8 bytes */ - - /* }---[ 8 byte boundary ]---{ */ - double skewX; /* skew about the X axis: 8 bytes */ - - /* }---[ 8 byte boundary ]---{ */ - double skewY; /* skew about the Y axis: 8 bytes */ - - /* }---[ 8 byte boundary ]--- */ - int32_t srid; /* Spatial reference id: 4 bytes */ - uint16_t width; /* pixel columns: 2 bytes */ - uint16_t height; /* pixel rows: 2 bytes */ -} rt_pgraster_header; - -typedef rt_pgraster_header rt_pgraster; +typedef struct rt_raster_serialized_t rt_pgraster; #endif /* RT_PG_H_INCLUDED */ diff --git a/raster/rt_pg/rtpostgis.sql.in.c b/raster/rt_pg/rtpostgis.sql.in.c index 3e88ed68d..abc6eae48 100644 --- a/raster/rt_pg/rtpostgis.sql.in.c +++ b/raster/rt_pg/rtpostgis.sql.in.c @@ -1660,6 +1660,31 @@ CREATE OR REPLACE FUNCTION st_mapalgebrafct(rast raster, userfunction regprocedu AS $$ SELECT st_mapalgebrafct($1, 1, NULL, $2, NULL) $$ LANGUAGE SQL; +----------------------------------------------------------------------- +-- Two Raster ST_MapAlgebra +----------------------------------------------------------------------- +CREATE OR REPLACE FUNCTION st_mapalgebra2expr( + rast1 raster, band1 integer, + rast2 raster, band2 integer, + expression text, + pixeltype text DEFAULT NULL, extenttype text DEFAULT 'INTERSECTION', + nodata1expr text DEFAULT NULL, nodata2expr text DEFAULT NULL, nodatanodataval double precision DEFAULT NULL +) + RETURNS raster + AS 'MODULE_PATHNAME', 'RASTER_mapAlgebra2Expr' + LANGUAGE 'C' IMMUTABLE; + +CREATE OR REPLACE FUNCTION st_mapalgebra2expr( + rast1 raster, + rast2 raster, + expression text, + pixeltype text DEFAULT NULL, extenttype text DEFAULT 'INTERSECTION', + nodata1expr text DEFAULT NULL, nodata2expr text DEFAULT NULL, nodatanodataval double precision DEFAULT NULL +) + RETURNS raster + AS $$ SELECT st_mapalgebra2expr($1, 1, $2, 1, $3, $4, $5, $6, $7, $8) $$ + LANGUAGE 'SQL' IMMUTABLE; + ----------------------------------------------------------------------- -- Get information about the raster ----------------------------------------------------------------------- @@ -2676,7 +2701,7 @@ CREATE OR REPLACE FUNCTION _st_intersects(geom geometry, rast raster, nband inte BEGIN convexhull := ST_ConvexHull(rast); IF nband IS NOT NULL THEN - SELECT bmd.hasnodata INTO hasnodata FROM ST_BandMetaData(rast, nband) bmd; + SELECT bmd.hasnodata INTO hasnodata FROM ST_BandMetaData(rast, nband) AS bmd; END IF; IF ST_Intersects(geom, convexhull) IS NOT TRUE THEN diff --git a/raster/scripts/plpgsql/st_mapalgebra.sql b/raster/scripts/plpgsql/st_mapalgebra.sql index 3756e5a7b..5c6086903 100644 --- a/raster/scripts/plpgsql/st_mapalgebra.sql +++ b/raster/scripts/plpgsql/st_mapalgebra.sql @@ -571,7 +571,7 @@ CREATE OR REPLACE FUNCTION ST_MapAlgebra(rast1 raster, newscalex := rast1scalex; newscaley := rast1scaley; newskewx := rast1skewx; - newskewy := rast1skewx; + newskewy := rast1skewy; --rastoffset is the offset of a raster in relatively to the first raster rast2offsetx1 := -st_world2rastercoordx(rast2, rast1ulx, rast1uly) + 1; diff --git a/raster/scripts/plpgsql/st_mapalgebra_optimized.sql b/raster/scripts/plpgsql/st_mapalgebra_optimized.sql index 44810a77c..9b72a188a 100644 --- a/raster/scripts/plpgsql/st_mapalgebra_optimized.sql +++ b/raster/scripts/plpgsql/st_mapalgebra_optimized.sql @@ -276,7 +276,7 @@ RAISE NOTICE 'ST_MapAlgebra2 000'; newscalex := rast1scalex; newscaley := rast1scaley; newskewx := rast1skewx; - newskewy := rast1skewx; + newskewy := rast1skewy; --r1x & r2x are the offset of each rasters relatively to global extent r1x := 0; diff --git a/raster/test/regress/Makefile.in b/raster/test/regress/Makefile.in index 1b6655954..09be1fb54 100644 --- a/raster/test/regress/Makefile.in +++ b/raster/test/regress/Makefile.in @@ -87,7 +87,8 @@ TEST_UTILITY = \ rt_utility.sql \ create_rt_mapalgebra_test.sql \ rt_mapalgebraexpr.sql \ - rt_mapalgebrafct.sql \ + rt_mapalgebrafct.sql \ + rt_mapalgebra2expr.sql \ rt_reclass.sql \ rt_resample.sql \ rt_asraster.sql \ diff --git a/raster/test/regress/rt_mapalgebra2expr.sql b/raster/test/regress/rt_mapalgebra2expr.sql new file mode 100644 index 000000000..64128a025 --- /dev/null +++ b/raster/test/regress/rt_mapalgebra2expr.sql @@ -0,0 +1,253 @@ +DROP TABLE IF EXISTS raster_mapalgebra; +CREATE TABLE raster_mapalgebra ( + rid integer, + rast raster +); +DROP TABLE IF EXISTS raster_mapalgebra_out; +CREATE TABLE raster_mapalgebra_out ( + rid1 integer, + rid2 integer, + extent varchar, + rast raster +); +CREATE OR REPLACE FUNCTION make_test_raster( + rid integer, + width integer DEFAULT 2, + height integer DEFAULT 2, + ul_x double precision DEFAULT 0, + ul_y double precision DEFAULT 0, + skew_x double precision DEFAULT 0, + skew_y double precision DEFAULT 0, + initvalue double precision DEFAULT 1, + nodataval double precision DEFAULT 0 +) + RETURNS void + AS $$ + DECLARE + x int; + y int; + rast raster; + BEGIN + rast := ST_MakeEmptyRaster(width, height, ul_x, ul_y, 1, 1, skew_x, skew_y, 0); + rast := ST_AddBand(rast, 1, '8BUI', initvalue, nodataval); + + + INSERT INTO raster_mapalgebra VALUES (rid, rast); + + RETURN; + END; + $$ LANGUAGE 'plpgsql'; +-- no skew +SELECT make_test_raster(0, 4, 4, -2, -2); +SELECT make_test_raster(1, 2, 2, 0, 0, 0, 0, 2); +SELECT make_test_raster(2, 2, 2, 1, -1, 0, 0, 3); +SELECT make_test_raster(3, 2, 2, 1, 1, 0, 0, 4); +SELECT make_test_raster(4, 2, 2, 2, 2, 0, 0, 5); + +-- skew +SELECT make_test_raster(10, 4, 4, -2, -2, 1, -1); +SELECT make_test_raster(11, 2, 2, 0, 0, 1, -1, 2); +SELECT make_test_raster(12, 2, 2, 1, -1, 1, -1, 3); +SELECT make_test_raster(13, 2, 2, 1, 1, 1, -1, 4); +SELECT make_test_raster(14, 2, 2, 2, 2, 1, -1, 5); + +DROP FUNCTION make_test_raster(integer, integer, integer, double precision, double precision, double precision, double precision, double precision, double precision); + +-- INTERSECTION +INSERT INTO raster_mapalgebra_out + (SELECT r1.rid, r2.rid, 'INTERSECTION', st_mapalgebra2expr( + r1.rast, r2.rast, 'rast1', '32BF', 'INTERSECTION' + ) + FROM raster_mapalgebra r1 + JOIN raster_mapalgebra r2 + ON r1.rid != r2.rid + WHERE r1.rid = 0 + AND r2.rid BETWEEN 1 AND 9 + ) UNION ALL ( + SELECT r1.rid, r2.rid, 'INTERSECTION', st_mapalgebra2expr( + r1.rast, r2.rast, 'rast1', '32BF', 'INTERSECTION' + ) + FROM raster_mapalgebra r1 + JOIN raster_mapalgebra r2 + ON r1.rid != r2.rid + WHERE r1.rid = 10 + AND r2.rid BETWEEN 11 AND 19) +; + +INSERT INTO raster_mapalgebra_out + SELECT NULL AS rid, rid, 'INTERSECTION', st_mapalgebra2expr( + NULL::raster, rast, 'rast1', '32BF', 'INTERSECTION' + ) + FROM raster_mapalgebra +; + +INSERT INTO raster_mapalgebra_out + SELECT rid, NULL AS rid, 'INTERSECTION', st_mapalgebra2expr( + rast, NULL::raster, 'rast1', '32BF', 'INTERSECTION' + ) + FROM raster_mapalgebra +; + +INSERT INTO raster_mapalgebra_out + SELECT NULL AS rid, NULL AS rid, 'INTERSECTION', st_mapalgebra2expr( + NULL::raster, NULL::raster, 'rast1', '32BF', 'INTERSECTION' + ) +; + +-- UNION +INSERT INTO raster_mapalgebra_out + (SELECT r1.rid, r2.rid, 'UNION', st_mapalgebra2expr( + r1.rast, r2.rast, '((rast1 + rast2)/2.)::numeric', '32BF', 'UNION', 'rast2', 'rast1', NULL + ) + FROM raster_mapalgebra r1 + JOIN raster_mapalgebra r2 + ON r1.rid != r2.rid + WHERE r1.rid = 0 + AND r2.rid BETWEEN 1 AND 9 + ) UNION ALL ( + SELECT r1.rid, r2.rid, 'UNION', st_mapalgebra2expr( + r1.rast, r2.rast, '((rast1 + rast2)/2.)::numeric', '32BF', 'UNION', 'rast2', 'rast1', NULL + ) + FROM raster_mapalgebra r1 + JOIN raster_mapalgebra r2 + ON r1.rid != r2.rid + WHERE r1.rid = 10 + AND r2.rid BETWEEN 11 AND 19) +; + +INSERT INTO raster_mapalgebra_out + SELECT NULL AS rid, rid, 'UNION', st_mapalgebra2expr( + NULL::raster, rast, '((rast1 + rast2)/2.)::numeric', '32BF', 'UNION', 'rast2', 'rast1', NULL + ) + FROM raster_mapalgebra +; + +INSERT INTO raster_mapalgebra_out + SELECT rid, NULL AS rid, 'UNION', st_mapalgebra2expr( + rast, NULL::raster, '((rast1 + rast2)/2.)::numeric', '32BF', 'UNION', 'rast2', 'rast1', NULL + ) + FROM raster_mapalgebra +; + +INSERT INTO raster_mapalgebra_out + SELECT NULL AS rid, NULL AS rid, 'UNION', st_mapalgebra2expr( + NULL::raster, NULL::raster, '((rast1 + rast2)/2.)::numeric', '32BF', 'UNION', 'rast2', 'rast1', NULL + ) +; + +-- FIRST +INSERT INTO raster_mapalgebra_out + (SELECT r1.rid, r2.rid, 'FIRST', st_mapalgebra2expr( + r1.rast, r2.rast, 'CASE WHEN rast2 IS NOT NULL THEN NULL ELSE rast1 END', '32BF', 'FIRST', NULL, 'rast1', NULL + ) + FROM raster_mapalgebra r1 + JOIN raster_mapalgebra r2 + ON r1.rid != r2.rid + WHERE r1.rid = 0 + AND r2.rid BETWEEN 1 AND 9 + ) UNION ALL ( + SELECT r1.rid, r2.rid, 'FIRST', st_mapalgebra2expr( + r1.rast, r2.rast, '((rast1 + rast2)/2.)::numeric', '32BF', 'UNION', 'rast2', 'rast1', NULL + ) + FROM raster_mapalgebra r1 + JOIN raster_mapalgebra r2 + ON r1.rid != r2.rid + WHERE r1.rid = 10 + AND r2.rid BETWEEN 11 AND 19) +; + +INSERT INTO raster_mapalgebra_out + SELECT NULL AS rid, rid, 'FIRST', st_mapalgebra2expr( + NULL::raster, rast, 'CASE WHEN rast1 IS NOT NULL THEN NULL ELSE rast2 END', '32BF', 'SECOND', 'rast2', NULL, NULL + ) + FROM raster_mapalgebra +; + +INSERT INTO raster_mapalgebra_out + SELECT rid, NULL AS rid, 'FIRST', st_mapalgebra2expr( + rast, NULL::raster, 'CASE WHEN rast2 IS NOT NULL THEN NULL ELSE rast1 END', '32BF', 'FIRST', NULL, 'rast1', NULL + ) + FROM raster_mapalgebra +; + +INSERT INTO raster_mapalgebra_out + SELECT NULL AS rid, NULL AS rid, 'FIRST', st_mapalgebra2expr( + NULL::raster, NULL::raster, 'CASE WHEN rast2 IS NOT NULL THEN NULL ELSE rast1 END', '32BF', 'FIRST', NULL, 'rast1', NULL + ) +; + +-- SECOND +INSERT INTO raster_mapalgebra_out + (SELECT r1.rid, r2.rid, 'SECOND', st_mapalgebra2expr( + r1.rast, r2.rast, 'CASE WHEN rast1 IS NOT NULL THEN NULL ELSE rast2 END', '32BF', 'SECOND', 'rast2', NULL, NULL + ) + FROM raster_mapalgebra r1 + JOIN raster_mapalgebra r2 + ON r1.rid != r2.rid + WHERE r1.rid = 0 + AND r2.rid BETWEEN 1 AND 9 + ) UNION ALL ( + SELECT r1.rid, r2.rid, 'SECOND', st_mapalgebra2expr( + r1.rast, r2.rast, 'CASE WHEN rast1 IS NOT NULL THEN NULL ELSE rast2 END', '32BF', 'SECOND', 'rast2', NULL, NULL + ) + FROM raster_mapalgebra r1 + JOIN raster_mapalgebra r2 + ON r1.rid != r2.rid + WHERE r1.rid = 10 + AND r2.rid BETWEEN 11 AND 19) +; + +INSERT INTO raster_mapalgebra_out + SELECT NULL AS rid, rid, 'SECOND', st_mapalgebra2expr( + NULL::raster, rast, 'CASE WHEN rast1 IS NOT NULL THEN NULL ELSE rast2 END', '32BF', 'SECOND', 'rast2', NULL, NULL + ) + FROM raster_mapalgebra +; + +INSERT INTO raster_mapalgebra_out + SELECT rid, NULL AS rid, 'SECOND', st_mapalgebra2expr( + rast, NULL::raster, 'CASE WHEN rast1 IS NOT NULL THEN NULL ELSE rast2 END', '32BF', 'SECOND', 'rast2', NULL, NULL + ) + FROM raster_mapalgebra +; + +INSERT INTO raster_mapalgebra_out + SELECT NULL AS rid, NULL AS rid, 'SECOND', st_mapalgebra2expr( + NULL::raster, NULL::raster, 'CASE WHEN rast1 IS NOT NULL THEN NULL ELSE rast2 END', '32BF', 'SECOND', 'rast2', NULL, NULL + ) +; + +-- output +SELECT + rid1, + rid2, + extent, + round(upperleftx::numeric, 3) AS upperleftx, + round(upperlefty::numeric, 3) AS upperlefty, + width, + height, + round(scalex::numeric, 3) AS scalex, + round(scaley::numeric, 3) AS scaley, + round(skewx::numeric, 3) AS skewx, + round(skewy::numeric, 3) AS skewy, + srid, + numbands, + pixeltype, + hasnodata, + round(nodatavalue::numeric, 3) AS nodatavalue, + round(firstvalue::numeric, 3) AS firstvalue, + round(lastvalue::numeric, 3) AS lastvalue +FROM ( + SELECT + rid1, + rid2, + extent, + (ST_Metadata(rast)).*, + (ST_BandMetadata(rast, 1)).*, + ST_Value(rast, 1, 1, 1) AS firstvalue, + ST_Value(rast, 1, ST_Width(rast), ST_Height(rast)) AS lastvalue + FROM raster_mapalgebra_out +) AS r; + +DROP TABLE IF EXISTS raster_mapalgebra; +DROP TABLE IF EXISTS raster_mapalgebra_out; diff --git a/raster/test/regress/rt_mapalgebra2expr_expected b/raster/test/regress/rt_mapalgebra2expr_expected new file mode 100644 index 000000000..dc96f5ab8 --- /dev/null +++ b/raster/test/regress/rt_mapalgebra2expr_expected @@ -0,0 +1,264 @@ +NOTICE: table "raster_mapalgebra" does not exist, skipping +NOTICE: table "raster_mapalgebra_out" does not exist, skipping +NOTICE: The two rasters provided have no intersection. Returning no band raster +NOTICE: The two rasters provided have no intersection. Returning no band raster +NOTICE: The two rasters provided have no intersection. Returning no band raster +NOTICE: The two rasters provided have no intersection. Returning no band raster +NOTICE: The two rasters provided have no intersection. Returning no band raster +NOTICE: The two rasters provided have no intersection. Returning no band raster +NOTICE: The two rasters provided have no intersection. Returning no band raster +NOTICE: The two rasters provided have no intersection. Returning no band raster +NOTICE: The two rasters provided have no intersection. Returning no band raster +NOTICE: The two rasters provided have no intersection. Returning no band raster +NOTICE: The two rasters provided have no intersection. Returning no band raster +NOTICE: The two rasters provided have no intersection. Returning no band raster +NOTICE: The two rasters provided have no intersection. Returning no band raster +NOTICE: The two rasters provided have no intersection. Returning no band raster +NOTICE: The two rasters provided have no intersection. Returning no band raster +NOTICE: The two rasters provided have no intersection. Returning no band raster +NOTICE: The two rasters provided have no intersection. Returning no band raster +NOTICE: The two rasters provided have no intersection. Returning no band raster +NOTICE: The two rasters provided have no intersection. Returning no band raster +NOTICE: The two rasters provided have no intersection. Returning no band raster +NOTICE: The two rasters provided have no intersection. Returning no band raster +NOTICE: The two rasters provided have no intersection. Returning no band raster +NOTICE: The two rasters provided are NULL. Returning NULL +NOTICE: The two rasters provided are NULL. Returning NULL +NOTICE: The two rasters provided are NULL. Returning NULL +NOTICE: The SECOND raster is NULL. Returning NULL +NOTICE: The SECOND raster is NULL. Returning NULL +NOTICE: The SECOND raster is NULL. Returning NULL +NOTICE: The SECOND raster is NULL. Returning NULL +NOTICE: The SECOND raster is NULL. Returning NULL +NOTICE: The SECOND raster is NULL. Returning NULL +NOTICE: The SECOND raster is NULL. Returning NULL +NOTICE: The SECOND raster is NULL. Returning NULL +NOTICE: The SECOND raster is NULL. Returning NULL +NOTICE: The SECOND raster is NULL. Returning NULL +NOTICE: The two rasters provided are NULL. Returning NULL +NOTICE: Raster provided has no bands +NOTICE: Raster provided has no bands +NOTICE: Raster provided has no bands +NOTICE: Could not find raster band of index 1 when getting pixel value. Returning NULL +NOTICE: Could not find raster band of index 1 when getting pixel value. Returning NULL +NOTICE: Raster provided has no bands +NOTICE: Raster provided has no bands +NOTICE: Raster provided has no bands +NOTICE: Could not find raster band of index 1 when getting pixel value. Returning NULL +NOTICE: Could not find raster band of index 1 when getting pixel value. Returning NULL +NOTICE: Raster provided has no bands +NOTICE: Raster provided has no bands +NOTICE: Raster provided has no bands +NOTICE: Could not find raster band of index 1 when getting pixel value. Returning NULL +NOTICE: Could not find raster band of index 1 when getting pixel value. Returning NULL +NOTICE: Raster provided has no bands +NOTICE: Raster provided has no bands +NOTICE: Raster provided has no bands +NOTICE: Could not find raster band of index 1 when getting pixel value. Returning NULL +NOTICE: Could not find raster band of index 1 when getting pixel value. Returning NULL +NOTICE: Raster provided has no bands +NOTICE: Raster provided has no bands +NOTICE: Raster provided has no bands +NOTICE: Could not find raster band of index 1 when getting pixel value. Returning NULL +NOTICE: Could not find raster band of index 1 when getting pixel value. Returning NULL +NOTICE: Raster provided has no bands +NOTICE: Raster provided has no bands +NOTICE: Raster provided has no bands +NOTICE: Could not find raster band of index 1 when getting pixel value. Returning NULL +NOTICE: Could not find raster band of index 1 when getting pixel value. Returning NULL +NOTICE: Raster provided has no bands +NOTICE: Raster provided has no bands +NOTICE: Raster provided has no bands +NOTICE: Could not find raster band of index 1 when getting pixel value. Returning NULL +NOTICE: Could not find raster band of index 1 when getting pixel value. Returning NULL +NOTICE: Raster provided has no bands +NOTICE: Raster provided has no bands +NOTICE: Raster provided has no bands +NOTICE: Could not find raster band of index 1 when getting pixel value. Returning NULL +NOTICE: Could not find raster band of index 1 when getting pixel value. Returning NULL +NOTICE: Raster provided has no bands +NOTICE: Raster provided has no bands +NOTICE: Raster provided has no bands +NOTICE: Could not find raster band of index 1 when getting pixel value. Returning NULL +NOTICE: Could not find raster band of index 1 when getting pixel value. Returning NULL +NOTICE: Raster provided has no bands +NOTICE: Raster provided has no bands +NOTICE: Raster provided has no bands +NOTICE: Could not find raster band of index 1 when getting pixel value. Returning NULL +NOTICE: Could not find raster band of index 1 when getting pixel value. Returning NULL +NOTICE: Raster provided has no bands +NOTICE: Raster provided has no bands +NOTICE: Raster provided has no bands +NOTICE: Could not find raster band of index 1 when getting pixel value. Returning NULL +NOTICE: Could not find raster band of index 1 when getting pixel value. Returning NULL +NOTICE: Raster provided has no bands +NOTICE: Raster provided has no bands +NOTICE: Raster provided has no bands +NOTICE: Could not find raster band of index 1 when getting pixel value. Returning NULL +NOTICE: Could not find raster band of index 1 when getting pixel value. Returning NULL +NOTICE: Raster provided has no bands +NOTICE: Raster provided has no bands +NOTICE: Raster provided has no bands +NOTICE: Could not find raster band of index 1 when getting pixel value. Returning NULL +NOTICE: Could not find raster band of index 1 when getting pixel value. Returning NULL +NOTICE: Raster provided has no bands +NOTICE: Raster provided has no bands +NOTICE: Raster provided has no bands +NOTICE: Could not find raster band of index 1 when getting pixel value. Returning NULL +NOTICE: Could not find raster band of index 1 when getting pixel value. Returning NULL +NOTICE: Raster provided has no bands +NOTICE: Raster provided has no bands +NOTICE: Raster provided has no bands +NOTICE: Could not find raster band of index 1 when getting pixel value. Returning NULL +NOTICE: Could not find raster band of index 1 when getting pixel value. Returning NULL +NOTICE: Raster provided has no bands +NOTICE: Raster provided has no bands +NOTICE: Raster provided has no bands +NOTICE: Could not find raster band of index 1 when getting pixel value. Returning NULL +NOTICE: Could not find raster band of index 1 when getting pixel value. Returning NULL +NOTICE: Raster provided has no bands +NOTICE: Raster provided has no bands +NOTICE: Raster provided has no bands +NOTICE: Could not find raster band of index 1 when getting pixel value. Returning NULL +NOTICE: Could not find raster band of index 1 when getting pixel value. Returning NULL +NOTICE: Raster provided has no bands +NOTICE: Raster provided has no bands +NOTICE: Raster provided has no bands +NOTICE: Could not find raster band of index 1 when getting pixel value. Returning NULL +NOTICE: Could not find raster band of index 1 when getting pixel value. Returning NULL +NOTICE: Raster provided has no bands +NOTICE: Raster provided has no bands +NOTICE: Raster provided has no bands +NOTICE: Could not find raster band of index 1 when getting pixel value. Returning NULL +NOTICE: Could not find raster band of index 1 when getting pixel value. Returning NULL +NOTICE: Raster provided has no bands +NOTICE: Raster provided has no bands +NOTICE: Raster provided has no bands +NOTICE: Could not find raster band of index 1 when getting pixel value. Returning NULL +NOTICE: Could not find raster band of index 1 when getting pixel value. Returning NULL +NOTICE: Raster provided has no bands +NOTICE: Raster provided has no bands +NOTICE: Raster provided has no bands +NOTICE: Could not find raster band of index 1 when getting pixel value. Returning NULL +NOTICE: Could not find raster band of index 1 when getting pixel value. Returning NULL +NOTICE: Raster provided has no bands +NOTICE: Raster provided has no bands +NOTICE: Raster provided has no bands +NOTICE: Could not find raster band of index 1 when getting pixel value. Returning NULL +NOTICE: Could not find raster band of index 1 when getting pixel value. Returning NULL +0|1|INTERSECTION|0.000|0.000|2|2|1.000|1.000|0.000|0.000|0|1|32BF|t|0.000|1.000|1.000 +0|2|INTERSECTION|1.000|-1.000|1|2|1.000|1.000|0.000|0.000|0|1|32BF|t|0.000|1.000|1.000 +0|3|INTERSECTION|1.000|1.000|1|1|1.000|1.000|0.000|0.000|0|1|32BF|t|0.000|1.000|1.000 +0|4|INTERSECTION|0.000|0.000|0|0|0.000|0.000|0.000|0.000|0|0||||| +10|11|INTERSECTION|0.000|0.000|2|2|1.000|1.000|1.000|-1.000|0|1|32BF|t|0.000|1.000|1.000 +10|12|INTERSECTION|1.000|-1.000|2|2|1.000|1.000|1.000|-1.000|0|1|32BF|t|0.000|1.000|1.000 +10|13|INTERSECTION|1.000|1.000|2|1|1.000|1.000|1.000|-1.000|0|1|32BF|t|0.000|1.000|1.000 +10|14|INTERSECTION|0.000|0.000|0|0|0.000|0.000|0.000|0.000|0|0||||| +|0|INTERSECTION|0.000|0.000|0|0|0.000|0.000|0.000|0.000|0|0||||| +|1|INTERSECTION|0.000|0.000|0|0|0.000|0.000|0.000|0.000|0|0||||| +|2|INTERSECTION|0.000|0.000|0|0|0.000|0.000|0.000|0.000|0|0||||| +|3|INTERSECTION|0.000|0.000|0|0|0.000|0.000|0.000|0.000|0|0||||| +|4|INTERSECTION|0.000|0.000|0|0|0.000|0.000|0.000|0.000|0|0||||| +|10|INTERSECTION|0.000|0.000|0|0|0.000|0.000|0.000|0.000|0|0||||| +|11|INTERSECTION|0.000|0.000|0|0|0.000|0.000|0.000|0.000|0|0||||| +|12|INTERSECTION|0.000|0.000|0|0|0.000|0.000|0.000|0.000|0|0||||| +|13|INTERSECTION|0.000|0.000|0|0|0.000|0.000|0.000|0.000|0|0||||| +|14|INTERSECTION|0.000|0.000|0|0|0.000|0.000|0.000|0.000|0|0||||| +0||INTERSECTION|0.000|0.000|0|0|0.000|0.000|0.000|0.000|0|0||||| +1||INTERSECTION|0.000|0.000|0|0|0.000|0.000|0.000|0.000|0|0||||| +2||INTERSECTION|0.000|0.000|0|0|0.000|0.000|0.000|0.000|0|0||||| +3||INTERSECTION|0.000|0.000|0|0|0.000|0.000|0.000|0.000|0|0||||| +4||INTERSECTION|0.000|0.000|0|0|0.000|0.000|0.000|0.000|0|0||||| +10||INTERSECTION|0.000|0.000|0|0|0.000|0.000|0.000|0.000|0|0||||| +11||INTERSECTION|0.000|0.000|0|0|0.000|0.000|0.000|0.000|0|0||||| +12||INTERSECTION|0.000|0.000|0|0|0.000|0.000|0.000|0.000|0|0||||| +13||INTERSECTION|0.000|0.000|0|0|0.000|0.000|0.000|0.000|0|0||||| +14||INTERSECTION|0.000|0.000|0|0|0.000|0.000|0.000|0.000|0|0||||| +||INTERSECTION||||||||||||||| +0|1|UNION|-2.000|-2.000|4|4|1.000|1.000|0.000|0.000|0|1|32BF|t|0.000|1.000|1.500 +0|2|UNION|-2.000|-2.000|5|4|1.000|1.000|0.000|0.000|0|1|32BF|t|0.000|1.000| +0|3|UNION|-2.000|-2.000|5|5|1.000|1.000|0.000|0.000|0|1|32BF|t|0.000|1.000|4.000 +0|4|UNION|-2.000|-2.000|6|6|1.000|1.000|0.000|0.000|0|1|32BF|t|0.000|1.000|5.000 +10|11|UNION|-2.000|-2.000|4|4|1.000|1.000|1.000|-1.000|0|1|32BF|t|0.000|1.000|1.000 +10|12|UNION|-2.000|-2.000|4|4|1.000|1.000|1.000|-1.000|0|1|32BF|t|0.000|1.000|1.000 +10|13|UNION|-2.000|-2.000|4|5|1.000|1.000|1.000|-1.000|0|1|32BF|t|0.000|1.000| +10|14|UNION|-2.000|-2.000|4|6|1.000|1.000|1.000|-1.000|0|1|32BF|t|0.000|1.000| +|0|UNION|-2.000|-2.000|4|4|1.000|1.000|0.000|0.000|0|1|32BF|t|0.000|1.000|1.000 +|1|UNION|0.000|0.000|2|2|1.000|1.000|0.000|0.000|0|1|32BF|t|0.000|2.000|2.000 +|2|UNION|1.000|-1.000|2|2|1.000|1.000|0.000|0.000|0|1|32BF|t|0.000|3.000|3.000 +|3|UNION|1.000|1.000|2|2|1.000|1.000|0.000|0.000|0|1|32BF|t|0.000|4.000|4.000 +|4|UNION|2.000|2.000|2|2|1.000|1.000|0.000|0.000|0|1|32BF|t|0.000|5.000|5.000 +|10|UNION|-2.000|-2.000|4|4|1.000|1.000|1.000|-1.000|0|1|32BF|t|0.000|1.000|1.000 +|11|UNION|0.000|0.000|2|2|1.000|1.000|1.000|-1.000|0|1|32BF|t|0.000|2.000|2.000 +|12|UNION|1.000|-1.000|2|2|1.000|1.000|1.000|-1.000|0|1|32BF|t|0.000|3.000|3.000 +|13|UNION|1.000|1.000|2|2|1.000|1.000|1.000|-1.000|0|1|32BF|t|0.000|4.000|4.000 +|14|UNION|2.000|2.000|2|2|1.000|1.000|1.000|-1.000|0|1|32BF|t|0.000|5.000|5.000 +0||UNION|-2.000|-2.000|4|4|1.000|1.000|0.000|0.000|0|1|32BF|t|0.000|1.000|1.000 +1||UNION|0.000|0.000|2|2|1.000|1.000|0.000|0.000|0|1|32BF|t|0.000|2.000|2.000 +2||UNION|1.000|-1.000|2|2|1.000|1.000|0.000|0.000|0|1|32BF|t|0.000|3.000|3.000 +3||UNION|1.000|1.000|2|2|1.000|1.000|0.000|0.000|0|1|32BF|t|0.000|4.000|4.000 +4||UNION|2.000|2.000|2|2|1.000|1.000|0.000|0.000|0|1|32BF|t|0.000|5.000|5.000 +10||UNION|-2.000|-2.000|4|4|1.000|1.000|1.000|-1.000|0|1|32BF|t|0.000|1.000|1.000 +11||UNION|0.000|0.000|2|2|1.000|1.000|1.000|-1.000|0|1|32BF|t|0.000|2.000|2.000 +12||UNION|1.000|-1.000|2|2|1.000|1.000|1.000|-1.000|0|1|32BF|t|0.000|3.000|3.000 +13||UNION|1.000|1.000|2|2|1.000|1.000|1.000|-1.000|0|1|32BF|t|0.000|4.000|4.000 +14||UNION|2.000|2.000|2|2|1.000|1.000|1.000|-1.000|0|1|32BF|t|0.000|5.000|5.000 +||UNION||||||||||||||| +0|1|FIRST|-2.000|-2.000|4|4|1.000|1.000|0.000|0.000|0|1|32BF|t|0.000|1.000| +0|2|FIRST|-2.000|-2.000|4|4|1.000|1.000|0.000|0.000|0|1|32BF|t|0.000|1.000|1.000 +0|3|FIRST|-2.000|-2.000|4|4|1.000|1.000|0.000|0.000|0|1|32BF|t|0.000|1.000| +0|4|FIRST|-2.000|-2.000|4|4|1.000|1.000|0.000|0.000|0|1|32BF|t|0.000|1.000|1.000 +10|11|FIRST|-2.000|-2.000|4|4|1.000|1.000|1.000|-1.000|0|1|32BF|t|0.000|1.000|1.000 +10|12|FIRST|-2.000|-2.000|4|4|1.000|1.000|1.000|-1.000|0|1|32BF|t|0.000|1.000|1.000 +10|13|FIRST|-2.000|-2.000|4|5|1.000|1.000|1.000|-1.000|0|1|32BF|t|0.000|1.000| +10|14|FIRST|-2.000|-2.000|4|6|1.000|1.000|1.000|-1.000|0|1|32BF|t|0.000|1.000| +|0|FIRST|-2.000|-2.000|4|4|1.000|1.000|0.000|0.000|0|1|32BF|t|0.000|1.000|1.000 +|1|FIRST|0.000|0.000|2|2|1.000|1.000|0.000|0.000|0|1|32BF|t|0.000|2.000|2.000 +|2|FIRST|1.000|-1.000|2|2|1.000|1.000|0.000|0.000|0|1|32BF|t|0.000|3.000|3.000 +|3|FIRST|1.000|1.000|2|2|1.000|1.000|0.000|0.000|0|1|32BF|t|0.000|4.000|4.000 +|4|FIRST|2.000|2.000|2|2|1.000|1.000|0.000|0.000|0|1|32BF|t|0.000|5.000|5.000 +|10|FIRST|-2.000|-2.000|4|4|1.000|1.000|1.000|-1.000|0|1|32BF|t|0.000|1.000|1.000 +|11|FIRST|0.000|0.000|2|2|1.000|1.000|1.000|-1.000|0|1|32BF|t|0.000|2.000|2.000 +|12|FIRST|1.000|-1.000|2|2|1.000|1.000|1.000|-1.000|0|1|32BF|t|0.000|3.000|3.000 +|13|FIRST|1.000|1.000|2|2|1.000|1.000|1.000|-1.000|0|1|32BF|t|0.000|4.000|4.000 +|14|FIRST|2.000|2.000|2|2|1.000|1.000|1.000|-1.000|0|1|32BF|t|0.000|5.000|5.000 +0||FIRST|-2.000|-2.000|4|4|1.000|1.000|0.000|0.000|0|1|32BF|t|0.000|1.000|1.000 +1||FIRST|0.000|0.000|2|2|1.000|1.000|0.000|0.000|0|1|32BF|t|0.000|2.000|2.000 +2||FIRST|1.000|-1.000|2|2|1.000|1.000|0.000|0.000|0|1|32BF|t|0.000|3.000|3.000 +3||FIRST|1.000|1.000|2|2|1.000|1.000|0.000|0.000|0|1|32BF|t|0.000|4.000|4.000 +4||FIRST|2.000|2.000|2|2|1.000|1.000|0.000|0.000|0|1|32BF|t|0.000|5.000|5.000 +10||FIRST|-2.000|-2.000|4|4|1.000|1.000|1.000|-1.000|0|1|32BF|t|0.000|1.000|1.000 +11||FIRST|0.000|0.000|2|2|1.000|1.000|1.000|-1.000|0|1|32BF|t|0.000|2.000|2.000 +12||FIRST|1.000|-1.000|2|2|1.000|1.000|1.000|-1.000|0|1|32BF|t|0.000|3.000|3.000 +13||FIRST|1.000|1.000|2|2|1.000|1.000|1.000|-1.000|0|1|32BF|t|0.000|4.000|4.000 +14||FIRST|2.000|2.000|2|2|1.000|1.000|1.000|-1.000|0|1|32BF|t|0.000|5.000|5.000 +||FIRST||||||||||||||| +0|1|SECOND|0.000|0.000|2|2|1.000|1.000|0.000|0.000|0|1|32BF|t|0.000|| +0|2|SECOND|1.000|-1.000|2|2|1.000|1.000|0.000|0.000|0|1|32BF|t|0.000||3.000 +0|3|SECOND|1.000|1.000|2|2|1.000|1.000|0.000|0.000|0|1|32BF|t|0.000||4.000 +0|4|SECOND|2.000|2.000|2|2|1.000|1.000|0.000|0.000|0|1|32BF|t|0.000|5.000|5.000 +10|11|SECOND|0.000|0.000|2|2|1.000|1.000|1.000|-1.000|0|1|32BF|t|0.000|| +10|12|SECOND|1.000|-1.000|2|2|1.000|1.000|1.000|-1.000|0|1|32BF|t|0.000|| +10|13|SECOND|1.000|1.000|2|2|1.000|1.000|1.000|-1.000|0|1|32BF|t|0.000||4.000 +10|14|SECOND|2.000|2.000|2|2|1.000|1.000|1.000|-1.000|0|1|32BF|t|0.000|5.000|5.000 +|0|SECOND|-2.000|-2.000|4|4|1.000|1.000|0.000|0.000|0|1|32BF|t|0.000|1.000|1.000 +|1|SECOND|0.000|0.000|2|2|1.000|1.000|0.000|0.000|0|1|32BF|t|0.000|2.000|2.000 +|2|SECOND|1.000|-1.000|2|2|1.000|1.000|0.000|0.000|0|1|32BF|t|0.000|3.000|3.000 +|3|SECOND|1.000|1.000|2|2|1.000|1.000|0.000|0.000|0|1|32BF|t|0.000|4.000|4.000 +|4|SECOND|2.000|2.000|2|2|1.000|1.000|0.000|0.000|0|1|32BF|t|0.000|5.000|5.000 +|10|SECOND|-2.000|-2.000|4|4|1.000|1.000|1.000|-1.000|0|1|32BF|t|0.000|1.000|1.000 +|11|SECOND|0.000|0.000|2|2|1.000|1.000|1.000|-1.000|0|1|32BF|t|0.000|2.000|2.000 +|12|SECOND|1.000|-1.000|2|2|1.000|1.000|1.000|-1.000|0|1|32BF|t|0.000|3.000|3.000 +|13|SECOND|1.000|1.000|2|2|1.000|1.000|1.000|-1.000|0|1|32BF|t|0.000|4.000|4.000 +|14|SECOND|2.000|2.000|2|2|1.000|1.000|1.000|-1.000|0|1|32BF|t|0.000|5.000|5.000 +0||SECOND||||||||||||||| +1||SECOND||||||||||||||| +2||SECOND||||||||||||||| +3||SECOND||||||||||||||| +4||SECOND||||||||||||||| +10||SECOND||||||||||||||| +11||SECOND||||||||||||||| +12||SECOND||||||||||||||| +13||SECOND||||||||||||||| +14||SECOND||||||||||||||| +||SECOND|||||||||||||||