Read an ArcGIS Feature dataset using GDAL

Question

I am new to GDAL in Python and would like to use it to create a list of populated layers in a ArcGIS mdb Feature datasets. I have the code below that provides a list of all the populated layers but it returns a list for all the datasets in the database. The database has Feature datasets for municipal water systems, sewer systems, transportation systems, etc. How can I return only the layer counts for one of these Feature datasets and all of its layers; i.e Feaure dataset = water_systems and returns layers of Water_lines, Fire_Hydrants, Water_Valves, etc.?

Thanks in advance.

import osgeo.ogr as ogr
driver = ogr.GetDriverByName('PGeo')
ds = driver.Open(database)
numLayers = ds.GetLayerCount()
for layerNum in range(numLayers):
    layer = ds.GetLayer(layerNum)
    if layer.GetFeatureCount() > 0:
        print '%s features: %s' % (layer.GetName(), layer.GetFeatureCount())

[Edit: Additional Info]

When running OGRINFO -al on the data base I get a lot of information but was able to pull out the following:

OGRFeature(GDB_Items):70
    ObjectID (Integer) = 70
    UUID (String) = DB50B2B0-9DC6-4492-9CAD-4162580231B6
    Type (String) = 70737809-852C-4A03-9E22-2CECEA5B9BFA
    Name (String) = WaterV
    PhysicalName (String) = WATERV
    Path (String) = \Water_Sys\WaterV
    DatasetSubtype1 (Integer) = 1
    DatasetInfo1 (String) = SHAPE
    ...

Further down the massive list was:

OGRFeature(GDB_Items):827
    ObjectID (Integer) = 827
    UUID (String) = 22CAD8B9-DBBE-44A6-9D8C-86896CC23B72
    Type (String) = 70737809-852C-4A03-9E22-2CECEA5B9BFA
    Name (String) = Signal_3
    PhysicalName (String) = SIGNAL_3
    Path (String) = \Trans\Signal_3
    DatasetSubtype1 (Integer) = 1
    DatasetInfo1 (String) = Shape
    ...

Note the Path parameter contains the ArcGIS Feature Dataset Name!!!

How can I read this using GDAL in Python?

Answer 1

It is a little difficult to understand your problem without hearing what your code outputs.

I don't know if it is a cut and pase error or not, but I am guessing that you want to indent the last three lines of your code so that they get executed within the loop.

import osgeo.ogr as ogr
driver = ogr.GetDriverByName('PGeo')
ds = driver.Open(database)
numLayers = ds.GetLayerCount()
for layerNum in range(numLayers):
    layer = ds.GetLayer(layerNum)
    if layer.GetName() == "water_systems":
        print '%s features: %s' % (layer.GetName(), layer.GetFeatureCount())

I don't have a personal geodatabase handy to test with, but I think that this should work. It assumes that you know what the 'water systems' layer name is and in this case, I used "water_systems".

[EDIT] I think that this example should give you what you are looking for. Since you can't apparently namespace based on the feature dataset, you will need to create a list of the layers that you are interested in and then use that to determine which layers to get feature counts from.

import osgeo.ogr as ogr
database = 'test.mdb'
layerList = ['water_lines','water_valves','fire_hydrants']
driver = ogr.GetDriverByName('PGeo')
ds = driver.Open(database)
numLayers = ds.GetLayerCount()
for layerNum in range(numLayers):
    layer = ds.GetLayer(layerNum)
    if layer.GetName() in layerList:
        print '%s features: %s' % (layer.GetName(), layer.GetFeatureCount())

I think that part of the issue is that the 'personal geodatabase' format is a proprietary ESRI format and the OGR driver doesn't/can't fully access all of the metadata and interface with it as well as ESRI's proprietary tools.

[Another Edit]

OK, it isn't documented in the OGR PGeo driver page, but based on your information about the GDB_Items Layer, this should give you what you want, if you have the same feature data set name across all of the .mdb files.

import osgeo.ogr as ogr
database = 'test.mdb'
layerList = ['GDB_Items']
fdLayerList = []
goodPathString = "water_systems"
driver = ogr.GetDriverByName('PGeo')
ds = driver.Open(database)
numLayers = ds.GetLayerCount()
metaLayer = ds.GetLayerByName("GDB_Items")

feat = metaLayer.GetNextFeature()
while feat:
    fdLayerList.append(feat.path[len(goodPathString)+2:])
    feat = metaLayer.GetNextFeature()

for layerNum in range(numLayers):
    layer = ds.GetLayer(layerNum)
    if layer.GetName() in fdLayerList:
        print '%s features: %s' % (layer.GetName(), layer.GetFeatureCount())

Answer 2

As I understand it, you're saying your mdb has multiple tables (e.g., a table for "municipal water systems", a table for "sewer systems", etc), and you want to get the feature count for a specific table.

If that's the case, you can just get the feature count of the specific table (dataset/layer/etc) name that you want to look at.

if layer.GetName() == '<name of table in question>':
    count = layer.GetFeatureCount()
    table = layer.GetName()
    print 'There are %s features in %s' %(count, table)

Also, since you are new to GDAL in Python, you might find it useful to look at something like Fiona, which bills itself as "OGR's neater API". I have not used it much myself yet, but have read great things - and its documentation seems much friendlier than the default GDAL/OGR Python docs.