[raster] raster2pgsql is unable to handle TIFF subdatasets

[turbothomas]

When loading TIFFs with subdatasets, only the first subdataset is loaded, all others are ignored.

[dustymugs]

Can you provide the output from running gdalinfo on the tiff file?

[turbothomas]

Replying to dustymugs:

Can you provide the output from running gdalinfo on the tiff file?

Driver: GTiff/GeoTIFF
Files: pathto/imagename.tif
Size is 611, 1
Coordinate System is `'
Metadata:
  TIFFTAG_IMAGEDESCRIPTION=CXmlTiff type 
  TIFFTAG_RESOLUSTIONUNIT=2 (pixels/inch)
  TIFFTAG_XRESOLUTION=100
  TIFFTAG_YRESOLUTION=100
ImageStructure Metadata:
  INTERLEAVE=BAND
Subdatasets:
  SUBDATASET_1_NAME=GTIFF_DIR:1:pathto/imagename.tif
  SUBDATASET_1_DESC=Page 1 (611P x 1L x 1B)
  SUBDATASET_2_NAME=GTIFF_DIR:2:pathto/imagename.tif
  SUBDATASET_2_DESC=Page 1 (392P x 514352L x 1B)
Corner Coordinaes:
Upper Left  (   0.0, 0.0)
Lower Left  (   0.0, 1.0)
Upper Right ( 611.0, 0.0)
Lower Right ( 611.0, 1.0)
Center      ( 305.5, 0.5)
Band 1 Block=611x1 Type=Byte, ColorInterp=Gray

[dustymugs]

Should be fixed in r9807. Please close this ticket once you've tested the subdataset. I don't have any GeoTIFF files with subdatasets to test with.

[robe]

Bborie,

Do you need to make the same change to trunk or is trunk not have the issue?

turbothomas,

Sounds like you are testing on Windows maybe. Let me know if you need a build to test with. I was going to get one out this week anyway and can probably expedite it a bit or at least give you a new raster2pgsql to test now before I release.

[dustymugs]

Yes, I do need to make the same change to trunk. Just holding off on it until I can resolve #1840.

[dustymugs]

Fixed for 2.1 in r9809

[rouault]

If I'm allowed to comment a bit about the code, I don't think that the approach taken in isSubDataset() is so great in term of exaustiveness and maintainability. I'm not sure you have captured all the GDAL drivers that support subdatasets (actually I know some others, and the list can potentially grow with newer GDAL version), and it should not be the job of PostGIS to know their syntax. I'd recommand using GDALIdentify() and if it succeeds, then it's a GDAL dataset or subdataset.

Another related remark about the code around line 2537 that uses fopen() to test if it is a raster, that is based on fopen(). GDAL supports opening files on exotic filesystems like /vsizip/ , /vsicurl/ or even things that are not files, like a PostgisRaster? dataset, or a GeoRaster? dataset... Same as above, GDALIdentify() would be better.

So all in all I think you could just drop isSubDataset() and use GDALIdentify().

(Sorry if I'm perhaps misunderstanding the code. My comments are only based on code reading.)

[dustymugs]

Thanks rouault! I'll definitely look into using GDALIdentify() instead.

[turbothomas]

Replying to robe:

yes, I am working with Windows 7 64-bit. If I could get a build, I can test on my datasets.

[dustymugs]

turbothomas,

Is any of the TIFF dataset publicly available? I'd like to be able to have an example of one for future testing...

[rouault]

try  http://svn.osgeo.org/gdal/trunk/autotest/gcore/data/twoimages.tif

or use tiffcp, like in : tiffcp src.tif src.tif dst.tif

[turbothomas]

If that doesn't do I can upload one on a file server tomorrow.

[robe]

turbothomas,

Are you using the 64-bit PostgreSQL or 32-bit. Sorry been busy this week and will be much of weekend trying to meet our book deadline.

Attached is the 64-bit raster2pgsql from 2.0.1 branch. Should be a drop in replace for what you've got already. Let me know if you need the 32-bit one.

[robe]

2.0.1SVN raster2pgsql 64-bit PostgreSQL windows

[turbothomas]

I am using 64-bit PostgreSQL. Thanks, I will test it.

[turbothomas]

Only got to test it now.

The tiling seems to work fine, now. For the subsets, I am still getting the same results. Only the first one is loaded, the others are ignored. Do I have to set a special flag?

I uploaded two files I am using. One with two subsets, the other one with several

 http://www.sendspace.com/file/8u3vzq  http://www.sendspace.com/file/lzqdog

[dustymugs]

What are you passing to raster2pgsql? If I'm correct based upon with working with other raster formats (HDF4, HDF5, NetCDF) with subdatasets, you can't pass the raster file in as that automatically assumes the first subdataset. Instead, you need to pass in something like:

SUBDATASET_1_NAME=GTIFF_DIR:1:pathto/imagename.tif

The above is from the gdalinfo output you posted in a prior comment. So, the following only gets the first subdataset.

raster2pgsql -s 0 -t 50x50 pathto/imagename.tif tablename

The following gets the specific subdataset

raster2pgsql -s 0 -t 50x50 SUBDATASET_1_NAME=GTIFF_DIR:1:pathto/imagename.tif tablename

The behavior is as such because a raster can have multiple subdatasets, each of which can have multiple bands. Actually, the behavior is exactly the same as with any of the gdal utilities.

[turbothomas]

ok, I didn't get that part with the name and subdataset.

Tested it now, passing the name as dustymugs described above.

Got back

ERROR: Unable to read raster file: SUBDATASET_2_NAME=GTIFF_DIR:2:myfile.tif

(same for the first subdataset)

[dustymugs]

Can you run the following?

gdalinfo SUBDATASET_2_NAME=GTIFF_DIR:2:myfile.tif

If that doesn't work, raster2pgsql won't work either.

Also, always post your raster2psql string as there is no way for anyone helping to replicate your situation.

[turbothomas]

sorry, before I was trying

raster2pgsql SUBDATASET_2_NAME=GTIFF_DIR:2:myfile.tif  tablename

Now I tried

gdalinfo SUBDATASET_2_NAME=GTIFF_DIR:2:myfile.tif 

and got

ERROR 4: 'SUBDATASET_2_NAME=GTIFF_DIR:2:myfile.tif ' does not exist in the file system, and is not recognised as a supported dataset name.

gdalinfo failed - unable to open 'SUBDATASET_2_NAME=GTIFF_DIR:2:myfile.tif '.

[dustymugs]

Now that I look at it, you're not supposed to use the "SUBDATASET_2_NAME=" part. Only use the "GTIFF_DIR:2:myfile.tif" after the equals sign. I see that one of my prior comments let you astray. Sorry about that. So...

raster2pgsql GTIFF_DIR:2:myfile.tif  tablename

[turbothomas]

ok,

gdalinfo GTIFF_DIR:2:myfile.tif

Driver GTIFF/GeoTIFF
Files: myfile.tif
Size is 349, 77520
Coordinate System is `´
Metadata:
  TIFFTAG_DOCUMENTNAME=MICROTEC3DFAST
  TIFFTAG_RESOLUTIONUNIT=2 (pixels/inch)
  TIFFTAG_XRESOLUTION=100
  TIFFTAG_YRESOLUTION=100
Image Structure Metadata:
  COMPRESSION=LZW
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left (   0.0,    0.0)
Lower Left (   0.0,77520.0)
Upper Right( 340.0.    0.0)
Lower Right( 340.0,77520.0)
Center     ( 170.0,38760.0)
Band 1 Block=340x8 Type=Int16, ColorInterp=Gray

That works :)

Than

raster2pgsql GTIFF:2:myfile.tif test > test.sql

Processing 1/1: GTIFF_DIR:2:myfile.tif
INFO: Using default geotransform matrix (0, 1, 0, 0, 0, -1) for raster: GTIFF_DIR:2:myfile.tif
ERROR: rt_band_set_pixel_line: Coordinates out of range
ERROR: rt_band_set_pixel_line: Coordinates out of range
ERROR: rt_band_set_pixel_line: Coordinates out of range
...

That doesn't work :(

[dustymugs]

Now that is strange. I'll have to do some testing with the same rasters you've provided.

[dustymugs]

Well, that was fun. PostGIS Raster has a maximum permitted width/height size of 65535 x 65535. Unfortunately, there was no obvious error message when you attempt to exceed that limit. I was testing your Eiche*.tif file and hit it using...

raster2pgsql GTIFF_DIR:2:Eiche_1_1_2011_08_17_768_M2_3D_5mm.tif blabla 

So, the lesson learned is that I need to add an explicit error message when someone hits this maximum and you really shouldn't call it the way I did above. Please tile your raster to be at most 65535 x 65535.

And before anyone asks why the maximum width and height is set to 65535 x 65535, imagine the size of the following...

65535 width * 65535 height * 8 bytes per pixel (64BF) = 34358689800 bytes (not counting overhead)

That would be a very monstrous number considering the maximum size permitted in a single field of PostgreSQL is 1 GB ( http://www.postgresql.org/about/) .

[turbothomas]

Well, for the testing I was just to lazy to type in the tiling flag. Normally I would have used it.

But it seems that we could identify another small bug that way.

Thanks for the help

Thomas

[dustymugs]

Added appropriate error message for 2.0 in r9839. For trunk (2.1) in r 9838.

[pracine]

I would have think the message to be in raster2pgsql also. No?

[dustymugs]

Nope. The error occurs in the C API, specifically in rt_raster_new().