[FWTools] Get projection from shp using GDAL and Python

[Susana Iraiis Delgado Rodriguez]

Hello list!

I designed a python module to get all the information I need from a
shapefile. I've been using GDAL to work with my script and I already
collected all the information I need. My application saves the values in an
excel sheet, but I haven't gotten the projection info, the code is:

import os, time, socket
from xlwt import Workbook
from osgeo import ogr,gdal,osr
from dbf import *
gdal.AllRegister()
file_list = []
folders = None
 # look in this (root) folder for files with specified extension
for root, folders, files in os.walk( "C:\\" ):
 file_list.extend(os.path.join(root,fi) for fi in files if
fi.endswith(".shp"))
wrkbk = Workbook()
wksht = wrkbk.add_sheet('shp')
wksht.row(0).write(0,'ruta')
wksht.row(0).write(1,'archivo')
wksht.row(0).write(2,'x_min')
wksht.row(0).write(3,'x_max')
wksht.row(0).write(4,'y_min')
wksht.row(0).write(5,'y_max')
wksht.row(0).write(6,'geometria')
wksht.row(0).write(7,'num_elem')
wksht.row(0).write(8,'proyeccion')
wksht.row(0).write(9,'prj')
wksht.row(0).write(10,'estructura bd')
wksht.row(0).write(11,'fecha_modificacion')
wksht.row(0).write(12,'maquina_host')
wksht.row(0).write(13,'usuario')

for row, filepath in enumerate(file_list, start=1):
 #Get path and filename
 wksht.row(row).write(0, filepath)
 #Get filename
 (ruta, filename) = os.path.split(filepath)
 wksht.row(row).write(1, filename)
 # Get extent
 shapeData = ogr.Open(filepath)
 layer = shapeData.GetLayer()
 feature = layer.GetNextFeature()
 x_min, x_max, y_min, y_max = layer.GetExtent()
 x_y = layer.GetExtent()
 wksht.row(row).write(2, x_y[0])
 wksht.row(row).write(3, x_y[1])
 wksht.row(row).write(4, x_y[2])
 wksht.row(row).write(5, x_y[3])
 #Get geometry type;1=Point, 2=LineString, 3=Polygon
 defn = layer.GetLayerDefn()
 wksht.row(row).write(6, defn.GetGeomType())
 #Get featurecount()
 wksht.row(row).write(7, layer.GetFeatureCount())
 #Verify prj from shp
 n = os.path.splitext(filepath)
 p = n[0]+'.prj'
 if os.path.lexists(p):
  wksht.row(row).write(9, 1)
 else:
  wksht.row(row).write(9, 0)
 #Get host from the file
 wksht.row(row).write(12, socket.gethostname())

 #Extarct database
 t = n[0]+'_bd.txt'
 d = n[0]+'.dbf'
 if os.path.lexists(d):
  a = open (t,"w+")
  dbf = Dbf(d,new=False)

  for fldName in dbf.fieldDefs:
   a.write(fldName.name)
   a.write(" || ")
   a.write(fldName.typeCode)
   a.write("\n")
  dbf.close()
  a.close()
  wksht.row(row).write(10, t)
 else:
  print "El archivo " +n[0]+".shp" " no tiene dbf"
  wksht.row(row).write(10, "Sin bd")

#Get projection, this is the part that fails
 print layer.GetSpatialRef()
 wksht.row(row).write(8, layer.GetSpatialRef())

 #Get file last modification
 t = time.strftime("%m/%d/%Y %I:%M:%S
%p",time.localtime(os.path.getmtime(filepath)))
 wksht.row(row).write(11, t)
 #Get username
 wksht.row(row).write(13,os.environ.get("USERNAME"))

wrkbk.save('shp.xls')

SEARCH_PATH = os.getcwd()
TARGET_FILE = os.path.realpath('shp.xls')
print "Buscando en", SEARCH_PATH, "and writing to", TARGET_FILE
print "Encontrando archivos..."
print "Escribiendo archivo de Excel..."
print "Listo."

I got the next output:
Python 2.6.6 (r266:84297, Aug 24 2010, 18:46:32) [MSC v.1500 32 bit (Intel)]
on
win32
Type "help", "copyright", "credits" or "license" for more information.
>>> import crawler_shp
None
None
PROJCS["WGS_1984_UTM_Zone_13N",
    GEOGCS["GCS_WGS_1984",
        DATUM["WGS_1984",
            SPHEROID["WGS_1984",6378137,298.257223563]],
        PRIMEM["Greenwich",0],
        UNIT["Degree",0.017453292519943295]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",-105],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["Meter",1]]
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "crawler_shp.py", line 85, in <module>
    wksht.row(row).write(8, layer.GetSpatialRef())
  File "C:\Python26\lib\site-packages\xlwt\Row.py", line 248, in write
    raise Exception("Unexpected data type %r" % type(label))
Exception: Unexpected data type <class 'osgeo.osr.SpatialReference'>
>>>
Any idea?
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://lists.maptools.org/pipermail/fwtools/attachments/20101025/4a24053e/attachment.htm