[gdal-dev] [WKTRaster] WKT Raster band to GDAL band

[Jorge Arevalo]

Hello,

I've created a WKTRasterBand-to-GDALRasterBand conversor, in the frame of
WKT Raster project. I'm not sure if the results are correct. You can get the
code from here
https://svn.osgeo.org/postgis/spike/wktraster/rt_core/rt_api.c, but
basically, the steps are:

1.- Register GDAL/OGR drivers

2.- Create empty OGR in-memory
Datasource: OGR_Dr_CreateDataSource(OGRGetDriverByName("Memory"),"", NULL);

3.- Create empty GDAL Dataset, taking width and height from WKT Raster Band
("band"
variable): GDALCreate(GDALGetDriverByName("MEM"), "",rt_band_get_width(ctx,
band), rt_band_get_height(ctx, band), 0, GDT_Byte, NULL);

4.- Add default geotransform to GDAL empty dataset: {0.0, 1.0, 0.0, 0.0,
0.0,1.0 }; (I can't get geotransform from WKT Raster Band)

5.- Add GDAL band, getting data from WKT Raster Band:

 sprintf(szDataPointer, "DATAPOINTER=%p", rt_band_get_data(ctx, band));
 apszOptions[0] = szDataPointer;
 GDALAddBand(memdataset, nPixelType, apszOptions)

(pixel type got from WKT Raster Band too)

6.- Create empty OGR layer, to store polygonized
raster: OGR_DS_CreateLayer(memdatasource, "Polygonized layer", NULL,
        wkbPolygon, NULL);

7.- Get GDAL RasterBand from GDAL Dataset, and polygonize it

 gdal_band = GDALGetRasterBand(memdataset, 1);
 GDALPolygonize(gdal_band, NULL, hLayer, -1, NULL, NULL, NULL);

8.- Get all the features from layer and export them to WKT format:

 for(j = 0; j < nFeatureCount; j++)
 {
        hFeature = OGR_L_GetFeature(hLayer, j);
        hGeom = OGR_F_GetGeometryRef(hFeature);
        OGR_G_ExportToWkt(hGeom, &pszSrcText);
 }

So, at the end, pszSrcText has each of the WKT polygons from original
raster. My problems are:

1.- I'm testing this code with really simple raster, and I'm confused. For
example:
* 2x2 32 bits signed int raster, with "1" in all cells. The result:
POLYGON((0 0, 0 2, 2 2, 2 0, 0 0)). I was expecting POLYGON((0 0, 0 1, 1 1,
1 0, 0 0))
* the raster from
http://trac.osgeo.org/postgis/attachment/wiki/WKTRaster/SpecificationWorking01/WKTRasterEnvelopeConvexHullAndShape.gif.
The result:
POLYGON ((3 1,3 2,2 2,2 3,1 3,1 6,2 6,2 7,3 7,3 8,5 8,5 6,3 6,3 3,4 3,5 3,5
1,3 1))
POLYGON ((3 3,3 6,6 6,6 3,3 3))
POLYGON ((5 1,5 3,6 3,6 6,5 6,5 8,6 8,6 7,7 7,7 6,8 6,8 3,7 3,7 2,6 2,6 1,5
1))
POLYGON ((0 0,0 9,9 9,9 0,0 0),(6 7,6 8,3 8,3 7,2 7,2 6,1 6,1 3,2 3,2 2,3
2,3 1,6 1,6 2,7 2,7 3,8 3,8 6,7 6,7 7,6 7))

Something strange here? Am I doing things right?

2.- I'd like to store the pixel value in each polygon. Is a good idea to
create a field in the layer to store pixel value?

3.- From a GDAL RasterBand, I can fetch the owning dataset handle
(GDALGetBandDataset).I think it could be useful to get the WKT Raster handle
from a WKT RasterBand too... Am I right?

Sorry, I know this post is really long. Thanks for reading it, anyway :-)

Best regards,
Jorge

http://www.gis4free.org/blog
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://lists.osgeo.org/pipermail/gdal-dev/attachments/20100306/52e73fde/attachment.html