Ticket #1839 (closed defect: fixed)

Opened 2 years ago

Last modified 2 years ago

[raster] raster2pgsql is unable to handle TIFF subdatasets

Reported by: turbothomas Owned by: dustymugs
Priority: medium Milestone: PostGIS 2.0.1
Component: raster Version: 2.0.x
Keywords: raster2pgsql, TIFF, subdataset, history Cc: rouault

Description

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

Attachments

raster2pgsqlx64_r9822.zip Download (442.3 KB) - added by robe 2 years ago.
2.0.1SVN raster2pgsql 64-bit PostgreSQL windows

Change History

  Changed 2 years ago by dustymugs

  • owner changed from mcayland to dustymugs
  • component changed from loader/dumper to raster
  • summary changed from raster2pgsql is unable to handle TIFF subdatasets to [raster] raster2pgsql is unable to handle TIFF subdatasets

follow-up: ↓ 3   Changed 2 years ago by dustymugs

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

in reply to: ↑ 2   Changed 2 years ago by 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

  Changed 2 years ago by dustymugs

  • status changed from new to assigned

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.

follow-up: ↓ 10   Changed 2 years ago by 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.

  Changed 2 years ago by dustymugs

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

  Changed 2 years ago by dustymugs

  • keywords subdataset, history added; subdataset removed

Fixed for 2.1 in r9809

  Changed 2 years ago by rouault

  • cc rouault added

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.)

  Changed 2 years ago by dustymugs

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

in reply to: ↑ 5   Changed 2 years ago by 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.

  Changed 2 years ago by 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...

  Changed 2 years ago by 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

  Changed 2 years ago by turbothomas

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

follow-ups: ↓ 15 ↓ 16   Changed 2 years ago by 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.

Changed 2 years ago by robe

2.0.1SVN raster2pgsql 64-bit PostgreSQL windows

in reply to: ↑ 14   Changed 2 years ago by turbothomas

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

in reply to: ↑ 14   Changed 2 years ago by 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

follow-up: ↓ 18   Changed 2 years ago by 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.

in reply to: ↑ 17   Changed 2 years ago by 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)

follow-up: ↓ 20   Changed 2 years ago by 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.

in reply to: ↑ 19   Changed 2 years ago by 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 '.

follow-up: ↓ 22   Changed 2 years ago by 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

in reply to: ↑ 21   Changed 2 years ago by 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 :(

  Changed 2 years ago by dustymugs

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

follow-up: ↓ 25   Changed 2 years ago by 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/) .

in reply to: ↑ 24   Changed 2 years ago by 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

  Changed 2 years ago by dustymugs

  • status changed from assigned to closed
  • resolution set to fixed

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

  Changed 2 years ago by pracine

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

  Changed 2 years ago by dustymugs

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

Note: See TracTickets for help on using tickets.