[postgis-users] Reference vrt datasets in postgis raster

Christian Yrrman yrrman at gmail.com
Fri Jul 11 07:33:05 PDT 2014

Hi All,

I'm a newbie to this group, so hello everyone and thanks in advance for
letting me participate :-)

After several tests and searching on the net I'm a little stuck with the
following problem:

I have downloaded a bunch of MODIS tiles (MOD13A1) from [0] that I want to
(only) reference in a postgis raster database. These tiles come in hdf
containers with several layers in it
- For a start,  I'm only interested in the EVI layer, so I created EVI vrt
files with command [1]
- This gives me several EVI*.vrt files referencing EVI with the original
sinusoidal projection. According to spatialreference.org, the closest SRS
match is sr-org6974 [2].
- As this is not in postgis, I've imported this into spatial_ref_sys [3]
- Next, I'm using the following command to insert the data into raster
table "modistest" with the following command:

for i in `find / -name "EVI_origSRS*.vrt"`; do
        raster2pgsql -s 96974 -R -F -f rast -a $i modistest;

"-s 96974" - use the freshly imported MODIS SRS from [2]
  -R  Register the raster as an out-of-db (filesystem) raster.  Provided
      raster should have absolute path to the file
  -F  Add a column with the filename of the raster.
  -f <column> Specify the name of the raster column
     -a  Appends raster into current table, must be
         exactly the same table schema.

I've built the index on rast in a separate step (not using switch -I)

This populates table modistest (using | psql...), so it looks fine. However
when I want to retrieve data from the raster, e.g. via

SELECT ST_AsTiff(rast) from public.modistest where id=1158;

I get an error message like this
ERROR:  rt_raster_from_gdal_dataset: Unable to get data from transformed
CONTEXT:  PL/pgSQL function "st_astiff" line 20 at RETURN

and in the postgresql log I get

16 days EVI' does not exist in the file system,
and is not recognised as a supported dataset name.

So it looks like the layers were not correctly referenced - however I
cannot see a mistake. I'm actually using "find / -name "EVI_origSRS*.vrt"
to get the full path for the loop (not "for i in *vrt")

Any help/advice is highly appreciated! In a next step, I would even like to
transform the sinusoidal tile afterwards with ST_transform to epsg:4326.
That would be fantastic, I'd have all tiles in original format in the file
system and let postgis do all the rest.. :-)

Thanks in advance and greetings!


[0] MODIS download website

[1] command to create vrt files:

for i in *hdf; do
        echo Processing $i
        j=`gdalinfo $i| grep SUBDATASET_2_NAME | cut -d '=' -f 2`
        echo Extracted $j
        gdal_translate -of VRT "$j" EVI_origSRS_$i.vrt


[2] SRS projeciton for MODIS

[3] SQL command to insert MODIS SRS

INSERT into spatial_ref_sys (srid, auth_name, auth_srid, proj4text, srtext)
values ( 96974, 'sr-org', 6974, '+proj=sinu +lon_0=0 +x_0=0 +y_0=0
+ellps=WGS84 +datum=WGS84 +units=m +no_defs ',
'PROJCS["MODIS Sinusoidal",
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/postgis-users/attachments/20140711/40de8e79/attachment.html>

More information about the postgis-users mailing list