Changeset 2446

Timestamp:

May 14, 2014, 3:21:37 AM (10 years ago)

Author:

rouault

Message:

add_esri_column.py: add an alternative method to get ESRI datum names based on FileGDB SDK

File:

• trunk/libgeotiff/csv/add_esri_column.py (modified) (5 diffs)

trunk/libgeotiff/csv/add_esri_column.py

@@ -1 +1 @@
 #!/usr/bin/env python
+# -*- coding: utf-8 -*-
 #******************************************************************************
 #  $Id$
@@ -31 +32 @@
 #******************************************************************************
 # 
-# $Log$
-# Revision 1.2  2004/03/20 07:52:44  warmerda
-# use local paths
-#
-# Revision 1.1  2003/06/20 21:28:24  warmerda
-# New
-#
+
+# This scritps depends on a set of ArcGIS .prj files that are usually only
+# available on Frank's machine. But it can
+# also work with an alternate method if you have OGR configured with FileGDB and
+# OpenFileGDB drivers (GDAL >= 1.11)
 
 import string
@@ -46 +45 @@
 import csv_tools
 
+
+epsg_prj_path = '/usr3/data/esri/prj/epsg'
+
+def get_esri_wkt_from_gcs_code_from_epsg_prj(gcs_code):
+    filename = '%s/%d.prj' % (epsg_prj_path, gcs_code)
+
+    try:
+        esri_gcs_wkt = open(filename).read()
+        return esri_gcs_wkt
+    except:
+        print 'Failed to find ', filename
+        return None
+
+def get_esri_wkt_from_gcs_code_from_filegdb(gcs_code):
+
+    import shutil
+    try:
+        shutil.rmtree('tmp.gdb')
+    except:
+        pass
+    from osgeo import ogr
+    ds = ogr.GetDriverByName('FileGDB').CreateDataSource('tmp.gdb')
+    srs = osr.SpatialReference()
+    # Fake WKT: the important thing is the EPSG code
+    srs.ImportFromWkt( """GEOGCS["BLA",DATUM["BLA",SPHEROID["BLA",1,0]],AUTHORITY["EPSG","%d"]]""" % gcs_code )
+    lyr = ds.CreateLayer('foo', geom_type = ogr.wkbPoint, srs = srs)
+    ds = None
+
+    # Read with OpenFileGDB driver !
+    ds = ogr.Open('tmp.gdb/a00000003.gdbtable')
+    lyr = ds.GetLayer(0)
+    # Skip first WKT that is always EPSG:4326
+    feat = lyr.GetNextFeature()
+    # Here's the interesting WKT
+    feat = lyr.GetNextFeature()
+    if feat is not None:
+        ret = feat.GetField('SRTEXT')
+    else:
+        ret = None
+    ds = None
+
+    try:
+        shutil.rmtree('tmp.gdb')
+    except:
+        pass
+
+    return ret
+
 gcs_table = csv_tools.CSVTable()
 gcs_table.read_from_csv( 'gcs.csv' )
@@ -52 +99 @@
 datum_table.read_from_csv( 'datum.csv' )
 
-print '%d GCS defined.', len(gcs_table.data.keys())
+print('%d GCS defined.' % len(gcs_table.data.keys()))
 
 esri_gcs_names = {}
 esri_datum_names = {}
 
-#
-# First try looking up all the EPSG GCS codes in the ESRI "epsg" directory,
-# and using those examples to correspond an EPSG datum code with an ESRI datum
-# name ... note that datums that aren't used in the EPSG GCS set will be missed.
-#
-for gcs_code in gcs_table.data.keys():
-    gcs_code = int(gcs_code)
-    filename = '/usr3/data/esri/prj/epsg/%d.prj' % gcs_code
+try:
+    os.stat(epsg_prj_path)
+    prj_epsg_exists = True
+except:
+    prj_epsg_exists = False
+    pass
 
-    try:
-        esri_gcs_wkt = open(filename).read()
-    except:
-        print 'Failed to find ', filename
-        continue
+filegdb_method = False
+if not prj_epsg_exists:
+    from osgeo import ogr
+    if ogr.GetDriverByName('FileGDB') is not None and ogr.GetDriverByName('OpenFileGDB') is not None:
+        filegdb_method = True
 
-    srs = osr.SpatialReference()
-    srs.ImportFromWkt( esri_gcs_wkt )
+if prj_epsg_exists or filegdb_method:
+    #
+    # First try looking up all the EPSG GCS codes in the ESRI "epsg" directory,
+    # and using those examples to correspond an EPSG datum code with an ESRI datum
+    # name ... note that datums that aren't used in the EPSG GCS set will be missed.
+    #
+    for gcs_code in gcs_table.data.keys():
+        gcs_code = int(gcs_code)
 
-    gcs_name = srs.GetAttrValue( 'GEOGCS' )
-    datum_name = srs.GetAttrValue( 'DATUM' )
+        if prj_epsg_exists:
+            esri_gcs_wkt = get_esri_wkt_from_gcs_code_from_epsg_prj(gcs_code)
+        else:
+            esri_gcs_wkt = get_esri_wkt_from_gcs_code_from_filegdb(gcs_code)
+        if esri_gcs_wkt is None:
+            continue
 
-    esri_gcs_names[gcs_code] = gcs_name
+        srs = osr.SpatialReference()
+        srs.ImportFromWkt( esri_gcs_wkt )
 
-    print 'GCS %d = %s, %s' % (gcs_code, gcs_name, datum_name)
+        gcs_name = srs.GetAttrValue( 'GEOGCS' )
+        datum_name = srs.GetAttrValue( 'DATUM' )
 
-    try:
-        gcs_rec = gcs_table.get_record( gcs_code )
-        datum_code = int(gcs_rec['DATUM_CODE'])
+        esri_gcs_names[gcs_code] = gcs_name
 
-        esri_datum_names[datum_code] = datum_name
-    except:
-        print 'Failed to get gcs record, or datum info'
+        print 'GCS %d = %s, %s' % (gcs_code, gcs_name, datum_name)
+
+        try:
+            gcs_rec = gcs_table.get_record( gcs_code )
+            datum_code = int(gcs_rec['DATUM_CODE'])
+
+            esri_datum_names[datum_code] = datum_name
+        except:
+            print 'Failed to get gcs record, or datum info'
+else:
+    print 'WARNING: ESRI_DATUM_NAME column will be empty !'
 
 #
@@ -119 +182 @@
 
 datum_table.write_to_csv( 'gdal_datum.csv' )
-
-
-
-    
-
-