Support for spatialite

[perrygeo]

Can we import spatialite tables directly as a vector data source? If so, document it. If not, make it so.

[perrygeo]

I think there is a bug in OGR; an off-by-one error when indexing features on Spatialite data sources. This causes the first feature of spatialite data sources (layer.GetFeature(0)) to return None which raises an Exception.

Code to reproduce it the OGR error is below. Aside from fixing the bug upstream, I don't think we should try to code around this. My advice to anyone using spatialite is to either:

1. load the data into a pandas GeoDataFrame (see examples)
2. Convert the data to another format using ogr2ogr
3. Use a custom sqlalchemy to write a generator function that reads features from spatialite and yields
   geojson-like mappings.

Until this is fixed in ogr, we'll have to resort to workarounds.

import ogr

# Two datasets with the exact same spatial data
# spatialite was created using QGIS > Save as > Spatialite
splite = "/Users/mperry/Desktop/countries.sqlite"
shp = "/Users/mperry/data/natearth/110m_cultural/ne_110m_admin_0_countries.shp"

# Open layers
splite_ds = ogr.Open(splite)
splite_lyr = splite_ds.GetLayer(0)
shp_ds = ogr.Open(shp)
shp_lyr = shp_ds.GetLayer(0)

# Confirm feature counts
assert splite_lyr.GetFeatureCount() == shp_lyr.GetFeatureCount()

# First feature, index zero
shp_feat = shp_lyr.GetFeature(0)
assert shp_feat.GetField(3) == 'Afghanistan'

# First feature with index zero returns None
splite_feat = splite_lyr.GetFeature(0)
assert splite_feat is not None  # fails

# Features appear to be 1-indexed in spatialite
splite_feat = splite_lyr.GetFeature(1)  # 'Afghanistan'

[perrygeo]

Based on this email exchange on gdal-dev, we should rewrite all OGR feature loops with GetNextFeature()

On Mon, Feb 23, 2015 at 12:46 PM, Even Rouault
<even.rouault@spatialys.com> wrote:
> Le lundi 23 février 2015 19:40:42, Matt Perry a écrit :
>> I'm running into some strange behavior with the Spatialite driver,
>> accessing features through the ogr python bindings.
>>
>>     layer.GetFeature(0) == None
>>     layer.GetFeature(1) == ...   # first feature
>>
>> for spatialite layers only; all other drivers appear to be 0-indexed.
>>
>> More details and code to reproduce:
>> https://github.com/perrygeo/python-raster-stats/issues/55#issuecomment-7560
>> 2253
>>
>> I'm using GDAL 1.11.1 on OS X via homebrew and the latest python
>> bindings from pypi. Is this a bug or undocumented behavior?
>
> Matt,
>
> I'd say driver specific behaviour. Drivers for DBMS should return a feature ID
> that is the value of the integer primary key of the table (generally starting
> at 1). For other drivers shch as shapefile etc, the convention may change
> according to the driver.
> Anyway in the general case there might be holes in numbering in case of
> deleted features. So the recommanded way of interating over features is to
> call GetNextFeature(), as GetFeature() might be slow in some drivers due to
> the lack of indexing or proper method to seek to an arbitrary feature.
>

[perrygeo]

fixed by b43e719