[FWTools] Get projection from shp using GDAL and Python
David Fawcett
david.fawcett at gmail.com
Mon Oct 25 11:10:40 EST 2010
My guess is that you are trying to write out the spatial reference
object, not the text representation of that object.
Instead of this line:
wksht.row(row).write(8, layer.GetSpatialRef())
Try something like:
srsObj = layer.GetSpatialRef()
srsText = srsObj.ExportToWkt()
wksht.row(row).write(8, layer.srsText)
This is a little verbose, but should give you the idea.
David.
On Mon, Oct 25, 2010 at 9:38 AM, Susana Iraiis Delgado Rodriguez
<susana.delgado_s at utzmg.edu.mx> wrote:
> 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?
>
> _______________________________________________
> FWTools mailing list
> FWTools at lists.maptools.org
> http://lists.maptools.org/mailman/listinfo/fwtools
> http://fwtools.maptools.org/
>
>
More information about the FWTools
mailing list