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.

dimanche 25 août 2013

A Linux sandbox for the benefit of GDAL/OGR binaries : seccomp_launcher

You probably already know that GDAL/OGR is a fantastic tool to read, convert and do other processing on raster and vector datasets.

But when you deal with more than 200 file formats, it is difficult (not to say impossible) to ensure that no defect exists in a code base of nearly 1 million lines (C and C++ files, empty and comment lines included), and I don't mention the sources of the load of libraries, open source or sometimes closed source, that GDAL/OGR might depend on. Especially defects that are normally not triggered by correct datasets.

If you have to process data that can come from untrusted sources, you could find yourself in a situation where an hostile party would submit a specially crafted dataset aimed at triggering a defect, with unfortunate consequences (e.g. arbitrary code execution, theft of data, ...). A page aimed at discussing security issues has been recently added on the Trac wiki to collect knowledge on that topic and provide a few recommendations (contributions from people having deployed GDAL and wishing to share the security measures that they have taken are welcome)

I have recently discovered an interesting and very elegant security mechanism provided by the Linux kernel : seccomp. The principle of that mechanism is very simple to understand : once an executable (more exactly a thread) has turned seccomp on, it can only run 4 (yes four) system calls : read(), write(), exit() and sigreturn(). System calls are the interface between a user program (e.g. a GDAL/OGR utility) and the Linux kernel. From the name, you probably figured that read() and write() are used to... read and write files (regular files, but also pipelines, network sockets). exit() is called at process termination, and sigreturn() is too obscure to be worth an explanation.

Reducing the number of system calls available to a binary considerably restricts what the user program can do, in good... or bad. In particular, once in seccomp mode ("strict seccomp", since there is also a relaxed and more customizable form of seccomp in newer Linux kernels), a program can no longer open files, create threads, initiate network connections, or even jump to an arbitrary position in a file (seek() operation) etc...

I have started recently an experiment, seccomp_launcher, to use that mechanism in order to provide a sandbox for the benefit of GDAL/OGR utilities.

Using seccomp_launcher is very simple. It is just a matter of writting "seccomp_launcher", with an optional acces mode, in front of the command. See the below examples :

$ ./seccomp_launcher gdalinfo some.tif -stats

$ ./seccomp_launcher -rw gdal_translate some.tif target.tif

$ ./seccomp_launcher -rw ogr2ogr -f filegdb out.gdb poly.gdb -progress

$ ./seccomp_launcher python swig/python/samples/gdalinfo.py some.tif
And now, a situation where it can prevent a confidential file (my private SSH key) from being accessed :

 $ ./seccomp_launcher gdal_translate hostile.vrt out.tif

INFO: in PR_SET_SECCOMP mode
AccCtrl: open(/home/even/.ssh/id_dsa,2,00) rejected. Not in white list
AccCtrl: open(/home/even/.ssh/id_dsa,0,00) rejected. Not in white list
ERROR 4: Unable to open /home/even/.ssh/id_dsa.
Permission denied
GDALOpen failed - 4
Unable to open /home/even/.ssh/id_dsa.
Permission denied

The software is made of two main parts :
  • the seccomp_launcher binary (source: seccomp_launcher.c), which as implied by its name, launches the user binary, and can run priviledge system calls (opening a file, etc...) on its behalf, after having checked that they are authorized.
  • the libseccomp_preload.so dynamic library (source: seccomp_preload.c) that is "injected" into the user binary (e.g. gdalinfo) before it starts, to force it to run in seccomp mode, and forward priviledged system calls to seccomp_launcher, by overriding some interesting entry points of the GNU libc library.
You will find more technical details in the README file.

How to build it ? (provided that you have a C compiler, gcc or clang, and make already installed)

 
What else to mention ?
  • This is still alpha / experimental code made available under the "release early, release often" motto. In particular, no independant audit of the seccomp_launcher.c code has been made yet, which is the critical place where the security checks and delegated system calls are done. So if you intend to use it in production, please take some time to review it. Please also take time to read the README carefully, in particular the intended scope of the software (in short: do not use it to protect against hostile binaries, but only against hostile input data)
  • It is available under the same X/MIT licence as the GDAL/OGR sources.
  • Contributions (testing, bug reports, code contributions) are of course welcome !


dimanche 30 décembre 2012

A geocoding client API in GDAL/OGR

A new API has been added to GDAL/OGR to use geocoding services, such as OpenStreetMap Nominatim, MapQuest Nominatim, Yahoo! PlaceFinder, GeoNames or Bing.

From C, this is as simple as :
/* Create a session with default options */
OGRGeocodingSessionH hSession = OGRGeocodeCreateSession(NULL);

/* Now query the service */
OGRLayerH hLayer = OGRGeocode(hSession, "Paris", NULL, NULL);

/* Use the result layer (the first feature is generally the most relevant) */

/* Cleanup */
OGRGeocodeFreeResult(hLayer);
OGRGeocodeDestroySession(hSession);
More interesting, you can use that capability in SQL with the ogr_geocode() function that has been added to the SQLite SQL dialect :
SELECT ST_Centroid(ogr_geocode('Paris'))

returns:

OGRFeature(SELECT):0
POINT (2.342878767069653 48.85661793020374)
or a bit more evolved, to add geocoded information to a CSV file :
ogrinfo cities.csv -dialect sqlite -sql \
   "SELECT *, ogr_geocode(city, 'country') AS country,
    ST_Centroid(ogr_geocode(city)) FROM cities" \
    --config OGR_GEOCODE_LANGUAGE en

returns:

OGRFeature(SELECT):0
id (Real) = 1
city (String) = Paris
country (String) = France
POINT (2.342878767069653 48.85661793020374)

OGRFeature(SELECT):1
id (Real) = 2
city (String) = London
country (String) = United Kingdom
POINT (-0.109369427546499 51.500506667319407)

[...]

OGRFeature(SELECT):6
id (Real) = 7
city (String) = Beijing
country (String) = People's Republic of China
POINT (116.391195 39.9064702)

For better efficiency, the results of geocoding queries are cached into a local datasource (by default, a simple SQLite database in the working directory, but the location can be overriden by a configuration option, and other formats - CSV or PostgreSQL - can be selected).

Reverse geocoding is also possible with OGRGeocodeReverse() or the ogr_geocode_reverse() SQL function :
SELECT ogr_geocode_reverse(-100,45,'country','zoom=3') AS country

returns :
OGRFeature(SELECT):0
  country (String) = United States of America
Note: the online services that you may use have generally Term of Uses that disallow bulk geocoding, so use the above capability with care.

dimanche 24 juin 2012

GDAL/OGR using Shapefile native .sbn spatial index

This is the end of a mistery and a frustration that have lasted for about 15 years. The Shapefile format had been documented since 1998, but the documentation was limited to the minimum core, that is to say the .shp file that contains the geometries and the .shx that is an index to the geometries. However the format of the .sbn file, that was known to contain spatial index (aimed at speeding up spatial filters), has never been published.

The FOSS community came with an alternate spatial index format, the .qix format, originally used by Shapelib, MapServer and GDAL/OGR, and whose use has propagated into other software stacks, such as GeoTools.

But Joel Lawhead, Marc Pfister and another developer, Francisco, have actively started since the end of 2011 the reverse engineering of the .sbn format, and their effort finally came to a conclusion. Many kudos to them for leading this great effort !

There is not yet a formal specification of the .sbn format, but with the help of the blog entries and the debug/testing code that they made available, it was relatively easy to put pieces together. So, as an application of their work and a new step to increase interoperability between FOSS and proprietary software, I've just added support for reading and using .sbn files in the latest revision of the developement version of GDAL/OGR (GDAL 2.0dev). You can refer to ticket 4719 for the details and the code.

~~~

As this code is rather new, crowd testing is of course much appreciated. Provided that you use the latest version of GDAL 2.0dev and have python-gdal ready, you can use the testsbn.py script to check if the spatial index is correctly used with your own datasets.

Let's suppose that you have a shapefile "some_shapefile.shp", with an accompaying .sbn file "some_shapefile.sbn", you can try :
python testsbn.py some_shapefile.shp
This script will iterate over all the shapes of the shapefile, and for each shape, define a spatial filter whose extent is the bounding box of the geometry and issue a spatial request. It will check that the result of the spatial request contains the shape that served to define the spatial filter.

~~~

My initial testing of .sbn vs .qix queries would tend to show that search in .sbn is faster than in .qix (more than twice on large datasets, and with a lot of spatial filter request such as in testsbn.py), while being smaller. An explaination would be that spatial requests against .qix files use floating-point comparisons, whereas equivalent requests against .sbn files only use integer comparisons.

As far as adding write support for .sbn files, I'm a bit ambivalent for now. According to Joel's blog, there are still a few details that have to be solved, in order to understand all the subtelties of the algorithm that dispatches shapes into tree nodes, and ensure full interoperability with the software stack of the vendor at the origin of the .sbn format.

vendredi 18 mai 2012

A new GDAL virtual file system to read streamed data (e.g. for OGR WFS)

GDAL/OGR can of course read data from regular file systems, but also from more exotic sources thanks to a "virtual file system" API.

Let's start with the /vsizip/ virtual file system. If you have a ZIP file, myzip.zip, that contains a shapefile myshape.shp (and its associated .shx and .dbf files), you can read with :
ogrinfo -ro /vsizip/myzip.zip/myshape.shp
(or more simply /vsizip/myzip.zip as it is considered as a directory, and a directory is a valid datasource for the Shapefile driver).

If your data is located on a HTTP/FTP server, you can use the /vsicurl/ virtual file system, like this :
ogrinfo -ro /vsicurl/htttp://example.com/myshape.shp
As a bonus, you can combine both to read inside a remote ZIP file :
ogrinfo -ro /vsizip/vsicurl/htttp://example.com/myzip.zip/myshape.shp
Developers or advanced users can be interested by the /vsimem/, to use a memory buffer as a datasource, or /vsisubfile/ to access a file located instead another one.

/vsicurl/ is convenient to access static files and provides random access to data inside them provided that the web server supports "range downloading", i.e. the capability of returning data in a range of offsets.

Unfortunately, in some circumstances, the file is dynamically generated at the time you request it, so range downloading isn't supported. One such example in the scope of GDAL/OGR is the GML document generated by a WFS GetFeature request. Currently, the OGR WFS driver fetches the document as a whole with the CPLHTTPFetch() API, and passes the buffer to the GML driver (as a /vsimem/ file).

This behaviour has at least 2 drawbacks :
  • Even if you need to read one single feature, the driver will fetch the whole WFS GetFeature response, which can be long.
  • If the WFS GetFeature response is too long, it might not fit into memory at all.
It was possible to mitigate that by using the paging capability of some WFS servers, that is a non-standardized extension for WFS 1.0 or 1.1 (now normalized in WFS 2.0 spec).

In the GDAL/OGR trunk (2.0dev >= r24460), you can find a /vsicurl_streaming/ virtual file system that can be used to read data from a streaming server. This works efficiently only if the access pattern to the data is linear, and not random access. The OGR GML driver already natively parses data as a stream, so it can work nicely with /vsicurl_streaming/ :
ogrinfo -ro -al "/vsicurl_streaming/http://testing.deegree.org/deegree-wfs/services?SERVICE=WFS&VERSION=1.1.0&REQUEST=GetFeature&TYPENAME=app:Springs"
Or, more simply, since the OGR WFS driver has been retrofitted to use it transparently :
ogrinfo -ro WFS:http://testing.deegree.org/deegree-wfs/services
app:Springs