1

I'm attempting to read the AADT_2014.shp file from the New York State Department of Transportation average traffic data (TDV_Shapefile_AADT_2014.zip specifically through fiona (by way of geopandas). Here's the equivalent minimum working example:

import fiona
import io
import zipfile
import requests

r = requests.get("https://www.dot.ny.gov/divisions/engineering/applications/traffic-data-viewer/traffic-data-viewer-repository/TDV_Shapefile_AADT_2014.zip", 
                 verify=False)
zipfile.ZipFile(io.BytesIO(r.content), 'r').extractall()
for feature in fiona.open("AADT_2014.shp"):
    print(feature)

geopandas wraps a simple fiona call for reading this data in. It successfully works out localizing the similar NYC Lion dataset, for example.

But with this data fiona raised a KeyError:

---------------------------------------------------------------------------
KeyError                                  Traceback (most recent call last)
<ipython-input-9-d9cf07a334ea> in <module>()
      1 import geopandas as gpd
----> 2 aadt = gpd.read_file("../data/aadt/AADT_2014.shp")

/Users/Honors/anaconda/envs/nyc-traffic-model/lib/python3.5/site-packages/geopandas/io/file.py in read_file(filename, **kwargs)
     26         else:
     27             f_filt = f
---> 28         gdf = GeoDataFrame.from_features(f, crs=crs)
     29 
     30     return gdf

/Users/Honors/anaconda/envs/nyc-traffic-model/lib/python3.5/site-packages/geopandas/geodataframe.py in from_features(cls, features, crs)
    185         """
    186         rows = []
--> 187         for f in features:
    188             if hasattr(f, "__geo_interface__"):
    189                 f = f.__geo_interface__

fiona/ogrext.pyx in fiona.ogrext.Iterator.__next__ (fiona/ogrext.c:19524)()

fiona/ogrext.pyx in fiona.ogrext.FeatureBuilder.build (fiona/ogrext.c:4795)()

fiona/_geometry.pyx in fiona._geometry.GeomBuilder.build (fiona/_geometry.c:3406)()

KeyError: 2002

Why is that?

2

Unfortunately, it's not that simple. The shapefile seems corrupted with many records without geometry. You can see it using osgeo.ogr:

from osgeo import ogr
import json
reader = ogr.Open("AADT_2014.shp")
layer = reader.GetLayer()
layer.GetFeatureCount()
54955
one = layer.GetFeature(0)
json.loads(a.ExportToJson())
{u'geometry': {u'type': u'LineString', u'coordinates': [[590128.580000001, 4520399.800000003, -1.3124998304e-05], [590125.9100000001, 4520411.200000001, 0.007286875001682], [590126.2800000012, 4520415.650000002, 0.010086875001889], [590127.9700000007, 4520422.22, 0.014286875001744], [590154.9900000021, 4520488.390000004, 0.058986875001665], [590158.330000001, 4520498.2, 0.06548687500208], [590208.5700000012, 4520617.61, 0.146386875001781], [590238.97, 4520690.15, 0.195486875001734], [590267.54, 4520760.33, 0.242886875001659], [590337.6700000009, 4520928.4, 0.356586875002904], [590393.200000002, 4521062.64, 0.447386875001712], [590446.2000000002, 4521190.080000002, 0.533586875001674], [590453.9900000021, 4521208.05, 0.545786875001795], [590494.0100000007, 4521300.69, 0.608886875001645], [590513.5199999996, 4521346.230000004, 0.639786875001619], [590530.9300000016, 4521384.3, 0.665986875001636], [590563.5700000012, 4521465.859999998, 0.720886875001725], [590601.5600000005, 4521558.26, 0.78328687500175], [590605.4500000002, 4521570.96, 0.791586875001713], [590608.7599999998, 4521596.87, 0.807886875001714], [590615.83, 4521659.91, 0.847486875002005], [590618.5099999998, 4521675.04, 0.857086875002778], [590623.3800000008, 4521690.530000001, 0.867186875001607], [590638.9800000023, 4521717.930000002, 0.886886875001665], [590698.8499999996, 4521827.139999999, 0.964686875001689], [590708.7800000003, 4521843.020000003, 0.976386875001708], [590716.580000001, 4521859.440000003, 0.987786875002712], [590723.5499999998, 4521875.510000002, 0.998686875001795], [590725.5999999996, 4521887.080000002, 1.006086875001756], [590730.5000000019, 4521927.89, 1.031786875001671], [590733.6399999997, 4521960.790000003, 1.052386875001957], [590738.2100000009, 4522020.8, 1.089986875001614]]}, u'type': u'Feature', u'properties': {u'AADT_YEAR': u'2014', u'TDV_ROUTE': u'OGDEN AVE', u'OBJECTID': 1, u'AADT_TYPE': u'FORECAST', u'SPEED_TDV': u'http://ftp.dot.ny.gov/tdv/YR2014/nospeed.pdf', u'LEN': 1.09, u'EP': 1.09, u'MUNI': u'City of Bronx', u'LOC_ERROR': u'NO ERROR', u'BMP': 0.0, u'ROADWAYTYP': u'Road', u'GISID': 25594001, u'AADT': u'7847', u'EndDesc': u'UNIVERSITY AVE', u'Shape_Leng': 1744.8616779, u'BegDesc': u'JEROME AVE', u'CLASS_TDV': u'http://ftp.dot.ny.gov/tdv/YR2014/noclass.pdf', u'VOL_TDV': u'http://ftp.dot.ny.gov/tdv/YR2011/R11/01_Bronx/01_2030.pdf', u'RC_ID': u'01_2030', u'DATA_TYPE': u'VOLUME', u'CC_ID': None}, u'id': 0}
# and
end  = layer.GetFeature(54954) # last record
json.loads(end.ExportToJson())
{u'geometry': None, u'type': u'Feature', u'properties': {u'AADT_YEAR': u'2014', u'TDV_ROUTE': u'NORRIS', u'OBJECTID': 54955, u'AADT_TYPE': u'FORECAST', u'SPEED_TDV': u'http://ftp.dot.ny.gov/tdv/YR2014/nospeed.pdf', u'LEN': 0.13, u'EP': 0.13, u'MUNI': u'Town of Candor', u'LOC_ERROR': u'ROUTE NOT FOUND', u'BMP': 0.0, u'ROADWAYTYP': u'Road', u'GISID': 17361601, u'AADT': u'26', u'EndDesc': u'DEAD END', u'Shape_Leng': 0.0, u'BegDesc': u'PROSPECT VALLE', u'CLASS_TDV': u'http://ftp.dot.ny.gov/tdv/YR2014/noclass.pdf', u'VOL_TDV': u'http://ftp.dot.ny.gov/tdv/YR2009/R09/97_Tioga/97_6191.pdf', u'RC_ID': u'97_6191', u'DATA_TYPE': u'VOLUME', u'CC_ID': None}, u'id': 54954}

And if I try a Python generator (as in Fiona)

def records(file):  
    # generator 
    reader = ogr.Open(file)
    layer = reader.GetLayer()
    for i in range(layer.GetFeatureCount()):
        feature = layer.GetFeature(i)
        yield json.loads(feature.ExportToJson())
features = records("AADT_2014.shp")
feature_id = [feature['id'] for feature in features if not feature['geometry']]
len(feature_id)
19761 # there are 19761 records without geometry

The importation of the shapefile in GRASS GIS (v.in.ogr) gives the same result

Importing 54955 features (OGR layer )...
WARNING: 19761 features without geometry skipped

That does not work with Fiona

import fiona
reader = fiona.open("AADT_2014.shp")
feature_id = [feature['id'] for feature in reader if not feature['geometry']]
Traceback (most recent call last):
....
KeyError: 2002

But if I "clean" the shapefile with ogr2ogr AADT_2014_fixed.shp AADT_2014.shp (ogr2ogr every Shapefile), it works with empty geometries:

import fiona
reader = fiona.open("AADT_2014_fixed.shp")
len(reader)
54955
feature_id = [feature['id'] for feature in reader if not feature['geometry']]
len(feature_id)
19761
features = list(reader)
features[54954]
{'geometry': None, 'type': 'Feature', 'id': '54954', 'properties': OrderedDict([(u'OBJECTID', 54955), (u'BegDesc', u'PROSPECT VALLE'), (u'EndDesc', u'DEAD END'), (u'RC_ID', u'97_6191'), (u'GISID', 17361601), (u'BMP', 0.0), (u'EP', 0.13), (u'LEN', 0.13), (u'LOC_ERROR', u'ROUTE NOT FOUND'), (u'TDV_ROUTE', u'NORRIS'), (u'MUNI', u'Town of Candor'), (u'AADT_TYPE', u'FORECAST'), (u'AADT_YEAR', u'2014'), (u'DATA_TYPE', u'VOLUME'), (u'VOL_TDV', u'http://ftp.dot.ny.gov/tdv/YR2009/R09/97_Tioga/97_6191.pdf'), (u'CLASS_TDV', u'http://ftp.dot.ny.gov/tdv/YR2014/noclass.pdf'), (u'SPEED_TDV', u'http://ftp.dot.ny.gov/tdv/YR2014/nospeed.pdf'), (u'CC_ID', None), (u'ROADWAYTYP', u'Road'), (u'Shape_Leng', 0.0), (u'AADT', u'26')])}

If I use the "topological clean" result of GRASS GIS (v.out.ogr, Grass eliminate all the redundant geometries and the records without geometries)

import fiona
reader = fiona.open("grass_aadt.shp")
len(reader)
35509
feature_id = [feature['id'] for feature in reader if not feature['geometry']]
len(feature_id)
0

Thus the original shapefiles can be read

  • with QGIS, OpenJUMP (signalling errors), GRASS GIS (correcting errors) but not gvSIG for example
  • with ogrinfo -al AADT_2014.shp (list of features)
  • with python and osgeo.ogr
  • with pyshp (shapefile) (pure Python module), with errors with records without geometry
 import shapefile
 reader = shapefile.Reader("AADT_2014.shp")
 features = list(reader.iterShapes())
 feature_id = [feature for feature in a if not feature.m]
 len(feature_id)
 19761
 # but
 features[0].points
 [(590128.580000001, 4520399.800000003), (590125.9100000001, 4520411.200000001), (590126.2800000012, 4520415.650000002), (590127.9700000007, 4520422.220000001), (590154.9900000021, 4520488.390000004), (590158.330000001, 4520498.199999999), (590208.5700000012, 4520617.609999999), (590238.9699999997, 4520690.15), (590267.54, 4520760.33), (590337.6700000009, 4520928.4), (590393.200000002, 4521062.640000001), (590446.2000000002, 4521190.080000002), (590453.9900000021, 4521208.050000001), (590494.0100000007, 4521300.6899999995), (590513.5199999996, 4521346.230000004), (590530.9300000016, 4521384.300000001), (590563.5700000012, 4521465.859999998), (590601.5600000005, 4521558.26), (590605.4500000002, 4521570.960000001), (590608.7599999998, 4521596.869999999), (590615.8300000001, 4521659.91), (590618.5099999998, 4521675.039999999), (590623.3800000008, 4521690.530000001), (590638.9800000023, 4521717.930000002), (590698.8499999996, 4521827.139999999), (590708.7800000003, 4521843.020000003), (590716.580000001, 4521859.440000003), (590723.5499999998, 4521875.510000002), (590725.5999999996, 4521887.080000002), (590730.5000000019, 4521927.890000001), (590733.6399999997, 4521960.790000003), (590738.2100000009, 4522020.800000001)]
 features[0].__geo_interface__  
{'type': 'LineString', 'coordinates': ((590128.580000001, 4520399.800000003), (590125.9100000001, 4520411.200000001), (590126.2800000012, 4520415.650000002), (590127.9700000007, 4520422.220000001), (590154.9900000021, 4520488.390000004), (590158.330000001, 4520498.199999999), (590208.5700000012, 4520617.609999999), (590238.9699999997, 4520690.15), (590267.54, 4520760.33), (590337.6700000009, 4520928.4), (590393.200000002, 4521062.640000001), (590446.2000000002, 4521190.080000002), (590453.9900000021, 4521208.050000001), (590494.0100000007, 4521300.6899999995), (590513.5199999996, 4521346.230000004), (590530.9300000016, 4521384.300000001), (590563.5700000012, 4521465.859999998), (590601.5600000005, 4521558.26), (590605.4500000002, 4521570.960000001), (590608.7599999998, 4521596.869999999), (590615.8300000001, 4521659.91), (590618.5099999998, 4521675.039999999), (590623.3800000008, 4521690.530000001), (590638.9800000023, 4521717.930000002), (590698.8499999996, 4521827.139999999), (590708.7800000003, 4521843.020000003), (590716.580000001, 4521859.440000003), (590723.5499999998, 4521875.510000002), (590725.5999999996, 4521887.080000002), (590730.5000000019, 4521927.890000001), (590733.6399999997, 4521960.790000003), (590738.2100000009, 4522020.800000001))}
# and
features[54954].points
[]
features[54954].__geo_interface__  
Traceback (most recent call last):
...
AttributeError: _Shape instance has no attribute 'parts'
1

This appears to perhaps be a core library issue, as the file opens without incident on QGIS. I've raised a GitHub issue for the matter.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.