[raster] Raster importation creates gaps and corrupted data

[molgor]

Hello everyone,

I have imported raster data to postgis 2.0 with raster2pgsql.

Loaded that table into postgis using the add raster postgis plugin and got an almost useless raster. Firstly I thought it could be a QGIS problem, but after quering the ST_Value of that dataset on a point where the gap occurs, it returns an empty value (not the nodata value).

I tried to import the raster with the 'out-db' option with similar results.

I exported the raster table on a png file and the gaps persists.

The png file is attached to this ticket.

The original image is 2.1GB in geotif format

[molgor]

Gaps made on raster loader raster2pgsql

[molgor]

Im sorry, where it says Loaded that table into postgis using the add raster postgis plugin It should say:

Loaded that table into QGIS using the Add a Postgis raster Layer plugin

[Bborie Park]

I'm confused as to how you loaded your raster into PostGIS. Which one of the following did you use to load that 2.1GB GeoTIFF file?

raster2pgsql

OR

QGIS raster loader plugin

If you used the command-line raster2pgsql, please provide the command line flags passed to raster2pgsql.

If you used the QGIS raster loader plugin, you should use the command-line raster2pgsql and see if it results in the same problem as the QGIS raster loader plugin isn't directly supported by the PostGIS raster devs.

Also, if you can post the output from running gdalinfo on your 2.1GB GeoTIFF file, it can help us recreate the spatial characteristics of the raster without knowing the data.

[molgor]

Yes, sure:

I used the command-line tool.

raster2pgsql -s 4326 -N -32767 -t 256x256 -I pendiente_percent_regional_250m_wgs84.tif mapas.pendiente_percnt_256 >> pendiente_perc_256.sql

And the gdalinfo exit for the loaded raster is this:

Driver: PostGISRaster/PostGIS Raster driver
Files: none associated
Size is 0, 0
Origin = (0.000000000000000,0.000000000000000)
Pixel Size = (1.000000000000000,1.000000000000000)
Subdatasets:
  SUBDATASET_1_NAME=PG:host=localhost port=5432 dbname=incendios user=cpr password=73  schema=mapas table=pendiente_percnt_256 column=rast where='rid = 1'
  SUBDATASET_1_DESC=PostGIS Raster at mapas.pendiente_percnt_256 (rast), rid = 1
  SUBDATASET_2_NAME=PG:host=localhost port=5432 dbname=incendios user=cpr password=73  schema=mapas table=pendiente_percnt_256 column=rast where='rid = 2'
  SUBDATASET_2_DESC=PostGIS Raster at mapas.pendiente_percnt_256 (rast), rid = 2
  SUBDATASET_3_NAME=PG:host=localhost port=5432 dbname=incendios user=cpr password=73 schema=mapas table=pendiente_percnt_256 column=rast where='rid = 3'
  SUBDATASET_3_DESC=PostGIS Raster at mapas.pendiente_percnt_256 (rast), rid = 3
  SUBDATASET_4_NAME=PG:host=localhost port=5432 dbname=incendios user=cpr password=73  schema=mapas table=pendiente_percnt_256 column=rast where='rid = 4'
  SUBDATASET_4_DESC=PostGIS Raster at mapas.pendiente_percnt_256 (rast), rid = 4
  SUBDATASET_5_NAME=PG:host=localhost port=5432 dbname=incendios user=cpr password=73  schema=mapas table=pendiente_percnt_256 column=rast where='rid = 5'
  SUBDATASET_5_DESC=PostGIS Raster at mapas.pendiente_percnt_256 (rast), rid = 5
  SUBDATASET_6_NAME=PG:host=localhost port=5432 dbname=incendios user=cpr password=73  schema=mapas table=pendiente_percnt_256 column=rast where='rid = 6'
  SUBDATASET_6_DESC=PostGIS Raster at mapas.pendiente_percnt_256 (rast), rid = 6
  SUBDATASET_7_NAME=PG:host=localhost port=5432 dbname=incendios user=cpr password=73  schema=mapas table=pendiente_percnt_256 column=rast where='rid = 7'
.....

 SUBDATASET_8658_NAME=PG:host=localhost port=5432 dbname=incendios user=cpr password=73  schema=mapas table=pendiente_percnt_256 column=rast where='rid = 8658'
  SUBDATASET_8658_DESC=PostGIS Raster at mapas.pendiente_percnt_256 (rast), rid = 8658
Corner Coordinates:
Upper Left  (   0.0000000,   0.0000000) 
Lower Left  (   0.0000000,   0.0000000) 
Upper Right (   0.0000000,   0.0000000) 
Lower Right (   0.0000000,   0.0000000) 
Center      (   0.0000000,   0.0000000) 

This is the gdalinfo exit for the original image:

Driver: GTiff/GeoTIFF
Files: pendiente_percent_regional_250m_wgs84.tif
Size is 28368, 19839
Coordinate System is:
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4326"]]
Origin = (-128.604911499087336,53.626968388905460)
Pixel Size = (0.002424431085498,-0.002424431085498)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (-128.6049115,  53.6269684) (128d36'17.68"W, 53d37'37.09"N)
Lower Left  (-128.6049115,   5.5286801) (128d36'17.68"W,  5d31'43.25"N)
Upper Right ( -59.8286505,  53.6269684) ( 59d49'43.14"W, 53d37'37.09"N)
Lower Right ( -59.8286505,   5.5286801) ( 59d49'43.14"W,  5d31'43.25"N)
Center      ( -94.2167810,  29.5778242) ( 94d13' 0.41"W, 29d34'40.17"N)
Band 1 Block=28368x1 Type=Float32, ColorInterp=Gray
  NoData Value=-32767

[Bborie Park]

What versions are you running for GDAL, PostgreSQL and PostGIS? Based upon your command-line flags for raster2pgsql, you shouldn't have any issues. Also, make sure your QGIS raster plugin is up-to-date.

I'll do some testing once I generate a raster of the same characteristics as your source raster.

[molgor]

gdal version is: 1.9.0

Postgres and Postgis versions are:

incendios=# select version();
                                           version                                           
 
---------------------------------------------------------------------------------------------
-
 PostgreSQL 9.1.3 on x86_64-unknown-linux-gnu, compiled by gcc (Debian 4.4.5-8) 4.4.5, 64-bit
(1 row)

incendios=# select postgis_version();
            postgis_version            
---------------------------------------
 2.0 USE_GEOS=1 USE_PROJ=1 USE_STATS=1
(1 row)

The three packages were build from source

Thank you for your quick response

[Bborie Park]

So, I created a test raster of the same characteristics as your source raster using the following.

gs -q -dNOPAUSE -sDEVICE=tiffgray -g28368x19839 -sOutputFile=temp.tif -c showpage quit

create contents of temp.tifw
--
0.002424431085498
0
0
-0.002424431085498
-128.603699283545
53.625756173363
--

gdal_translate -ot Float32 -a_srs EPSG:4326 -a_nodata -32767 temp.tif 1808.tif

I imported the raster 1808.tif into a test database using…

raster2pgsql -s 4326 -t 256x256 -I -C 1808.tif test_1808 | psql -U postgres -d test

I did not use the -N flag as your source raster already had a NODATA value as per the gdalinfo output.

Visualizing the imported raster against the source raster in QGIS shows no problems. I'll see about striping the test raster to see if there is any tearing or gaps as you're seeing.

[molgor]

I made what you proposed.

• Using the original raster it gives me same results i.e tearings and gaps.

• With the test raster, built with the same parameters as yours, it gives me no problems. It can be visualized on QGIS without the stripes

[Bborie Park]

So, I've striped the test raster using the python script attached to the file. I'm now seeing what you're seeing.

[Bborie Park]

VRT Tile from test raster found to cause problems

[Bborie Park]

TIFF version of VRT tile tile1114.vrt

[Bborie Park]

If I import tile1114.tif using raster2pgsql

raster2pgsql -s 4326 -I -C tile1114.tif tile1114 | psql -U postgres -d test

And then view the table tile1114 in QGIS, it gives me a crazy output and strange pixel values. But, when looking at the pixels through the following query…

WITH foo AS (
	SELECT ST_PixelAsPolygons(rast) AS gv FROM tile1114
)
SELECT
	(gv).val,
	count(*)
FROM foo
GROUP BY 1

There are no strange values and everything looks correct. I'll test some more tiles to see if I can reproduce this.

[molgor]

I've striped the test raster also and didn't find anything strange. It appears to be correct. I've never managed vrt files. I'll appreciate if you could post the command line to build that file (I presume it's done with gdalbuildvrt).

Also, I make a gdal_translate with the -ot Byte flag. The nodata value was set to 0 although I specified it explicitly with -a_nodata -32767. The result was a regular striped tif with no seemingly random patterns like before. Could it be something within the struc or nodata management?

[Bborie Park]

The vrt that I attached is one generated by raster2pgsql when it is processing the source raster. Normally it is only in memory (not a file) but I tweaked the raster2pgsql loader code to spit out the VRT as a file so that I could inspect it.

After doing a debug build of PostGIS and testing the loader, I can't find any problems with the entire loading process. Testing other tiles…

[Bborie Park]

After examining other tiles in the same manner described before and running a ST_SummaryStats on the table containing the original test raster

SELECT ST_SummaryStats('test_1808', 'rast');

The summary stats are as expected

(281410560,71759692800,255,0,255,255)

where the above is (count, sum, mean, std dev, min, max). This is exactly as expected of the test raster.

Are you able to share your source raster?

[molgor]

Hi dusty mugs I can share the file. How and where do I upload it?

[Bborie Park]

Good question. Can you send me an email to bkpark@… and I'll send instructions to upload the raster on one of my servers.

[Bborie Park]

fill raster with psuedo-random floating point values

[Bborie Park]

I updated the fill_1808.py script to populate the test raster with psuedo-random floating point values between 0 and 100. I also wrote a pixel-for-pixel comparison script to validate that what is loaded into PostGIS with raster2pgsql matches that of the source test raster. The comparison script requires the latest revision of postgis-trunk due to a bug I found in 2.0.

Due to the massive size of the source test raster, I had my initial comparisons only checking every N pixels (113 or some other prime number). All those tests ran fine. I'm running a final test checking every pixel.

[Bborie Park]

requires the latest revision from postgis-trunk

[Bborie Park]

From everything I have done to test for this issue, I am unable to reproduce the issue. Even the massive 28368x19839 test raster had no issues when doing pixel by pixel comparisons. As such, I suspect there is a problem outside PostGIS.

I did notice issues with QGIS's display of the raster having strange pixel values but this would have something to do with QGIS instead.

So if I don't hear back from the reporter in the next week or two, I plan on closing this ticket with "works for me". If anything, all the testing done in this ticket reinforced my confidence in the loader and PostGIS Raster code.

[Bborie Park]

Awaiting receipt of the raster that is causing this issue. Test rasters are unable to duplicate the problem.

[Bborie Park]

On the postgis-users mailing list, a different user described a similar (same?) issue regarding gaps and corrupted data.

​http://postgis.refractions.net/pipermail/postgis-users/2012-July/034710.html

The user was able to privately provide the raster in question for testing. The raster was loaded on a Linux 64-bit dev box running PostgreSQL 9.1.2 and PostGIS trunk r10035 with the the following:

raster2pgsql -s 4326 -t 10x10 -I -C sft00.tif test_1808 | psql -d test

The output of gdalinfo on sft00.tif is

Driver: GTiff/GeoTIFF
Files: sft00.tif
       sft00.tfw
       sft00.tif.aux.xml
Size is 1760, 880
Coordinate System is:
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4326"]]
Origin = (-180.102272599999992,89.946211599999998)
Pixel Size = (0.204545041250000,-0.204422967727273)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (-180.1022726,  89.9462116) (180d 6' 8.18"W, 89d56'46.36"N)
Lower Left  (-180.1022726, -89.9460000) (180d 6' 8.18"W, 89d56'45.60"S)
Upper Right ( 179.8970000,  89.9462116) (179d53'49.20"E, 89d56'46.36"N)
Lower Right ( 179.8970000, -89.9460000) (179d53'49.20"E, 89d56'45.60"S)
Center      (  -0.1026363,   0.0001058) (  0d 6' 9.49"W,  0d 0' 0.38"N)
Band 1 Block=1760x1 Type=Float32, ColorInterp=Gray
  Min=-78.197 Max=53.275 
  Minimum=-78.197, Maximum=53.275, Mean=8.062, StdDev=22.105
  Metadata:
    STATISTICS_MAXIMUM=53.275016784668
    STATISTICS_MEAN=8.0621312493133
    STATISTICS_MINIMUM=-78.196998596191
    STATISTICS_STDDEV=22.104605706292

To validate that the raster loaded into PostGIS was clean, a pixel by pixel validation was done using compare_new.sh (attached to ticket). Upon completion of the validation, the script found that every pixel value in PostGIS matched that of the source raster (as sampled using gdallocationinfo)

As the user specified, the raster DOES look corrupted in QGIS 1.7.3 when loaded with DB Manager and PostGIS Raster plugin.

The user is also running compare_new.sh (modified for PostGIS 2.0) on the user's machine to ensure that the comparison also completes cleanly. Assuming that the comparison is successful, a ticket will be filed with QGIS. QGIS 1.8 should be tested as well.

[Bborie Park]

a new pixel vlaue validation script. script used successfully on -trunk r10035

[Bborie Park]

The problem described in the second comment before this comment was due to a floating point issue in GDAL. The fix was committed to GDAL trunk in r24672 and documented in the ticket ​#4736. The raster now looks as it should in QGIS 1.7.

As for the original poster's problem, it is hoped that the poster can provide the source raster.

[robe]

sounds like nothing to do here and this may be a dupe.