ST_AsGDALRaster Missing Spatial Reference

[digimap2k]

I have a table full of rasters imported through the raster2pgsql script. They all index fine, they report metadatam and SRID correctly, they work as expected with intersect statements, they transform into other SRS projections just fine. I can export them using ST_AsPNG and they look perfect.

The problem is when I export to GTiff using….

ST_AsGDALRaster(rast, 'GTiff')

I write the resulting file out to disk and run gdalinfo on it. The co-ordinates are all fine BUT the projection is written out as an unnamed LOCAL_CS rather than the expected PROJ_CS corresponding to the SRID of the raster. Is this expected behaviour?

[dustymugs]

What is the output of

SELECT postgis_full_version()

[digimap2k]

POSTGIS="2.0.1 r9979" GEOS="3.3.5-CAPI-1.7.5" PROJ="Rel. 4.8.0, 6 March 2012" GDAL="GDAL 1.9.1, released 2012/05/15" LIBXML="2.7.8" LIBJSON="UNKNOWN" RASTER

Digging a little deeper it looks like the raster constraints may be fouled up or missing in the raster columns view. I'll sort those out then repost if this fixes the issue.

[dustymugs]

Also try ST_AsTiff()

[digimap2k]

No joy, sorted out the data import and raster columns are now showing correctly but output files are still missing the projection information either through ST_AsTiff or ST_AsGDALRaster (any format).

[digimap2k]

OK so for anyone reading this looking for a solution my conclusion are the following:

It looks like the construction of output files from rasters through ST_AsPng|Tiff|GDALRaster don't include the spatial reference system. This in itself is perhaps not that suprising as most folk would most likely be constructing some sort of transformed image as the endpoint in the processing, perhaps for a map server. The information is clearly in there as gdalinfo on the entire table does correctly extract and report the PROJ_CS.

So the workaround:

First up if you don't care about the PROJ_CS then do nothing. If you want to pull rasters as GeoTiffs for input into another GIS processing stage then you need to push these into the files manually. What I do to workaround this is in Qt C++ using GDAL/OGR is…

1) Grab the files using ST_AsGDALRaster(rast,'GTiff') 2) Grab the SRID of the raster in it's EPSG integer code from the database using ST_SRID. 3) Make the WKT for the spatial reference using OGR, I guess you could also do this using the spatial ref sys table in postgres. 4) Construct a GDALDatasetH from the binary raster data. 5) Apply the WKT spatial reference. 6) Save to file to check it and gdalinfo is all as should be.

Clearly the fix needed in postGIS is to push the spatial reference directly into the file in the first place. Pretty sure it's not me and generally failing to do this as gdalinfo on the bzipped results from the postgis test runs also shows no valid PROJ_CS in any of the test files.

Hope that helps someone

[robe]

How would it include spatial ref sys ST_AsPNG. I've tried for ST_GDALASRaster and at least for 2.1, I know it includes it for at least DTEM since I exported a raster I generated in USGS DTEM and saw the spatial refsys correctly displayed in QGIS.

As far as ST_AsPNG, ST_AsJPEG, the spatial refsys is kept in a separate file, so it wouldn't be able to since it only exports the binary part (not the metadata file)

[digimap2k]

Interesting, I tried using QGIS to recreate this and when loading my raster I get a warning from QGIS requesting that I choose a spatial reference system as none was found in the raster. I did wonder if QGIS might assume WGS84 with no SRS and accept your USGS DEM without error but seems not.

Fair comment on the PNG and JPEG, they do of course only output the binary image data and not the associated geo referencing files so couldn't be expected to work. I guess to move this forward I really need to suck it up and go for postgis 2.1 in case I'm just chasing something which has been sorted. I've also got some elevation data models which I'll be trying out today in case that's relevant to SRS export.

FWIW, I streamlined my workaround a little to left join the spatial ref sys table to the raster query so it returns the WKT for the projection along with the raster which makes the fix a lot simpler in Qt…

SELECT ST_AsGDALRaster(rast, 'GTiff'), t2.srtext FROM mapping_rasters.os_roadmap t1 LEFT JOIN spatial_ref_sys t2 ON t2.srid = ST_SRID(t1.rast) WHERE foo = bar

then in C++ it's just a simple SetProjection() on the GDALDataset using the t2.srtext string and everything is happy again.

[digimap2k]

Good news, I've migrated from postgis 2.01 —> 2.03 and the problem has gone away. I did this through the stackbuilder utility on Windows x86 Postgres install.

Apologies for the distraction, happy to close this ticket.

[robe]

Great

[pracine]

My 2 cents: Normally ST_AsGDALRaster(rast, 'GTiff') and ST_AsTiff() should return the georeference embedded in the raster. If you query for ST_AsJPEG() you should also query for ST_Georeference() which is a world file.