Work done

This week:

• I implemented access to overviews, adding overviews tables as subdatasets of the general dataset, and adding metadata info to SUBDATASETS domain.
• I implemented raster inplace update, bt overwriting the IWriteBlock method of RasterBand.
• I reorganized some parts of the code by creating RasterWrapper and RasterBandWrapper classes. In these classes, I do all the tasks related with raster format: parse hexwkb format, swap words with needed, get and set raster attributes, etc.  The WKTRasterWrapper takes the hexwkb string fetched from database, and fill raster properties by parsing it. The most useful thing, in my opinion, is that in any moment, you can get the updated hexwkb representation of the raster again.  The WKTRasterBandWrapper allows setting new data for the band.
• I implemented a TODO task of the gdal2wktraster script: support from out-db rasters. My version of the script has been sent for review to PostGIS devel list.

The two last tasks will help me in the implementation of support for out-db rasters in the driver

Planned work for next week

• Finish support for out-db rasters.  I asked Mateusz for this option in the loader script, and he told me to implement it, if I wanted. I need this option for adding out-db support in the driver.
• Correct errors. This is a transversal task. This week I found a couple of important errors in several parts of the code.
• Update RASTER_COLUMNS table with the information fetched in IReadBlock, if needed. This was an old doubt, and I think that’s the correct point, thanks to a Frank’s comment.
• Start with block caching. In each call to IReadBlock, I fetch all the raster from database.

Problems found

• The spatial query seems to work, but I’m not sure if is the most appropiate:
  • With GIST index:

    SELECT rast  FROM table WHERE rast ~ ST_SetSRID(ST_MakeBox2D(

    ST_Point(lowerLeftX, lowerLeftY) ,ST_Point(upperRightX, upperRightY)),

    srid);

  • Without GIST index:

    SELECT rast  FROM table WHERE _ST_Contains(rast, ST_SetSRID(ST_MakeBox2D(

    ST_Point(lowerLeftX, lowerLeftY) ,ST_Point(upperRightX, upperRightY)),

    srid));

• When calculating a block  Bounding Box based on a block offset and the block size in IReadBlock, I crossed two values (X-Y), and, of course, it failed for non-square blocks. I spent a couple of days with this error
• DOUBT: What happens if the bb of a block offset in IReadBlock doesn’t match any real block? Missing-tiles raster?
• DOUBT: Where should I fit WKTRaster “16BF” datatype? In GDT_Uint16?, in a single float?

Work done

This week I finished extending the basic GDAL WKT Raster driver code, and continued working with overviews and out-db rasters.

More specifically:

• I solved the problem with “Invalid angle” changing srid of the data to 26711, the correct one.
• I used a spatial query to fetch the correct block from database.
• I created my own PQ-style array parser, to extract each element as a string.
• I allowed raster with several bands and all pixel types, taking endianess into account

After this, I started working in adding support for reading overviews. The line executed to load raster overviews:

gdal2wktraster.py -r utm.tif -t table -s 26711 -b 1 -k 100x100 -l 2 -V -I -M -v

Where table is the name of my raster table.

Just now, I’m working finishing the overviews code and working in the support for out-db rasters.

Planned work for next week

Finish the support for overviews and out-db rasters, and try to implement the raster inplace update.

Problems found

• The spatial query seems to work, but I’m not sure if is the most appropiate:
  • With GIST index:

    SELECT rast  FROM table WHERE rast ~ ST_SetSRID(ST_MakeBox2D(

    ST_Point(lowerLeftX, lowerLeftY) ,ST_Point(upperRightX, upperRightY)),

    srid);

  • Without GIST index:

    SELECT rast  FROM table WHERE _ST_Contains(rast, ST_SetSRID(ST_MakeBox2D(

    ST_Point(lowerLeftX, lowerLeftY) ,ST_Point(upperRightX, upperRightY)),

    srid));

• I spent a lot of time implementing two methods that were implemented in GDAL . These are CPLHexToBinary, in cpl_string.cpp and GDALSwapWords, from rasterio.cpp. Actually, I changed my point several times before finding these useful functions. It was a really stupid waste of time…
• In previous versions, I transformed to binary only the raster data fetched from database, instead of the complete raster representation. This caused problems detecting the end of the data with several bands. Once changed, it worked.
• Each band has its own pixel type, and I was considered the same pixel type for all the bands of a given raster. It worked because all the bands actually had the same pixel type, but I realized on the error testing the code and I changed it.

Work done

This week I’ve written the code of an almost-finished-basic-version of the read-only driver.

Why almost? Because I have to perform a spatial query, looking for all raster blocks that contains one point (the center of each block). In case of regularly-blocked rasters, I’ll get only one block. Instead of this spatial query, I have a testing one, always reading the same block of data.

And why basic? Because the driver only works for rasters with 8 bits/pixel, and one band of grey colour (colour statically returned, from now)

Apart from this, the driver seems to work. Here, the exit for command

gdalinfo -mm -stats -checksum "PG:host='localhost' dbname='gsoc09_test'
user='postgres' password='postgres' table='usa_mountain_one_band'
where='rid > 0'"

Yes, there is an error. And it’s related with the fact that the Coordinate System description doesn’t have a part with UNIT["metre",1,AUTHORITY["EPSG","9001"]]. So, gdalinfo thinks that the coordinate units are degrees, and take the UTM coordinates (the original image coordinates) as invalid angle values. I suppose that the error was to choose an invalid srid (4267) when loading the data on PostGIS, but I’m not sure yet.

Oh, I sent the midterm evaluation this week, almost forget it…

Planned work for next week

1. Fix the problem with “Invalid angle”.
2. Change the testing query for a real one.
3. Continue working to extend the basic driver to a more general one.

Problems found

Basically, some segmentation faults and  memory leaks when moving over the HEXWKB representation of the raster. Solved, from now.

Work done

I have updated the project trac page by gathering all the mails, post comments, related docs and new ideas together: http://trac.osgeo.org/gdal/wiki/WKTRasterDriver. So, I defined a new project plan, available for comments, of course.

I looked for more testing data, and I found ftp://ftp.remotesensing.org/geotiff/samples/. Useful.

Then, I continued working on the Dataset/RasterBand to achieve the Objective 1 (see project plan in trac). Still working…

Planned work for next week

Finish the basic version of the read-only driver: reading support for one-band raster, without overviews neither outdb support, from now. First objective.

Problems found

I loaded this tif file for testing: ftp://ftp.remotesensing.org/geotiff/samples/gdal_eg/cea.tif

gdal2wktraster.py -r cea.tif -t mountain_one_band -s 4267 -b 1 -k -I -O -V

And I got an assertion error:

File "/home/jorge/gsoc09/src/wktraster/scripts/gdal2wktraster.py", line 855, in wkblify_raster
 assert band_ov is not None

This error is in wkblify_raster function, when trying to get Overview 0 from RasterBand 1.The function “calculate_overviews” returns  0 overviews, but the loop is executed anyway (look at the commented part):

for nov in range(0, 1): # noverviews):

So, I decided to avoid the overview creation (not in my objective 1), load the raster without this option and investigate a bit more. My first checking was if  I was using an old version of gdal2wktraster script. Yes, I was. Time to svn update.

Once updated, I realized on some syntax changes. Basically:

• The “-O” option now doesn’t exist. A new option “-l OVERVIEW_LEVEL” is used (How do you know the limit to OVERVIEW_LEVEL?)
• New “-M” option, to execute VACUUM command against created tables. Good idea
• The “-k option” now takes an additional param, the block size desired. Then, we don’t use the block size given by GDAL (why?)
• As a consequence of previous change, the “-m blocksize” option has been deleted