hi, thanks to all who replied, i'm familiar with the command-line tools from gdal and ogr, just getting used to this API. plus, i'm just switching from perl to python to take advantage of the greater number of GIS libs.<br>
for future reference (and in case anyone has pointers on how to improve this code) here's a script that will read a point file translate it into long/lat projection and print out the new points and one Attirbute column ('ATTRIB' in this case).
<br><br>#!/usr/bin/python -w<br><br>import gdal<br>import osr<br>import ogr<br><br>fd = ogr.Open('treetrail.shp');<br>lyr = fd.GetLayer(0);<br>osrs = lyr.GetSpatialRef();<br><br>latlong = osr.SpatialReference();<br>latlong.ImportFromProj4
('+proj=longlat +ellps=GRS80 +datum=NAD83 +no_defs');<br>#print latlong.ExportToProj4();<br><br>transform = osr.CoordinateTransformation(osrs,latlong)<br><br>nlayers = lyr.GetFeatureCount();<br>i = 0<br><br>field_of_interest = 'ATTRIB';
<br><br>for i in range(0,nlayers - 1):<br> feat = lyr.GetNextFeature()<br><br> info = feat.GetField(field_of_interest)<br> fgeom = feat.GetGeometryRef()<br><br> x = fgeom.GetX()<br> y = fgeom.GetY()<br>
<br> # must be a way to use TransformPoints() ...<br> [tx,ty,z] = transform.TransformPoint(x,y);<br> print "%f,%f,%s" % (tx,ty,info)<br><br><br><br><br><div><span class="gmail_quote">On 3/6/06, <b class="gmail_sendername">
Matthew Perry</b> <<a href="mailto:perrygeo@gmail.com">perrygeo@gmail.com</a>> wrote:</span><blockquote class="gmail_quote" style="border-left: 1px solid rgb(204, 204, 204); margin: 0pt 0pt 0pt 0.8ex; padding-left: 1ex;">
Brent,<br><br> Once you open the dataset, you have to get the layer then the spatial<br>reference for the layer. From there you can export the spatial<br>reference to a variety of formats as needed (the "exportTo*" methods
<br>at <a href="http://gdal.maptools.org/ogr/classOGRSpatialReference.html">http://gdal.maptools.org/ogr/classOGRSpatialReference.html</a>)<br><br>>>> import ogr<br>>>> ds = ogr.Open('ca_counties.shp')<br>
>>> layer = ds.GetLayer()<br>>>> sr = layer.GetSpatialRef()<br>>>> sr<br><osr.SpatialReference instance at 0xb7b5510c><br>>>> sr.ExportToProj4()<br>'+proj=laea +lat_0=45 +lon_0=-100 +x_0=0 +y_0=0 +a=6370997 +b=6370997
<br>+units=m +no_defs '<br>>>> sr.ExportToWkt()<br>'PROJCS["US National Atlas Equal Area",GEOGCS["Unspecified datum based<br>upon the GRS 1980 Authalic<br>Sphere",DATUM["Not_specified_based_on_GRS_1980_Authalic_Sphere",SPHEROID["GRS_1980_Authalic_Sphere",6370997,0]],PRIMEM["Greenwich",0],UNIT["Degree",
0.017453292519943295]],PROJECTION["Lambert_Azimuthal_Equal_Area"],PARAMETER["latitude_of_center",45],PARAMETER["longitude_of_center",-100],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["Meter",1]]'
<br><br>- matt<br><br>On 3/6/06, Brent Pedersen <<a href="mailto:bpederse@gmail.com">bpederse@gmail.com</a>> wrote:<br>> hi, i'd like to be able to have a user upload a file (.shapefile for now). i<br>> want to then use python/gdal to get the projection information associated
<br>> with that file. where do i start?<br>><br>> i have ogr.Open(filename) and the simple stuff like that, but how do i get<br>> the projection? i am having trouble seeing how ogr fits in with stuff like<br>>
osr.SpatialReference.ImportFromESRI ...<br>><br>> thanks for any pointers,<br>> -brent<br>><br>> _______________________________________________<br>> Gdal-dev mailing list<br>> <a href="mailto:Gdal-dev@lists.maptools.org">
Gdal-dev@lists.maptools.org</a><br>> <a href="http://lists.maptools.org/mailman/listinfo/gdal-dev">http://lists.maptools.org/mailman/listinfo/gdal-dev</a><br>><br>><br><br><br>--<br>Matt Perry<br><a href="mailto:perrygeo@gmail.com">
perrygeo@gmail.com</a><br><a href="http://www.perrygeo.net">http://www.perrygeo.net</a><br></blockquote></div><br>