Cidree / rpostgis

rpostgis: R Interface to a 'PostGIS' Database
http://cidree.github.io/rpostgis/
77 stars 14 forks source link

Lines in raster after importing from Database #19

Closed BernoGreyling closed 4 years ago

BernoGreyling commented 4 years ago

Hi,

I have a large dem stored on a remote DB and when importing with pgGetRast, lines appear on seamingly random locations :

image

When using QGIS or the readGDAL function in R, there are no issues. Any idea where this could be originating from? When importing smaller dems, there are no such artifacts.

Thanks!

dnbucklin commented 4 years ago

If you use the boundary argument in pgGetRast to subset what you're loading into R, do you still see these lines, in the same locations?

Also can you tell what raster value is being applied to these lines? It appears they are not NA(?)

BernoGreyling commented 4 years ago

Hey @dnbucklin

Here I selected 2 different boundaries. It seems as though the lines stay in the same place image

image

And then when I pull the DEM from the DB into QGIS

image

The lines have a value of 0 as shown below. The screenshot doesnt show the mouse ofcourse but it is hovering over the line: image

Thanks!

dnbucklin commented 4 years ago

Thanks @BernoGreyling ; these lines look like they could be along raster tile boundaries. If that's the case, and you pull just one of the problem tiles with from the raster with pgGetRast (use clauses to select one row), do you still see this stripe along an edge?

Also, are the raster properties (origin, resolution, extent, dimensions) in the imported R raster as expected (i.e. same as in PostGIS/QGIS)?

BernoGreyling commented 4 years ago

Hi @dnbucklin ,

It seems that almost all the tiles have these lines if I pull them row by row. And as far as I check, always at the bottom right. The rest of the properties come out perfectly, so no issues with that.

image

Thanks!

dnbucklin commented 4 years ago

Hmm, if the extent/dimensions match the tile in PostGIS, those data might actually be in the PostGIS raster. Possibly they're supposed to be nodata, and they're not getting imported correctly. You should be able to tell if that's the case with this query for one of those tiles:

select st_summary(rast), st_bandnodatavalue(rast) nodata_value, 
st_count(rast, 1, true) data_pixels, st_count(rast, 1, false) all_pixels, 
ST_ValueCount(rast,1,0) zero_value_pixels from 
{schema.table} where rid = {XX};
BernoGreyling commented 4 years ago

Running the command above I get the following:

image

dnbucklin commented 4 years ago

Based on this information, it looks like the zero values are also in the PostGIS raster, and the pgGetRast is performing as expected. It may be the case that those values are in overlap areas between tiles, but if they are not marked as nodata, they will not be ignored on import.

To get these rasters to import correctly with pgGetRast may require reclassifying those values as nodata in the PostGIS raster, and/or re-importing the raster so the extra line is not there in the first place.