<div dir="ltr"><div><div><div><div><div><div>Hi All,<br><br><br></div>I'm a newbie to this group, so hello everyone and thanks in advance for letting me participate :-)<br><br></div>After several tests and searching on the net I'm a little stuck with the following problem:<br>
<br></div>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<br>. <br>- For a start,  I'm only interested in the EVI layer, so I created EVI vrt files with command [1]<br>
</div><div>- This gives me several EVI*.vrt files referencing EVI with the original sinusoidal projection. According to <a href="http://spatialreference.org">spatialreference.org</a>, the closest SRS match is sr-org6974 [2]. <br>
- As this is not in postgis, I've imported this into spatial_ref_sys [3]<br></div><div>- Next, I'm using the following command to insert the data into raster table "modistest" with the following command:<br>
<br></div><div>for i in `find / -name "EVI_origSRS*.vrt"`; do<br>        raster2pgsql -s 96974 -R -F -f rast -a $i modistest;<br>done<br><br></div><div>where:<br>"-s 96974" - use the freshly imported MODIS SRS from [2]<br>
</div><div>  -R  Register the raster as an out-of-db (filesystem) raster.  Provided<br>      raster should have absolute path to the file<br>  -F  Add a column with the filename of the raster.<br>  -f <column> Specify the name of the raster column<br>
     -a  Appends raster into current table, must be<br>         exactly the same table schema.<br><br></div><div>I've built the index on rast in a separate step (not using switch -I)<br></div><div><br></div><div>This populates table modistest (using | psql...), so it looks fine. However when I want to retrieve data from the raster, e.g. via<br>
<br>SELECT ST_AsTiff(rast) from public.modistest where id=1158;<br><br></div><div>I get an error message like this<br>ERROR:  rt_raster_from_gdal_dataset: Unable to get data from transformed raster<br>CONTEXT:  PL/pgSQL function "st_astiff" line 20 at RETURN<br>
<br></div><div>and in the postgresql log I get<br><br></div><div>ERROR 4: `HDF4_EOS:EOS_GRID:"MOD13A1.A2012065.h35v10.005.2012082112322.hdf":MODIS_Grid_16DAY_500m_VI:500m 16 days EVI' does not exist in the file system,<br>
and is not recognised as a supported dataset name.<br><br><br></div><div>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")<br>
<br><br></div><div>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.. :-)<br>
<br><br><br></div><div>Thanks in advance and greetings!<br><br><br><br></div><div>Chris<br></div><div><br><br></div><div><br><br><br></div><div>[0] MODIS download website<br><a href="http://e4ftl01.cr.usgs.gov/MOLT/MOD13A1.005/">http://e4ftl01.cr.usgs.gov/MOLT/MOD13A1.005/</a><br>
</div><div><br></div><div>[1] command to create vrt files:<br></div><div>#!/bin/bash<br><br>for i in *hdf; do<br>        echo Processing $i<br>        j=`gdalinfo $i| grep SUBDATASET_2_NAME | cut -d '=' -f 2`<br>        echo Extracted $j<br>
        gdal_translate -of VRT "$j" EVI_origSRS_$i.vrt<br><br>done</div><br></div>[2] SRS projeciton for MODIS<br><a href="http://spatialreference.org/ref/sr-org/6974/">http://spatialreference.org/ref/sr-org/6974/</a><br>
<br></div>[3] SQL command to insert MODIS SRS<br><pre>INSERT into spatial_ref_sys (srid, auth_name, auth_srid, proj4text, srtext) <br>values ( 96974, 'sr-org', 6974, '+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs ', <br>
'PROJCS["MODIS Sinusoidal",<br>GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],<br>AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],<br>
UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],<br>PROJECTION["Sinusoidal"],<br>PARAMETER["false_easting",0.0],<br>
PARAMETER["false_northing",0.0],<br>PARAMETER["central_meridian",0.0],<br>PARAMETER["semi_major",6371007.181],<br>PARAMETER["semi_minor",6371007.181],<br>UNIT["m",1.0],AUTHORITY["SR-ORG","6974"]]');</pre>
<br><div><br></div></div>