lundi 24 mars 2014

Draft GDAL/OGR class hierarchy for GDAL 2.0

As a result of the first day of the OSGeo Code Sprint 2014 Vienna, I just wanted to share the outcome of my thoughts for a possible re-organisation of the GDAL/OGR class hierarchy, to achieve the mythical "Grand Unification". This is really work in progress, and I'm not even sure I will stick with it tomorrow morning... But here we go...

I have identified two principal aims :
  • adding support for metadata to OGR driver, datasource and layers (validation of creation options, etc...). That one is easy : make derive the 3 base classes from the GDALMajorObject class
  • more difficult and ambitious: making it possible to have a Dataset that contains both raster and vector data. You just open the data container once and can get both raster and vector data. Possible use cases: GeoPackage, PCIDSK, Spatialite/Rasterlite, Postgis/Postgis raster, ...
And one major constraint : avoid rewriting each of the existing 211 drivers... which represent 1.2 million lines of C/C++ code, 140 000 lines of code of Python autotests...

The class hierarchy of GDAL/OGR 1.X versions is quite simple :


So definitely 2 seperate worlds.

To achieve the first aim of getting metadata into OGR, you just have to do :



And now let's consider where the second aim could lead us :






Another way of presenting it with more details is the following pseudo-code :
/* Interface for major object */
class GDALIMajorObject
{
    public:
        virtual char      **GetMetadataDomainList() = 0;
        virtual char      **GetMetadata( const char * pszDomain = "" );
        [other methods go here]
};

/* Implementation of major object, and base class for dataset, bands, layers, etc.. */
class GDALMajorObject: public GDALIMajorObject
{
     /* existing code of GDALMajorObject */
};

/* Interface for raster functions */
class GDALIRasterDataset: public GDALIMajorObject
{
    public:
        virtual int         HandleRasterData() = 0;

        virtual int         GetRasterXSize( void ) = 0;
        virtual int         GetRasterYSize( void ) = 0;
        virtual int         GetRasterCount( void ) = 0;
        virtual GDALRasterBand *GetRasterBand( int ) = 0;
        [other methods go here]
};

/* Interface vor vector functions*/
class GDALIVectorDataset: public GDALIMajorObject
{
    public:
        virtual int         HandleVectorData() = 0;
       
        virtual int         GetLayerCount() = 0;
        virtual OGRLayer    *GetLayer(int) = 0;
        [other methods go here]
};

/* Convenience interface for both raster and vector functions */
class GDALIDataset: public GDALIRasterDataset, public GDALIVectorDataset
{
    public:
        /* That's all ! */
};

/* Partial implementation of GDALIRasterDataset */
class GDALAbstractRasterDataset : public GDALIRasterDataset, public GDALMajorObject
{
    /* Current code of GDALDataset GDAL v1 goes here */
   
    public:
        virtual int         HandleRasterData() { return TRUE; }
};

/* Convenience class used by vector only drivers */
class GDALEmptyRasterDataset : GDALAbstractRasterDataset
{
    public:
        virtual int         HandleRasterData() { return FALSE; }
};

/* Partial implementation of GDALIVectorDataset */
class GDALAbstractVectorDataset : public GDALIVectorDataset, public GDALMajorObject
{
    /* Current code of OGRDatasource GDAL v1 goes here*/
   
    public:
        virtual int         HandleVectorData() { return TRUE; }
};

/* Convenience class used by raster only drivers */
class GDALEmptyVectorDataset : public GDALIVectorDataset
{
    public:
        virtual int         HandleVectorData() { return FALSE; }
       
        virtual int         GetLayerCount() { return 0; }
        virtual OGRLayer    *GetLayer(int) { return NULL; }
        [other methods go here]
};

/* Equivalent of GDALDataset GDAL v1 (plus dummy vector interface). Existing GDAL drivers would derive from it. */
class GDALRasterDataset : public GDALAbstractRasterDataset, public virtual GDALEmptyVectorDataset, public virtual GDALIDataset
{
};

/* Equivalent of OGRDatasource GDAL v1 (plus dummy raster interface). Existing OGR drivers would derive from it. */
class GDALVectorDataset : public GDALAbstractVectorDataset, public virtual GDALEmptyRasterDataset, public virtual GDALIDataset
{
};

/* GDAL v2 base class for new drivers that need both raster and vector data */
class GDALDataset : public GDALAbstractRasterDataset, public virtual GDALAbstractVectorDataset, public virtual GDALIDataset
{
};

The impact on the existing code base is :
  • Current GDAL drivers must replace mentions of GDALDataset by GDALRasterDataset (automatic conversion)
  • Current OGR drivers must replace mentions of OGRDatasource by GDALVectorDataset (automatic conversion)
  • OGRDatasourceH becomes an alias of GDALDatasetH
  • C methods cast the opaque dataset pointer GDALDatasetH to GDALIDataset before invoking the C++ methods.
 Of course, all the above is just "nice" theory (rather complicated admitedly), and I should really try to go to the practice part of it to test if it can actually work...

jeudi 16 janvier 2014

OGR OpenFileGDB driver

Last year, I blogged about the reverse-engineering of the ESRI file geodatabase format.

What's new ?

That work, thanks to funding, has lead to the writing of a OGR OpenFileGDB driver, now available in GDAL/OGR trunk repository.

Since the initial phase of the reverse-engineering, significant advances have been made and documented :
  • .gdbtablx files (that contain offsets in the .gdbtable files of features) can use optimization when the feature ids are sparsed (if a series of 1024 consecutive feature ids do not exist)
  • .gdbindexes files contain the list of fields that have an attribute index, and the filename of that index (only the format for ArcGIS v10 geodatabases for now. The v9 format is different and more complicated)
  • .atx files are such attribute index files. The OpenFileGDB driver can use those attribute indexes to speed-up simple WHERE clauses of SQL requests, or attribute filters (SetAttributeFilter())
  • .spx files are partially deciphered, but not yet to the point of being usable. They share exactly the same base structure as attribute indexes. For each feature (or group of features in non-final pages of an index of depth greater than one), the value indexed is a 8 byte structure that describe the spatial footprint. The first 4 bytes seem to be about the y coordinate, and the next 4 bytes the x coordinate. Rather simple, no ? Except that the encoding of those bytes is still not fully understood (I think I've captured the logic for point geometries at the final page level, but non-final pages or non-point geometries are still mysterious). So, for now, we use the minimum bounding rectangle found at the beginning of geometry blob to speed-up spatial queries. And during the first full scan, we build an in-memory spatial index that is used for later spatial queries in the same session.
FileGDB vs OpenFileGDB

In the current state, if we do a comparison of the OpenFileGDB driver with the FileGDB driver using FileGDB API SDK v1.3, we find :

On the plus side :
  • Can read ArcGIS 9.X Geodatabases, and not only 10 or above.
  • Can open layers with any spatial reference system.
  • Thread-safe (i.e. datasources can be processed in parallel).
  • Uses the VSI Virtual File API, enabling the user to read a Geodatabase in a ZIP file or stored on a HTTP server.
  • Faster on databases with a big number of fields.
  • Does not depend on a third-party library. Available on any platform supported by GDAL/OGR
  • Robust against corrupted Geodatabase files.

On the minus side :
  • Read-only.
  • Cannot use spatial indexes.
And now ?

In the testing process, I discovered datasets that have layers compressed using a so-called Smart Data Compression. Which is a completely different beast from standard GDB tables. Not sure how "smart" that compression is, but the end result is particularly cryptic. The only thing that can be recognized is the field names... Those .gdbtable.sdc are neither supported by the OpenFileGDB driver, and guess what?, nor by the FileGDB API. So we are on (non-)feature parity...

I've encountered a few raster File Geodatabase datasets
(apparently tiled), and a quick inspection of the tables makes me believe
that a raster driver would be doable.

That's all for now. Testers appreciated as usual ! Windows users can for example download the builds tagged "-development" kindly provided by Tamas Szekeres on gisinternals.

mercredi 9 octobre 2013

FileGDB format reverse-engineered

For those who cannot wait, you can rush directly to the resulting specification. Caution: this is work-in-progress !

Now, the introduction. FileGDB is the ESRI File Geodatabase format used natively by ArcGIS to store datasets in a file system directory.

Since the version 1.9 of GDAL/OGR, a FileGDB driver exists and provides read/update/creation support for FileGDB datasources, but it has a few limitations :
  • it relies on a free (as in beer) but closed-source library, the FileGDB API. Not ideal philosophically and practically.
  • the FileGDB API is limited to opening FileGDB datasources created by ArcGIS 10 or later, but not the ones created by earlier versions.
  • the FileGDB API has some bugs that prevents it from opening valid FileGDB layers, e.g. if the SRS is a custom projection, and users can just wait for a version from the vendor that will eventually fix them.
  • it seems to have (unnecessary) quadratic performance when the number of fields grow, e.g. with some US Census datasources that have more than 1800 fields, which makes it reading them very slow.
So I decided to open my favorite hexadecimal editor and have a closer look at what there is in the guts of the files found in a .gdb directory. From the extensions, it is obvious that the .gdbtable and .gdbtablx files should be the most interesting.

A .gdbtable matches a layer / table, and contains the description of the fields (name, type, width, etc..), geospatial information ( type of geometries, SRS, extent ) as well as the content of the rows / features. This is the equivalent of a .shp and .dbf files of shapefiles. The .gdbtablx is an index that contains the offset to each row of the .gdbtable. This is an equivalent of the .shx file of shapefiles.

To my own surprise, the process of reverse engineering went rather fast. Generating very trivial layers with the FileGDB API, with small variations, and analyzing the differences helped a lot. Most datatypes can be guessed in an obvious way with an hexadecimal editor (little-endian int16/int32/IEEE 754 float64/UTF-16 strings). A few non-obvious technical details :
  • the use of a variable length encoding for integers (which AFAICS is identical to Protocol buffer base 128 varints), mostly for coordinates (and as well to specify the length of strings). The first coordinate tuple of a geometry is encoded in "absolute" form (that must be later offseted and scaled by constants described in the geometry field definition). All following coordinates are encoded as the difference with previous coordinates.
  • understanding how datetimes are encoded took me an astonishingly long time, compared to the outcome : it is just the number of days (possibly with decimals) since 1899/30/12 00:00:00, encoded as a float64.
  • how the flags indicating the absence or presence of fields worked. The tricky part was to understand that only nullable fields are represented in the bit field.
I've developed a small Python script that dumps the content of a .gdbtable file (using its .gdbtablx companion file) and, up to now, it manages to successfully dump all the .gdbtable files I've tried, both technical tables (GDB_xxxxxx ) and user tables (vector layers), including from GDB datasources that cannot be read by the FileGDB API (old v9.X datasources, or datasources with custom projections).

Just an example on a sample layer distributed with the FileGDB API (possibly only sexy to the eyes of the fans of command line utilities) :

$ python dump_gdbtable.py /home/even/FileGDB_API/samples/data/Shapes.gdb/a00000013.gdbtable
nfeaturesx = 233
nfeatures = 233
header_offset = 40
header_length = 627
layer_geom_type = 3
polyline
nfields = 6

nbcar = 8
name = OBJECTID
nbcar_alias = 0
alias =
type = 6 (objectid)
magic1 = 4
magic2 = 2
nullable = 0

nbcar = 5
name = Shape
nbcar_alias = 0
alias =
type = 7 (geometry)
magic1 = 0
magic2 = 7
wkt = GEOGCS["GCS_North_American_1983",DATUM["D_North_American_1983",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]]
magic3 = 7
xorigin = -400.000000000000000
yorigin = -400.000000000000000
xyscale = 11258999068426.238281250000000
zorigin = -10000.000000000000000
zscale = 10000.000000000000000
morigin = -10000.000000000000000
mscale = 10000.000000000000000
xytolerance = 0.000000008983153
ztolerance = 0.001000000000000
mtolerance = 0.001000000000000
xmin = -158.090073924418533
ymin = 21.277505248718398
xmax = -67.781176946816174
ymax = 62.145206966737987
magic4 = 3
25.6301473254
79.4534567087
246.305715797
nullable = 1

nbcar = 9
name = ROUTE_NUM
nbcar_alias = 0
alias =
type = 4 (string)
width = 8
flag = 5
magic = 0
nullable = 1

nbcar = 10
name = DIST_MILES
nbcar_alias = 0
alias =
type = 2 (float32)
width = 4
flag = 5
magic = 0
nullable = 1

nbcar = 7
name = DIST_KM
nbcar_alias = 0
alias =
type = 2 (float32)
width = 4
flag = 5
magic = 0
nullable = 1

nbcar = 12
name = Shape_Length
nbcar_alias = 0
alias =
type = 3 (float64)
width = 8
flag = 3
magic = 0
nullable = 1


FID = 1
feature_offset = 671
blob_len = 17880
flags = [224]
geom_len = 17856
geom_type = 3
polyline
nb_total_points: 1530
nb_geoms: 3
minx = -118.485959167123966
miny = 29.394035111730691
maxx = -81.682130704462821
maxy = 34.087306274650871
nb_points[0] = 35
nb_points[1] = 907
[1] -118.485959167123966 34.014739115964872
[2] -118.475133048781629 34.022692037809506
[...]
[34] -118.226133325056821 34.028554888549152
[35] -118.220945200338008 34.033902601472484

[1] -118.213980138991928 34.055149589953267
[2] -118.206975896201286 34.053569604619412
[...]
[906] -95.292373751724142 29.777536210780305
[907] -95.284004512870197 29.777996158193261

[1] -95.258536827702414 29.774727145411369
[2] -95.237242253618717 29.773808093496442
[...]
[587] -81.690798069370828 30.320550259572705
[588] -81.682130704462821 30.321330304058421

Field ROUTE_NUM : "I10"
Field DIST_MILES : 2449.120117
Field DIST_KM : 3941.479980
Field Shape_Length : 40.473531

[...]


Next steps ?
  • A GDAL/OGR driver implementing this specification, without any third-party dependency, would be cool, wouldn't it ? Probably read-only in the current state of knowledge. But it would likely solve all the limitations of the existing driver, and offer the following benefits : free and open source, FileGDB v9.X compatibility, no SRS limitations. Hint : funding would be welcome to help developing it.
  • Understanding the meaning of some "magic" fields that might have a practical importance.
  • I've encountered a datasource with raster data. The dump utility would need some extra love to deal with binary fields, that I did not investigate further. But perhaps raster FileGDB datasets could be read as well.
  • For brave people : decyphering the format of the other categories of files that have been left aside : .gdbindexes, .freelist, .TablesByName.atx, .CatItemsByPhysicalName.atx, .CatItemsByType.atx, .FDO_UUID.atx, .spx. Various indexes must be in there, some potentially for creation/update operations. This might be much more difficult to guess than the .gdbtable itself, if we take into account how difficult was the reverse engineering the shapefile .sbn spatial index.