I figured it out... with all of your help. I had the coordtrans = osr.CoordinateTransformation switched. I wanted to input lat/lon (WGS84) and transform it to lcc (source) and I also had to set the source prj to lcc. I think you all said this stuff already, just took me a little while. Thanks for all the help.<br>
<br><br>import os<br>try:<br> from osgeo import ogr<br> from osgeo import osr<br><br>except:<br> import ogr<br> import osr<br><br>#~ start driver<br>driver = ogr.GetDriverByName('ESRI Shapefile')<br>#~ set workspace<br>
ws = 'C:\\Users\\deadpickle\\Desktop\\case study\\verify\\'<br>os.chdir(ws)<br><br>#~ check if it exists<br>if os.path.exists('test.shp'):<br> driver.DeleteDataSource('test.shp')<br><br>#~ create a new shp<br>
ds = driver.CreateDataSource('test.shp')<br><br>#~ create spatref for new shp<br>source = osr.SpatialReference()<br>target = osr.SpatialReference()<br>source.SetFromUserInput('+proj=lcc +lat_1=33.000000 +lat_2=45.000000 +lat_0=39.000000 +lon_0=-96.000000 +x_0=0.0 +y_0=0.0 +datum=NAD83')<br>
target.SetFromUserInput('WGS84')<br>coordtrans = osr.CoordinateTransformation(target, source)<br><br>#~ create a new data layer<br>layer = ds.CreateLayer('test',geom_type=ogr.wkbLineString,srs=source)<br>
<br>
#~ add an id field<br>fieldDefn = ogr.FieldDefn('id', ogr.OFTInteger)<br>layer.CreateField(fieldDefn)<br><br>#~ create a geometry<br>line = ogr.Geometry(ogr.wkbLineString)<br><br>#~ get the def of the layer, scheme<br>
featureDefn = layer.GetLayerDefn()<br><br># create a new feature<br>feature = ogr.Feature(featureDefn)<br><br>#~ Read 2006track file<br>file = open("2006track.csv","r")<br>file.readline()<br>lnlist = file.readlines()<br>
for lines in lnlist:<br> value = lines.split(",")<br> if float(value[1]) == 2:<br> print float(value[3]),float(value[2])<br> line.AddPoint(float(value[3]),float(value[2]))<br> #~ set the field 'ID'<br>
feature.SetField('id', float(value[1]))<br> <br>#~ trans from lat/lon input to lcc<br>line.Transform(coordtrans)<br><br>#~ set the geometry for this shp<br>feature.SetGeometry(line)<br><br># add the feature to the output layer<br>
layer.CreateFeature(feature)<br><br># destroy the geometry and feature and close the data source<br>ds.Destroy()<br>line.Destroy()<br>feature.Destroy()<br>file.close()<br><br><br><div class="gmail_quote">On Mon, Mar 1, 2010 at 12:54 PM, Jamie Lahowetz <span dir="ltr"><<a href="mailto:deadpickle@gmail.com">deadpickle@gmail.com</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin: 0pt 0pt 0pt 0.8ex; border-left: 1px solid rgb(204, 204, 204); padding-left: 1ex;"><div><div></div><div class="h5">I'm new to GDAL and am usingthe python bindings. I want to transform a
newly created shapefile from WGS84 to a Lambert Conical Conformal so
that I can calculate the length of the contained polylines. The script
exits with no errors but when I look at the shapefile in arcCatalog I
see that the projection is still WGS84. Any ideas? Also, is there a way
to calculate the length directly from GDAL? Thanks.<br>
<br>import os<br>try:<br> from osgeo import ogr<br> from osgeo
import osr<br><br>except:<br> import ogr<br> import osr<br><br>#~
start driver<br>driver = ogr.GetDriverByName('ESRI Shapefile')<br>#~ set
workspace<br>
ws = 'C:\\Users\\deadpickle\\<div>Desktop\\case study\\verify\\'<br>os.chdir(ws)<br><br>#~
check if it exists<br>if os.path.exists('test.shp'):<br>
driver.DeleteDataSource('test.shp')<br><br>#~ create a new shp<br>
ds = driver.CreateDataSource('test.shp')<br><br>#~ create spatref
for new shp<br>source = osr.SpatialReference()<br>source.SetFromUserInput('WGS84')<br><br>#~
create a new data layer<br>layer = ds.CreateLayer('test',geom_type=ogr.wkbLineString,srs=source)<br>
<br>#~ add an id field<br>fieldDefn = ogr.FieldDefn('id',
ogr.OFTInteger)<br>layer.CreateField(fieldDefn)<br><br>#~ create a
geometry<br>line = ogr.Geometry(ogr.wkbLineString)<br><br>#~ get
the def of the layer, scheme<br>
featureDefn = layer.GetLayerDefn()<br><br># create a new feature<br>feature
= ogr.Feature(featureDefn)<br><br>#~ Read 2006track file<br>file =
open("2006track.csv","r")<br>file.readline()<br>lnlist =
file.readlines()<br>
for lines in lnlist:<br> value = lines.split(",")<br> if
float(value[1]) == 2:<br> print value<br>
line.AddPoint(float(value[3]),float(value[2]))<br> #~ set
the field 'ID'<br> feature.SetField('id', float(value[1]))<br>
<br>#~ set the geometry for this shp<br>feature.SetGeometry(line)<br><br>#
add the feature to the output layer<br>layer.CreateFeature(feature)<br><br>#~
feature = layer.GetFeature(0)<br>geom = feature.GetGeometryRef()<br>
<br>#~ change projection to lcc<br>target = osr.SpatialReference()<br>target.ImportFromProj4('+proj=lcc
+lat_1= 33.000000 +lat_2=45.000000 +lat_0=39.000000 +lon_0=-96.000000
+x_0=0.0 +y_0=0.0')<br>coordtrans = osr.CoordinateTransformation(source,
target)<br>
<br>geom.Transform(coordtrans)<br><br># destroy the geometry and feature
and close the data source<br><br>ds.Destroy()<br>line.Destroy()<br>feature.Destroy()<br>file.close()</div>
</div></div></blockquote></div><br><br clear="all"><br>-- <br>Jamie Ryan Lahowetz<br>University of Nebraska - Lincoln<br>Graduate Student - Geosciences<br>402.304.0766<br><a href="mailto:jlahowetz@gmail.com">jlahowetz@gmail.com</a><br>