Wrong results from the Rasterlite2 SQL function RL2_GetPixelFromRasterByPoint

[emiliano gregori]

Hi everyone, greetings from Italy.
My first post here.

Sandro, I'm getting wrong results from the Rasterlite2 SQL function RL2_GetPixelFromRasterByPoint().

In my tests I first imported a single DEM raster in a single-coverage single-section database with rl2tool, then compared the results of this kind of query:

ATTACH DATABASE 'dem.db' AS 'dem';

SELECT RL2_GetPixelValue( RL2_GetPixelFromRasterByPoint('dem','demcoverage', MakePoint(${x}, ${y}, ${srid}), 0), 0);

(where ${} is Bash syntax for parameter substitution) with the results of GDAL gdallocationinfo on the original DEM file, at the same point (points given in the same SRID of the raster).

The results differed.

I tested several GeoTIFF DEMs from TINITALY and ASTER: files from TINITALY are single-band float rasters with SRID 32632, while ASTER ones are int16 with SRID 4326, carefully checking for missing data.

I also converted some files to Esri ASCII Grid format, to import in the coverage database from it instead of GeoTIFF, and nothing changed.

In the end I crafted several GeoTIFFs with known integer pixel values.

Also tested several combinations of different Spatialite and Rasterlite2 releases, both with system packages and Rasterlite2 compiled from sources using the latest commit.

Long story short, the SQL query results were systematically different from gdallocationinfo ones.

So I started debugging the code, and I may have found several bugs in a few functions called by rl2_pixel_from_raster_by_point() that will affect the results.

1. Wrong pointer arithmetic in do_update_pixel8() and similar functions.

Pixel address is calculated as:

in += (rst->width * y * rst->nBands) + (rst->nBands + x);

which instead should be:

in += (rst->width * y * rst->nBands) + (rst->nBands * x);

This is clearly a typo.

The same correction applies to the other functions do_update_pixel16(), do_update_pixel32(), do_update_pixelflt() and do_update_pixeldbl().

2. Wrong calculation of resolution in rl2_raster_georeference_frame()

This function is called by rl2_load_cached_raster() to set georeference parameters in the decoded raster of a tile. 

Tile's minimum and maximum x and y coordinates are passed as parameters, while horizontal and vertical resolution are calculated as:

    rst->minX = min_x;

    rst->minY = min_y;
    rst->maxX = max_x;
    rst->maxY = max_y;
    hExt = max_x - min_x;
    vExt = max_y - min_y;
    horz_res = hExt / (double) (rst->width);
    vert_res = vExt / (double) (rst->height);
    rst->hResolution = horz_res;
    rst->vResolution = vert_res;

This does not work when tiles are only partially filled, for example the rightmost or lowermost tiles in sections whose width or height is not an integer multiple of tile size.

In that case "partially populated Tiles will have a reduced BBOX simply representing the portion being effectively populated".

So when hExt or vExt are calculated from a reduced BBOX they will underestimate the corresponding horz_res or vert_res values.

This will affect the following calculation of dx or dy in rl2_pixel_from_raster_by_point(), and then the calculation of pixel address in do_update_pixel8() and similar functions.

How to solve this?

It seems that rl2_raster_georeference_frame() is called only by rl2_load_cached_raster() and fnct_ImportSectionRawPixels().

So maybe the simplest way would be to pass horizontal and vertical resolutions as additional parameters of rl2_raster_georeference_frame() ?

Inside rl2_load_cached_raster() we could query the resolutions from the tables: first inspect table raster_coverages to know whether the policy is strictResolution or not, then query the resolution at the right pyramidal level from the appropriate table "xxx_levels" or "xxx_section_levels".

This should be doable.

Instead I have some doubts with the fnct_ImportSectionRawPixels() function.

This is the implementation of SQL function RL2_ImportSectionRawPixels() if I understand well.

In this case, if the already existing coverage is strict-resolution, we could get the base resolution from the raster_coverages table.

But what happens for mixedResolutions? In that case the resolution(s) should be passed as additional parameter(s) to the SQL function, thus breaking the current SQL API...?

Many thanks for your outstanding work.

Best regards,

Emiliano Gregori