[raster] ST_AsGDALRaster and the helper functions ST_getGDALDrivers and ST_srtext

[Bborie Park]

What follows is a proposal on the specifics of ST_AsGDALRaster and two supporting functions.

This is a generic interface to outputting a supported and installed GDAL raster:

ST_AsGDALRaster(rast raster, format text, options text[], srs text) → bytea

This is the most generic and GDAL-specific method to convert a raster to a GDAL raster. All other version of ST_AsGDALRaster and other format specific functions (ST_AsJPEG, ST_AsTIFF and ST_AsPNG) are all wrappers around this function. Reference information for the format and options arguments of a particular format are specified at: ​http://gdal.org/formats_list.html. The arguments specified are:

format: the GDAL format code. e.g. GTiff, JPEG, PNG

options: the GDAL creation options found in the Creation Options section of a specified format. e.g. COMPRESS=JPEG, JPEG_QUALITY=90

srs: the user-specified OGC WKT or the proj4 text for a spatial reference to embed in the GDAL raster. Not all formats support embedding this information. e.g. the non-empty value for the srtext or proj4text column from the spatial_ref_sys table.

ST_AsGDALRaster(rast, 'GTiff', ARRAY['COMPRESS=JPEG', 'JPEG_QUALITY=90'], '+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')

ST_AsGDALRaster(rast, 'GTiff', ARRAY['COMPRESS=JPEG', 'JPEG_QUALITY=90'], '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]]')

ST_AsGDALRaster(rast raster, format text, options text[]) → bytea

This one removes the user-specified srs argument. The output GDAL raster's spatial reference will be set to the same as the input raster, if possible.

ST_AsGDALRaster(rast, 'JPEG', ARRAY!['QUALITY=50'])

ST_AsGDALRaster(rast, 'PNG', ARRAY!['ZLEVEL=7'])

ST_AsGDALRaster(rast raster, format text) → bytea

The simplest implementation of this function. Since the options argument has been removed, the output GDAL raster will be created with default options. Like the prior function, the spatial reference of the GDAL raster will be set to the same as the input raster.

ST_AsGDALRaster(rast, 'JPEG')

The two helper functions are:

ST_getGDALDrivers() → set of record

As each GDAL installation may be different and ST_AsGDALRaster can be used to support formats other than GTiff, JPEG and PNG, a method is needed to expose to the end user the possible GDAL formats capable of being exported. This function will output the following columns.

idx: the internal GDAL index number

short_name: the GDAL format code. This is the value to pass to the format paramenter of ST_AsGDALRaster

long_name: the full name of the GDAL format

create_options: the creation options available for the format as an XML string.

SELECT * FROM ST_getGDALDrivers()

The formats outputted from ST_getGDALDrivers have been filtered to only those that the GDAL capabilities CreateCopy and Virtual IO support.

Should the GDAL raster process be capable of supporting the GDAL capability Create? As the GDAL raster process writes nothing to a file in the filesystem (via Virtual IO), should there be support for writing the output GDAL raster temporarily to the filesystem? If so, how is it done in other PostgreSQL extensions in a secure manner?

ST_srtext(rast raster) → text

A helper function to get the value of column srtext or proj4text for a raster with an SRID. By default, the srtext is returned. If srtext is not available but proj4text is, the proj4text is returned.

SELECT ST_srtext(rast);

This function may be removed based upon the capabilities of SPI. It may not be possible to remove this function as the srs function argument of ST_AsGDALRaster can be NULL, thereby instructing the function to not embed any spatial reference information into the output GDAL raster.

[pracine]

I'm wondering about the speed of any ST_srtext(rast) query…

Would it not be wise to add an index on the spatial_ref_sys srid column to make all those queries faster?

[robe]

Huh? There already is one. srid is the primary key of the table therefore it automatically has a unique key index associated with it.

[Bborie Park]

Adds ST_AsGDALRaster, ST_GDALDrivers and ST_srtext functions

[Bborie Park]

Attached incremental patch for adding ST_AsGDALRaster and helper functions ST_getGDALDrivers and ST_srtext. Patch can be applied with the following in the base postgis source directory.

patch -p1 < st_asgdalraster.patch

Patches for ST_Band, ST_MinMax and ST_Reclass should be applied first in the order listed.

[Bborie Park]

separate code for creating GDAL MEM raster in function rt_raster_to_gdal to separate function rt_raster_to_gdal_mem for future use

[Bborie Park]

No changes. New patch as to reflect changes in patches that this one depends on

[Bborie Park]

No changes. New patch to reflect line shifts in patches that this one depends on

[Bborie Park]

No changes made. Update diff due to incremental changes in ST_MinMax.

[Bborie Park]

Refactored for memory context changes.

[Bborie Park]

Attached latest ST_AsGDALRaster. Merges cleanly against r7145.

The following patches must be merged first for this patch to merge cleanly:

1. ST_Band

2. ST_SummaryStats

3. ST_Mean

4. ST_StdDev

5. ST_MinMax

6. ST_Histogram

7. ST_Quantile

8. ST_Reclass

[Bborie Park]

Added in r7155

[robe]

Bborie,

Does this work for you?

SELECT * FROM ST_GdalDrivers();

I get an error almost like the order of the OUT arguments is wrong. Can you confirm it works for you:

SELECT * FROM st_gdaldrivers();

Gives me error:

ERROR: invalid input syntax for integer: ""

This is running on windows.

[robe]

If I change your function to this:

CREATE OR REPLACE FUNCTION st_gdaldrivers(OUT idx text, OUT short_name text, OUT long_name text, OUT create_options text)
  RETURNS SETOF record AS
'$libdir/rtpostgis-2.0', 'RASTER_getGDALDrivers'
  LANGUAGE c IMMUTABLE STRICT
  COST 1
  ROWS 1000;

It works and gives this output (all the idx are "" or null it seems);

        VRT	Virtual Raster	
	NITF	National Imagery Transmission Format	<lots of xml>

:
:

[Bborie Park]

Replying to robe:

Bborie,

Does this work for you?

SELECT * FROM ST_GdalDrivers();

I get an error almost like the order of the OUT arguments is wrong. Can you confirm it works for you:

SELECT * FROM st_gdaldrivers();

Gives me error:

ERROR: invalid input syntax for integer: ""

This is running on windows.

Regina,

The function works fine for me in Linux. I'll do some poking at the what could be causing that error.

[Bborie Park]

Replying to robe:

If I change your function to this:

CREATE OR REPLACE FUNCTION st_gdaldrivers(OUT idx text, OUT short_name text, OUT long_name text, OUT create_options text)
  RETURNS SETOF record AS
'$libdir/rtpostgis-2.0', 'RASTER_getGDALDrivers'
  LANGUAGE c IMMUTABLE STRICT
  COST 1
  ROWS 1000;

It works and gives this output (all the idx are "" or null it seems);

        VRT	Virtual Raster	
	NITF	National Imagery Transmission Format	<lots of xml>

:
:

Now that is really strange. Everything but the lack of values for idx is correct.

[Bborie Park]

Regina,

The problem basically comes from my use of the snprintf function. I'm guessing your PostGIS windows build is on Visual Studio. I'll do what the postgis side uses once I figure out what they did.

[robe]

My PostGIS build is on MingW, but running in a VC++ PostgreSQL 9.1 beta. I'll try my 9.0 install just to rule out 9.1 beta as being the culprit.

Anyrate it is nice. I documented it (after fudging by changing my definition)

​http://www.postgis.org/documentation/manual-svn/RT_ST_GDALDrivers.html

One question I did have — for the ST_AsGDALRaster, does that always output band 1 or all bands. Was just tesitng with some samples so didn't know.

[robe]

Forgot to ask, which version of PostgreSQL are you testing with?

[Bborie Park]

I'm going to have to find a resolution for the snprintf issue ASAP as I use that function quite a bit.

ST_AsGDALRaster is a generic interface to creating an GDAL raster. Thus, it does no format-specific checking unlike ST_AsTIFF, ST_AsJPEG and ST_AsPNG. So, if you pass a multiband raster to ST_AsGDALRaster, the function will try to add all the bands to the GDAL raster.

Therefore, most users should use ST_AsTIFF, ST_AsJPEG and ST_AsPNG. Those users who don't need hand-holding (basically have read the appropriate GDAL format page) should be able to do everything that GDAL raster output supports.

Also, my development environment is on PostgreSQL 9.0.4. PostGIS output is…

"POSTGIS="2.0.0SVN" GEOS="3.2.2-CAPI-1.6.2" PROJ="Rel. 4.7.1, 23 September 2009" LIBXML="2.7.6" USE_STATS"

All have been built using GCC 4.4.4 running on Slackware 13.1.

[robe]

I checked my 9.0 and 8.4 builds and they have the same issue. You sure its snprintf at fault. It could be my GDAL driver. I'm sitll having issues with compiling hte latest GDAL 1.8 so my GDAL is an earlier 1.8. Though perhaps I should try again. It's been a couple of months since I tried last.

Which GDAL are you running with?

As far as GCC I'm still running an embarrasing gcc 3.4.5 (mingw-vista special r3) becuase ran into too many issues with the gcc 4.. versions.

[robe]

BTW: I tested the ST_AsPNG etc to coonfirm they didn't error out, but haven't output them yet to make sure I can see the output. That's next on my list after I've provided base documentation for each of the functions.

Wow you committed a lot of functions .

[Bborie Park]

My GDAL version is the latest stable version of 1.8.0.

​http://download.osgeo.org/gdal/gdal-1.8.0.tar.gz

I'm pretty sure that the problem has to do with how I use snprintf. Based upon all other uses of snprintf in postgis, my use of the function is non-standard. I use it to determine the size of any variable I'm converting to a char string to allocate exactly the amount of memory needed to store the char string. Evidently, Windows completely differently in its return value from snprintf.

I did add a very basic client in raster/scripts/python/rtgdalraster.py to test ST_AsGDALRaster. You'll want to check the source as I think I added an example in there. This is how I tested to make sure that the output GDAL rasters are valid. My final batch of tests were done using the pele image from your book and pele came out cleanly in TIFF, JPEG and PNG.

Yes, I posted a lot of functions. I think I may have waited too long to start providing patchsets. I plan on adding new features at regular intervals in the future… starting with ST_Count and then ST_Transform.

[robe]

Well on the bright-side I can compile the latest stable gdal 1.8.0. Swapping out didn't help though so you are probably right.. Anyrate in the middle of recompiling PostGIS, though I doubt that will help.

I wonder if its the same issue I have on windows as described here: #748

To overcome, I just used Mark's hack of remarking out the nls_build line in the loader/Makefile.in

Though don't see a similar setting in raster.

[Bborie Park]

Can you try r7168? I've changed all my usages of snprintf to use explicit buffer sizes instead of attempting to automatically compute the sizes on-the-fly.

I don't believe #748 is the same issue.

[robe]

Yap that works.