[Gdal-dev] Loop through Microstation and project

[Rob McCulley]

Thanks Frank!

I added in the ResetReading() call before the loop, and then added an 'if not feature: break' check line after the call to GetNextFeature().  That solved the looping problem.

Your instructions for the coordinate transformation were perfect.  I had some trouble getting the proj declaration for 10TM to work, and then discovered that ArcView has a nice Wkt declaration of 10TM, which worked fine using ImportFromWkt().  If anyone else is stuck with some 10TM data, the Wkt for it is:

PROJCS["N83_10TM",
    GEOGCS["GCS_North_American_1983",
        DATUM["D_North_American_1983",
            SPHEROID["GRS_1980",6378137.0,298.257222101]],
        PRIMEM["Greenwich",0.0],
        UNIT["Degree",0.0174532925199433]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["False_Easting",0.0],
    PARAMETER["False_Northing",0.0],
    PARAMETER["Central_Meridian",-115.0],
    PARAMETER["Scale_Factor",0.9992],
    PARAMETER["Latitude_Of_Origin",0.0],
    UNIT["Meter",1.0]]

About six monthes ago I had written python script using ESRI's geoprocessing module.  This OGR script runs in about 20 minutes (I have 75 microstation files with ~2500 features per file).  The same script using the ESRI functions takes nearly 2 hours to run.  The best part is, OGR is much more flexible, I will be able to add a few more steps to the script that I can't add using ESRI's functions.

Kudo's to an amazing library!

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


-----Original Message-----
From: fwarmerdam at gmail.com [mailto:fwarmerdam at gmail.com]On Behalf Of
Frank Warmerdam
Sent: Wednesday, December 14, 2005 12:46 PM
To: Rob McCulley
Cc: GDAL (E-mail)
Subject: Re: [Gdal-dev] Loop through Microstation and project


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