[Gdal-dev] Loop through Microstation and project

[Frank Warmerdam]

On 12/14/05, Rob McCulley <RMcCulley at county24.com> wrote:
> 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:

Rob,

Hmm, it would seem to indicate that you got less features from
GetNextFeature() than indicated by GetFeatureCount().  It might be
because you haven't called ResetRead() before doing your first
GetNextFeature() call, or it might be a counting error in the driver.

I would suggest calling ResetReading() before callin the first
GetNextFeature() (before the feature loop).

If that doesn't help, then restructure you code to keep calling
GetNextFeature() till it returns None.

> 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.

This will be a bit involved.

I think you need to do something like the following
to setup a transformation object before looping through features.

    src_cs = osr.SpatialReference()
    src_cs.ImportFromProj4( 'proj=tmerc +lon_0=-115 +k_0=0.9992
+x_0=500000 +datum=NAD27' )

    dst_cs = osr.SpatialReference()
    dst_cs.ImportFromEPSG( 2152 )

    ct = osr.CoordinateTransformation( src_cs, dst_cs )

Then use the "ct" to transform the geometries.

    newGeom.transform( ct )

Note that PROJ.4 (and OGR) don't know anything spectial
about the NAD83(CSR98) datum.  They will treat EPSG:2152
effectively as a simple UTM 12 WGS84 coordinate system.
So, the source coordinate system should go through NAD27
to NAD83 grid shifting, and then the lat/long just reprojected
back into UTM 12 on that ellipsoid.

Best regards,
--
---------------------------------------+--------------------------------------
I set the clouds in motion - turn up   | Frank Warmerdam, warmerdam at pobox.com
light and sound - activate the windows | http://pobox.com/~warmerdam
and watch the world go round - Rush    | Geospatial Programmer for Rent