From 5018f99cb5ab50d735d789702bde8d616105b47c Mon Sep 17 00:00:00 2001 From: Sandro Santilli Date: Wed, 1 Feb 2012 09:20:44 +0000 Subject: [PATCH] Implement ST_PixelAsPolygon in C, provide a core API entry for it The API entry point is expected to be useful for #1519 git-svn-id: http://svn.osgeo.org/postgis/trunk@8993 b70326c6-7e19-0410-871a-916f4a2858ee --- raster/rt_core/rt_api.c | 47 ++++++++++++++++++++++++++++ raster/rt_core/rt_api.h | 16 ++++++++++ raster/rt_pg/rt_pg.c | 55 +++++++++++++++++++++++++++++++++ raster/rt_pg/rtpostgis.sql.in.c | 54 +++++--------------------------- 4 files changed, 126 insertions(+), 46 deletions(-) diff --git a/raster/rt_core/rt_api.c b/raster/rt_core/rt_api.c index 5895577f1..8898580e4 100644 --- a/raster/rt_core/rt_api.c +++ b/raster/rt_core/rt_api.c @@ -9652,3 +9652,50 @@ rt_raster_from_two_rasters( *err = 1; return raster; } + + +LWPOLY* +rt_raster_pixel_as_polygon(rt_raster rast, int x, int y) +{ + double scale_x, scale_y; + double skew_x, skew_y; + double ul_x, ul_y; + int srid; + POINTARRAY **points; + POINT4D p, p0; + LWPOLY *poly; + + scale_x = rt_raster_get_x_scale(rast); + scale_y = rt_raster_get_y_scale(rast); + skew_x = rt_raster_get_x_skew(rast); + skew_y = rt_raster_get_y_skew(rast); + ul_x = rt_raster_get_x_offset(rast); + ul_y = rt_raster_get_y_offset(rast); + srid = rt_raster_get_srid(rast); + + points = rtalloc(sizeof(POINTARRAY *)*1); + points[0] = ptarray_construct(0, 0, 5); + + p0.x = scale_x * x + skew_x * y + ul_x; + p0.y = scale_y * y + skew_y * x + ul_y; + ptarray_set_point4d(points[0], 0, &p0); + + p.x = p0.x + scale_x; + p.y = p0.y + skew_y; + ptarray_set_point4d(points[0], 1, &p); + + p.x = p0.x + scale_x + skew_x; + p.y = p0.y + scale_y + skew_y; + ptarray_set_point4d(points[0], 2, &p); + + p.x = p0.x + skew_x; + p.y = p0.y + scale_y; + ptarray_set_point4d(points[0], 3, &p); + + /* close it */ + ptarray_set_point4d(points[0], 4, &p0); + + poly = lwpoly_construct(srid, NULL, 1, points); + + return poly; +} diff --git a/raster/rt_core/rt_api.h b/raster/rt_core/rt_api.h index 9d6ab6ad9..863faf926 100644 --- a/raster/rt_core/rt_api.h +++ b/raster/rt_core/rt_api.h @@ -949,6 +949,22 @@ int rt_raster_geopoint_to_cell(rt_raster raster, */ LWPOLY* rt_raster_get_convex_hull(rt_raster raster); +/** + * Get a raster pixel as a polygon. + * + * The pixel shape is a 4 vertices (5 to be closed) single + * ring polygon bearing the raster's rotation + * and using projection coordinates + * + * @param raster : the raster to get pixel from + * @param x : the column number + * @param y : the row number + * + * @return the pixel polygon, or NULL on error. + * + */ +LWPOLY* rt_raster_pixel_as_polygon(rt_raster raster, int x, int y); + /** * Returns a set of "geomval" value, one for each group of pixel diff --git a/raster/rt_pg/rt_pg.c b/raster/rt_pg/rt_pg.c index fc22dc340..64acf7e6c 100644 --- a/raster/rt_pg/rt_pg.c +++ b/raster/rt_pg/rt_pg.c @@ -198,6 +198,9 @@ Datum RASTER_getPixelValue(PG_FUNCTION_ARGS); /* Set pixel value */ Datum RASTER_setPixelValue(PG_FUNCTION_ARGS); +/* Get pixel geographical shape */ +Datum RASTER_getPixelPolygon(PG_FUNCTION_ARGS); + /* Raster and band creation */ Datum RASTER_makeEmpty(PG_FUNCTION_ARGS); Datum RASTER_addband(PG_FUNCTION_ARGS); @@ -2209,6 +2212,58 @@ Datum RASTER_setPixelValue(PG_FUNCTION_ARGS) PG_RETURN_POINTER(pgraster); } +/** + * Return the geographical shape of a single pixel. + * Pixel location is specified by + * X,Y coordinates (X <= RT_Width(raster) and Y <= RT_Height(raster)). + * + */ +PG_FUNCTION_INFO_V1(RASTER_getPixelPolygon); +Datum RASTER_getPixelPolygon(PG_FUNCTION_ARGS) +{ + rt_pgraster *pgraster = NULL; + rt_raster raster = NULL; + LWPOLY *poly; + int32_t x = 0; + int32_t y = 0; + GSERIALIZED* gser = NULL; + + + /* Deserialize raster */ + pgraster = (rt_pgraster *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); + + x = PG_GETARG_INT32(1); + + y = PG_GETARG_INT32(2); + + POSTGIS_RT_DEBUGF(3, "Pixel coordinates (%d, %d)", x, y); + + raster = rt_raster_deserialize(pgraster, FALSE); + if (!raster) { + elog(ERROR, "RASTER_getPixelValue: Could not deserialize raster"); + PG_RETURN_NULL(); + } + + /* Fetch pixel shape */ + poly = rt_raster_pixel_as_polygon(raster, x - 1, y - 1); + if (!poly) { + elog(ERROR, "RASTER_getPixelPolygon: could not get raster's pixel polygon"); + rt_raster_destroy(raster); + PG_RETURN_NULL(); + } + + { + size_t gser_size; + gser = gserialized_from_lwgeom(lwpoly_as_lwgeom(poly), 0, &gser_size); + SET_VARSIZE(gser, gser_size); + } + + rt_raster_destroy(raster); + lwfree(poly); + + PG_RETURN_POINTER(gser); +} + /** * Add a band to the given raster at the given position. */ diff --git a/raster/rt_pg/rtpostgis.sql.in.c b/raster/rt_pg/rtpostgis.sql.in.c index f000debed..3c517b84d 100644 --- a/raster/rt_pg/rtpostgis.sql.in.c +++ b/raster/rt_pg/rtpostgis.sql.in.c @@ -2419,54 +2419,16 @@ CREATE OR REPLACE FUNCTION st_polygon(rast raster, band integer DEFAULT 1) $$ LANGUAGE 'SQL' IMMUTABLE STRICT; -CREATE OR REPLACE FUNCTION st_pixelaspolygon(rast raster, band integer, x integer, y integer) - RETURNS geometry AS - $$ - DECLARE - w int; - h int; - scale_x float8; - scale_y float8; - skew_x float8; - skew_y float8; - ul_x float8; - ul_y float8; - sr_id int; - x1 float8; - y1 float8; - x2 float8; - y2 float8; - x3 float8; - y3 float8; - x4 float8; - y4 float8; - BEGIN - SELECT scalex, scaley, skewx, skewy, upperleftx, upperlefty, srid INTO scale_x, scale_y, skew_x, skew_y, ul_x, ul_y, sr_id FROM ST_Metadata(rast); - x1 := scale_x * (x - 1) + skew_x * (y - 1) + ul_x; - y1 := scale_y * (y - 1) + skew_y * (x - 1) + ul_y; - x2 := x1 + scale_x; - y2 := y1 + skew_y; - x3 := x1 + scale_x + skew_x; - y3 := y1 + scale_y + skew_y; - x4 := x1 + skew_x; - y4 := y1 + scale_y; - RETURN st_setsrid(st_makepolygon(st_makeline(ARRAY[st_makepoint(x1, y1), - st_makepoint(x2, y2), - st_makepoint(x3, y3), - st_makepoint(x4, y4), - st_makepoint(x1, y1)] - ) - ), - sr_id - ); - END; - $$ - LANGUAGE 'plpgsql' IMMUTABLE STRICT; - CREATE OR REPLACE FUNCTION st_pixelaspolygon(rast raster, x integer, y integer) + RETURNS geometry + AS 'MODULE_PATHNAME','RASTER_getPixelPolygon' + LANGUAGE 'C' IMMUTABLE STRICT; + +-- TODO: deprecate +CREATE OR REPLACE FUNCTION st_pixelaspolygon(rast raster, band integer, x integer, y integer) RETURNS geometry AS - $$ SELECT st_pixelaspolygon($1, 1, $2, $3) $$ - LANGUAGE SQL IMMUTABLE STRICT; + $$ SELECT st_pixelaspolygon($1, $3, $4); $$ + LANGUAGE 'sql' IMMUTABLE STRICT; ----------------------------------------------------------------------- -- ST_PixelAsPolygons -- 2.40.0