[Gdal-dev] Loop through Microstation and project

[Rob McCulley]

Hello All,

I'm fairly new to GDAL/OGR, and Python, so forgive my ignorance.

I'm trying to loop through a bunch of microstation files, and copy the features into a single shapefile.  Everything works fine at first, but at the end of the first microstation file I get the following error:

Traceback (most recent call last):
  File "S:\Rob\Python_Scripts\ProcessAltalis.py", line 26, in ?
    oldGeom = feature.GetGeometryRef()
AttributeError: 'NoneType' object has no attribute 'GetGeometryRef'

I assume it is a result of how I iterate through the microstation features.  Here is my code:

import os
import glob
from gdal import ogr
list = glob.glob('C:/GIS/Altalis/20051121/Cadastral/*.dgn')
driver = ogr.GetDriverByName('ESRI Shapefile')
cadastral_shp = driver.CreateDataSource('C:/GIS/Altalis/20051121/Cadastral/Cadastral.shp')
cadastral = cadastral_shp.CreateLayer('Cadastral',geom_type=ogr.wkbLineString)
field = ogr.FieldDefn('Level',ogr.OFTInteger)
field.SetWidth(8)
cadastral.CreateField(field)
field = ogr.FieldDefn('File',ogr.OFTString)
field.SetPrecision(8)
cadastral.CreateField(field)
for file in list:
    microstation = ogr.Open(file)
    layer = microstation.GetLayer(0)
    layer.ResetReading()
    i = 0
    j = layer.GetFeatureCount()
    while i < j:
        feature = layer.GetNextFeature()
        oldGeom = feature.GetGeometryRef()
        geomType = oldGeom.GetGeometryType()
        if geomType == ogr.wkbLineString:
            newFeature = ogr.Feature(cadastral.GetLayerDefn())
            gml = oldGeom.ExportToGML()
            newGeom = ogr.CreateGeometryFromGML(gml)
            newFeature.SetGeometryDirectly(newGeom)
            level = feature.GetFieldAsInteger(1)
            newFeature.SetField(0, level)
            newFeature.SetField(1, os.path.basename(file))
            cadastral.CreateFeature(newFeature)
        feature.Destroy()
    microstation.Destroy()
cadastral_shp.Destroy()

The second question I have is, how do I go about converting the features from one projection to another during this process.  I know what my projections are.  The source data is in a 10TM projection (proj: +proj=tmerc +lon_0=-115 +k_0=0.9992 +x_0=500000 +datum=NAD27) and the target projection is EPSG:2152 (NAD83 CSRS98 12N).  I looked through the Projection Tutorial, but I can't figure out how to do that using python.

Thank You,
Rob McCulley
GIS Coordinator
County of Vermilion River No. 24
(780) 846-2244
rmcculley at county24.com