From f3122464c11f8cb000048c155b38b22a713a18db Mon Sep 17 00:00:00 2001 From: Bborie Park Date: Mon, 25 Mar 2013 03:14:49 +0000 Subject: [PATCH] Addition of ST_FromGDALRaster(). This is what happens on a long airplane flight! Ticket #2123. git-svn-id: http://svn.osgeo.org/postgis/trunk@11202 b70326c6-7e19-0410-871a-916f4a2858ee --- NEWS | 1 + doc/reference_raster.xml | 73 ++++++++++++++- raster/rt_core/rt_api.c | 12 ++- raster/rt_pg/rt_pg.c | 90 +++++++++++++++++++ raster/rt_pg/rtpostgis.sql.in | 8 ++ raster/test/regress/Makefile.in | 1 + raster/test/regress/rt_fromgdalraster.sql | 77 ++++++++++++++++ .../test/regress/rt_fromgdalraster_expected | 3 + 8 files changed, 260 insertions(+), 5 deletions(-) create mode 100644 raster/test/regress/rt_fromgdalraster.sql create mode 100644 raster/test/regress/rt_fromgdalraster_expected diff --git a/NEWS b/NEWS index 5180f487a..2c52e32be 100644 --- a/NEWS +++ b/NEWS @@ -74,6 +74,7 @@ PostGIS 2.1.0 - lwgeom_from_geojson in liblwgeom (Sandro Santilli / Vizzuality) - #1687, ST_Simplify for TopoGeometry (Sandro Santilli / Vizzuality) - #2228, TopoJSON output for TopoGeometry (Sandro Santilli / Vizzuality) + - #2123, ST_FromGDALRaster * Enhancements * - #823, tiger geocoder: Make loader_generate_script download portion diff --git a/doc/reference_raster.xml b/doc/reference_raster.xml index bc04452b6..cb5ed90ae 100644 --- a/doc/reference_raster.xml +++ b/doc/reference_raster.xml @@ -1772,7 +1772,78 @@ FROM baz; - + + + + ST_FromGDALRaster + Returns a raster from a supported GDAL raster file. + + + + + + raster ST_FromGDALRaster + bytea gdaldata + integer srid=NULL + + + + + + + Description + + + Returns a raster from a supported GDAL raster file. gdaldata is of type bytea and should be the contents of the GDAL raster file. + + + + If srid is NULL, the function will try to autmatically assign the SRID from the GDAL raster. If srid is provided, the value provided will override any automatically assigned SRID. + + + Availability: 2.1.0 + + + + Examples + + +WITH foo AS ( + SELECT ST_AsPNG(ST_AddBand(ST_AddBand(ST_AddBand(ST_MakeEmptyRaster(2, 2, 0, 0, 0.1, -0.1, 0, 0, 4326), 1, '8BUI', 1, 0), 2, '8BUI', 2, 0), 3, '8BUI', 3, 0)) AS png +), +bar AS ( + SELECT 1 AS rid, ST_FromGDALRaster(png) AS rast FROM foo + UNION ALL + SELECT 2 AS rid, ST_FromGDALRaster(png, 3310) AS rast FROM foo +) +SELECT + rid, + ST_Metadata(rast) AS metadata, + ST_SummaryStats(rast, 1) AS stats1, + ST_SummaryStats(rast, 2) AS stats2, + ST_SummaryStats(rast, 3) AS stats3 +FROM bar +ORDER BY rid; + + rid | metadata | stats1 | stats2 | stats3 +-----+---------------------------+---------------+---------------+---------------- + 1 | (0,0,2,2,1,-1,0,0,0,3) | (4,4,1,0,1,1) | (4,8,2,0,2,2) | (4,12,3,0,3,3) + 2 | (0,0,2,2,1,-1,0,0,3310,3) | (4,4,1,0,1,1) | (4,8,2,0,2,2) | (4,12,3,0,3,3) +(2 rows) + + + + + + See Also + + + + + + + + Raster Accessors diff --git a/raster/rt_core/rt_api.c b/raster/rt_core/rt_api.c index 34b3aa4fd..195b8e720 100644 --- a/raster/rt_core/rt_api.c +++ b/raster/rt_core/rt_api.c @@ -9138,10 +9138,14 @@ rt_raster_from_gdal_dataset(GDALDatasetH ds) { /* get raster attributes */ cplerr = GDALGetGeoTransform(ds, gt); - if (cplerr != CE_None) { - rterror("rt_raster_from_gdal_dataset: Unable to get geotransformation"); - rt_raster_destroy(rast); - return NULL; + if (GDALGetGeoTransform(ds, gt) != CE_None) { + RASTER_DEBUG(4, "Using default geotransform matrix (0, 1, 0, 0, 0, -1)"); + gt[0] = 0; + gt[1] = 1; + gt[2] = 0; + gt[3] = 0; + gt[4] = 0; + gt[5] = -1; } /* apply raster attributes */ diff --git a/raster/rt_pg/rt_pg.c b/raster/rt_pg/rt_pg.c index 7445578f3..ee11e2f38 100644 --- a/raster/rt_pg/rt_pg.c +++ b/raster/rt_pg/rt_pg.c @@ -310,6 +310,9 @@ Datum RASTER_valueCountCoverage(PG_FUNCTION_ARGS); /* reclassify specified bands of a raster */ Datum RASTER_reclass(PG_FUNCTION_ARGS); +/* convert GDAL raster to raster */ +Datum RASTER_fromGDALRaster(PG_FUNCTION_ARGS); + /* convert raster to GDAL raster */ Datum RASTER_asGDALRaster(PG_FUNCTION_ARGS); Datum RASTER_getGDALDrivers(PG_FUNCTION_ARGS); @@ -10377,6 +10380,93 @@ Datum RASTER_reclass(PG_FUNCTION_ARGS) { PG_RETURN_POINTER(pgrtn); } +/** + * Returns raster from GDAL raster; + */ +PG_FUNCTION_INFO_V1(RASTER_fromGDALRaster); +Datum RASTER_fromGDALRaster(PG_FUNCTION_ARGS) +{ + bytea *bytea_data; + uint8_t *data; + int data_len = 0; + VSILFILE *vsifp = NULL; + GDALDatasetH hdsSrc; + int srid = -1; /* -1 for NULL */ + + rt_pgraster *pgraster = NULL; + rt_raster raster; + + /* NULL if NULL */ + if (PG_ARGISNULL(0)) + PG_RETURN_NULL(); + + /* get data */ + bytea_data = (bytea *) PG_GETARG_BYTEA_P(0); + data = (uint8_t *) VARDATA(bytea_data); + data_len = VARSIZE(bytea_data) - VARHDRSZ; + + /* process srid */ + /* NULL srid means try to determine SRID from bytea */ + if (!PG_ARGISNULL(1)) + srid = clamp_srid(PG_GETARG_INT32(1)); + + /* create memory "file" */ + vsifp = VSIFileFromMemBuffer("/vsimem/in.dat", data, data_len, FALSE); + if (vsifp == NULL) { + elog(ERROR, "RASTER_fromGDALRaster: Could not load bytea into memory file for use by GDAL"); + PG_FREE_IF_COPY(bytea_data, 0); + PG_RETURN_NULL(); + } + + /* register all GDAL drivers */ + rt_util_gdal_register_all(); + + /* open GDAL raster */ + hdsSrc = GDALOpenShared("/vsimem/in.dat", GA_ReadOnly); + if (hdsSrc == NULL) { + elog(ERROR, "RASTER_fromGDALRaster: Could not open bytea with GDAL. Check that the bytea is of a GDAL supported format"); + VSIFCloseL(vsifp); + PG_FREE_IF_COPY(bytea_data, 0); + PG_RETURN_NULL(); + } + +#if POSTGIS_DEBUG_LEVEL > 3 + { + GDALDriverH hdrv = GDALGetDatasetDriver(hdsSrc); + + POSTGIS_RT_DEBUGF(4, "Input GDAL Raster info: %s, (%d x %d)", + GDALGetDriverShortName(hdrv), + GDALGetRasterXSize(hdsSrc), + GDALGetRasterYSize(hdsSrc) + ); + } +#endif + + /* convert GDAL raster to raster */ + raster = rt_raster_from_gdal_dataset(hdsSrc); + + GDALClose(hdsSrc); + VSIFCloseL(vsifp); + PG_FREE_IF_COPY(bytea_data, 0); + + if (raster == NULL) { + elog(ERROR, "RASTER_fromGDALRaster: Could not convert GDAL raster to raster"); + PG_RETURN_NULL(); + } + + /* apply SRID if set */ + if (srid != -1) + rt_raster_set_srid(raster, srid); + + pgraster = rt_raster_serialize(raster); + rt_raster_destroy(raster); + if (!pgraster) + PG_RETURN_NULL(); + + SET_VARSIZE(pgraster, pgraster->size); + PG_RETURN_POINTER(pgraster); +} + /** * Returns formatted GDAL raster as bytea object of raster */ diff --git a/raster/rt_pg/rtpostgis.sql.in b/raster/rt_pg/rtpostgis.sql.in index f3af59f72..de9ae3c1a 100644 --- a/raster/rt_pg/rtpostgis.sql.in +++ b/raster/rt_pg/rtpostgis.sql.in @@ -1510,6 +1510,14 @@ CREATE OR REPLACE FUNCTION st_reclass(rast raster, reclassexpr text, pixeltype t AS $$ SELECT st_reclass($1, ROW(1, $2, $3, NULL)) $$ LANGUAGE 'sql' IMMUTABLE STRICT; +----------------------------------------------------------------------- +-- ST_FromGDALRaster +----------------------------------------------------------------------- +CREATE OR REPLACE FUNCTION st_fromgdalraster(gdaldata bytea, srid integer DEFAULT NULL) + RETURNS raster + AS 'MODULE_PATHNAME', 'RASTER_fromGDALRaster' + LANGUAGE 'c' IMMUTABLE; + ----------------------------------------------------------------------- -- ST_AsGDALRaster and supporting functions ----------------------------------------------------------------------- diff --git a/raster/test/regress/Makefile.in b/raster/test/regress/Makefile.in index 3df8d2ab3..a003bde63 100644 --- a/raster/test/regress/Makefile.in +++ b/raster/test/regress/Makefile.in @@ -83,6 +83,7 @@ TEST_BANDPROPS = \ TEST_UTILITY = \ rt_utility \ + rt_fromgdalraster \ rt_asgdalraster \ rt_astiff \ rt_asjpeg \ diff --git a/raster/test/regress/rt_fromgdalraster.sql b/raster/test/regress/rt_fromgdalraster.sql new file mode 100644 index 000000000..96b049d33 --- /dev/null +++ b/raster/test/regress/rt_fromgdalraster.sql @@ -0,0 +1,77 @@ +/* +Driver: PNG/Portable Network Graphics +Files: testraster.png +Size is 90, 90 +Coordinate System is `' +Metadata: + Comment=Created with GIMP +Image Structure Metadata: + INTERLEAVE=PIXEL +Corner Coordinates: +Upper Left ( 0.0, 0.0) +Lower Left ( 0.0, 90.0) +Upper Right ( 90.0, 0.0) +Lower Right ( 90.0, 90.0) +Center ( 45.0, 45.0) +Band 1 Block=90x1 Type=Byte, ColorInterp=Red +Band 2 Block=90x1 Type=Byte, ColorInterp=Green +Band 3 Block=90x1 Type=Byte, ColorInterp=Blue +*/ +WITH foo AS ( + SELECT ST_FromGDALRaster(E'\\x89504e470d0a1a0a0000000d494844520000005a0000005a0802000000b7cad655000000017352474200aece1ce9000000097048597300000b1300000b1301009a9c180000000774494d4507dc020e16332b734b11150000001d69545874436f6d6d656e7400000000004372656174656420776974682047494d50642e6507000000c14944415478daedd8410ac030080440edffff9cfe607b6a0c61bc8a974161b157c55aa9dfdd71f4c4d9aa34fb94c28103070e1c3870e018a8feca70b60387c28103070e1c3870fc9f4a737b5df72bcdb3b603070e1c3870e0c081632695fa95da0e1c3870e0c08103078ef9549adb7ea58e45e1c08103070e1c3876a452bf52db8103070e1c3870e0c0319f4a73dbafd4b1281c3870e0c08103c78e54ea576a3b70e0c08103070e1c38e65369c550baeabe5fa9ed702c3870e0c08103078ee3ea057386489dc63798000000000049454e44ae426082'::bytea) AS rast +) +SELECT + ST_Metadata(rast), + ST_SummaryStats(rast, 1), + ST_SummaryStats(rast, 2), + ST_SummaryStats(rast, 3) +FROM foo; + +/* +Driver: GTiff/GeoTIFF +Files: test4269.tif +Size is 10, 10 +Coordinate System is: +GEOGCS["NAD83", + DATUM["North_American_Datum_1983", + SPHEROID["GRS 1980",6378137,298.2572221010002, + AUTHORITY["EPSG","7019"]], + TOWGS84[0,0,0,0,0,0,0], + AUTHORITY["EPSG","6269"]], + PRIMEM["Greenwich",0], + UNIT["degree",0.0174532925199433], + AUTHORITY["EPSG","4269"]] +Origin = (-168.000000000000000,85.000000000000000) +Pixel Size = (0.083000000000000,-0.083000000000000) +Metadata: + AREA_OR_POINT=Area + TIFFTAG_DATETIME=2012:03:02 09:59:31 + TIFFTAG_RESOLUTIONUNIT=2 (pixels/inch) + TIFFTAG_SOFTWARE=Adobe Photoshop CS Windows + TIFFTAG_XRESOLUTION=96 + TIFFTAG_YRESOLUTION=96 +Image Structure Metadata: + INTERLEAVE=PIXEL +Corner Coordinates: +Upper Left (-168.0000000, 85.0000000) (168d 0' 0.00"W, 85d 0' 0.00"N) +Lower Left (-168.0000000, 84.1700000) (168d 0' 0.00"W, 84d10'12.00"N) +Upper Right (-167.1700000, 85.0000000) (167d10'12.00"W, 85d 0' 0.00"N) +Lower Right (-167.1700000, 84.1700000) (167d10'12.00"W, 84d10'12.00"N) +Center (-167.5850000, 84.5850000) (167d35' 6.00"W, 84d35' 6.00"N) +Band 1 Block=10x10 Type=Byte, ColorInterp=Red +Band 2 Block=10x10 Type=Byte, ColorInterp=Green +Band 3 Block=10x10 Type=Byte, ColorInterp=Blue +*/ +WITH foo AS ( + SELECT 1 AS rid, ST_FromGDALRaster(E'\\x49492a0008000000150000010300010000000a00000001010300010000000a00000002010300030000001a01000003010300010000000100000006010300010000000200000011010400010000001502000015010300010000000300000016010300010000000a00000017010400010000002c0100001a010500010000000a0100001b01050001000000120100001c0103000100000001000000280103000100000002000000310102001b0000003a0100003201020014000000260100005301030003000000200100000e830c00030000005601000082840c00060000006e010000af870300240000009e010000b0870c0005000000e6010000b1870200070000000e0200000000000060000000010000006000000001000000080008000800010001000100323031323a30333a30322030393a35393a33310041646f62652050686f746f73686f702043532057696e646f77730000736891ed7c3fb53f736891ed7c3fb53f000000000000000000000000000000000000000000000000000000000000000000000000000065c000000000004055400000000000000000010001000000080000040000010002000104000001000100000800000100ad100108b187060000000608000001008e230908b087010001000b08b087010000000e08b08703000200a8f9eb941da4724000000040a65458410000000000000000000000000000000000000000000000004e414438337c00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000'::bytea) AS rast + UNION ALL + SELECT 2 AS rid, ST_FromGDALRaster(E'\\x49492a0008000000150000010300010000000a00000001010300010000000a00000002010300030000001a01000003010300010000000100000006010300010000000200000011010400010000001502000015010300010000000300000016010300010000000a00000017010400010000002c0100001a010500010000000a0100001b01050001000000120100001c0103000100000001000000280103000100000002000000310102001b0000003a0100003201020014000000260100005301030003000000200100000e830c00030000005601000082840c00060000006e010000af870300240000009e010000b0870c0005000000e6010000b1870200070000000e0200000000000060000000010000006000000001000000080008000800010001000100323031323a30333a30322030393a35393a33310041646f62652050686f746f73686f702043532057696e646f77730000736891ed7c3fb53f736891ed7c3fb53f000000000000000000000000000000000000000000000000000000000000000000000000000065c000000000004055400000000000000000010001000000080000040000010002000104000001000100000800000100ad100108b187060000000608000001008e230908b087010001000b08b087010000000e08b08703000200a8f9eb941da4724000000040a65458410000000000000000000000000000000000000000000000004e414438337c00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000'::bytea, 4326) AS rast +) +SELECT + rid, + ST_Metadata(rast), + ST_SummaryStats(rast, 1), + ST_SummaryStats(rast, 2), + ST_SummaryStats(rast, 3) +FROM foo +ORDER BY rid; diff --git a/raster/test/regress/rt_fromgdalraster_expected b/raster/test/regress/rt_fromgdalraster_expected new file mode 100644 index 000000000..e0a010211 --- /dev/null +++ b/raster/test/regress/rt_fromgdalraster_expected @@ -0,0 +1,3 @@ +(0,0,90,90,1,-1,0,0,0,3)|(8100,1479000,182.592592592593,114.982851945091,0,255)|(8100,1479000,182.592592592593,114.982851945091,0,255)|(8100,1453500,179.444444444444,116.438931167192,0,255) +1|(-168,85,10,10,0.083,-0.083,0,0,4269,3)|(100,0,0,0,0,0)|(100,0,0,0,0,0)|(100,0,0,0,0,0) +2|(-168,85,10,10,0.083,-0.083,0,0,4326,3)|(100,0,0,0,0,0)|(100,0,0,0,0,0)|(100,0,0,0,0,0) -- 2.50.1