[gdal-dev] Setting CS for FileGDB layers

Jeremy Palmer JPalmer at linz.govt.nz
Sun Sep 30 00:52:56 PDT 2012


Hi All,

When creating an Esri FileGDB layer the spatial reference system has to be defined using the "Esri" WKT definition. Otherwise you get an "General function failure" error.

In particular the WKT coordinate system code name needs to match the name within the Esri WKT database which is complied into the FileGDB API library. However when using the EPSG coordinate definitions there are many WKT coordinate systems which do not have the matching Esri coordinate system code name - even after MorphToEsri has been called on the CS. There are some 75 coordinate systems in New Zealand, both geographic and projection, which do not morph to the required Esri WKT. Look at EPSG 2193 as an example:

EPSG WKT:

PROJCS["NZGD2000 / New Zealand Transverse Mercator 2000",
    GEOGCS["NZGD2000",
        DATUM["New_Zealand_Geodetic_Datum_2000",
            SPHEROID["GRS 1980",6378137,298.257222101,
                AUTHORITY["EPSG","7019"]],
            TOWGS84[0,0,0,0,0,0,0],
            AUTHORITY["EPSG","6167"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4167"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",173],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",1600000],
    PARAMETER["false_northing",10000000],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AXIS["Northing",NORTH],
    AXIS["Easting",EAST],
    AUTHORITY["EPSG","2193"]]

Morphed WKT:

PROJCS["NZGD2000_New_Zealand_Transverse_Mercator_2000",
    GEOGCS["GCS_NZGD2000",
        DATUM["D_NZGD_2000",
            SPHEROID["GRS_1980",6378137,298.257222101]],
        PRIMEM["Greenwich",0],
        UNIT["Degree",0.017453292519943295]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",173],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",1600000],
    PARAMETER["false_northing",10000000],
    UNIT["Meter",1]]

Esri WKT:

PROJCS["NZGD_2000_New_Zealand_Transverse_Mercator",
    GEOGCS["GCS_NZGD_2000",
        DATUM["D_NZGD_2000",
        SPHEROID["GRS_1980",6378137.0,298.257222101]],
        PRIMEM["Greenwich",0.0],
        UNIT["Degree",0.0174532925199433]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["False_Easting",1600000.0],
    PARAMETER["False_Northing",10000000.0],
    PARAMETER["Central_Meridian",173.0],
    PARAMETER["Scale_Factor",0.9996],
    PARAMETER["Latitude_Of_Origin",0.0],
    UNIT["Meter",1.0],
    AUTHORITY["EPSG",2193]]

My question: is it possible to improve the morphToEsri logic to explicitly set the coordinate system code name? Having to set the coordinate system from an Esri WKT external definition is very inconvenient. I see the the Esri pcs and gcs values have already been compiled into a CSV files here: http://trac.osgeo.org/gdal/browser/trunk/autotest/osr/data

Cheers,
Jeremy

P.S here's a script that tests the nz coordinate systems (nz_cs.csv here http://dl.dropbox.com/u/30623980/nz_cs.csv.tar.gz):

import os
import cs

from osgeo import gdal
from osgeo import ogr
from osgeo import osr

fgdb_drv = None
fgdb_path = "cs_test.gdb"
nz_cs_file = "nz_cs.csv"

try:
    fgdb_drv = ogr.GetDriverByName('FileGDB')
except:
    print( "Can't load FileGDB driver" )
    exit( 1 )

if os.path.exists( fgdb_path ):
    fgdb_drv.DeleteDataSource( fgdb_path )
fgdb_ds = fgdb_drv.CreateDataSource( fgdb_path )

failed = 0;
f = open( nz_cs_file , 'rt' )
csv_reader = csv.DictReader( f, delimiter = ',' )
for line in csv_reader:
    srs = osr.SpatialReference()
    srs.ImportFromEPSG( int( line['EPSG'] ) )
    srs.MorphToESRI()
    # set cs code to esri value, so filegdb does not fail
    srs.SetAttrValue( line['TYPE'], line['COORD_REF_SYS_CODE'] )
    lyr = fgdb_ds.CreateLayer( line['COORD_REF_SYS_CODE'] + '_TEST', srs, ogr.wkbPoint )
    if lyr is None:
        failed += 1;
f.close()

print "Number of fail %d" % failed

exit( 0 )

 

This message contains information, which is confidential and may be subject to legal privilege. If you are not the intended recipient, you must not peruse, use, disseminate, distribute or copy this message. If you have received this message in error, please notify us immediately (Phone 0800 665 463 or info at linz.govt.nz) and destroy the original message. LINZ accepts no responsibility for changes to this email, or for any attachments, after its transmission from LINZ. Thank You.


More information about the gdal-dev mailing list