Jamie, Try transforming the line geometry before adding it to the feature. That should do it.<br><br>chris<br><br><div class="gmail_quote">On Mon, Mar 1, 2010 at 4:24 PM, deadpickle <span dir="ltr"><<a href="mailto:deadpickle@gmail.com" target="_blank">deadpickle@gmail.com</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="border-left: 1px solid rgb(204, 204, 204); margin: 0pt 0pt 0pt 0.8ex; padding-left: 1ex;">Not much of a response going on here. In order to change the shapefile<br>
projection from geographic to projected do I need to reopen the file,<br>
get each feature, and transform the points individually? If so what<br>
calls do I need to invoke?<br>
<br>
On Mar 1, 12:54 pm, Jamie Lahowetz <<a href="mailto:deadpic...@gmail.com" target="_blank">deadpic...@gmail.com</a>> wrote:<br>
> I'm new to GDAL and am usingthe python bindings. I want to transform a newly<br>
> created shapefile from WGS84 to a Lambert Conical Conformal so that I can<br>
> calculate the length of the contained polylines. The script exits with no<br>
> errors but when I look at the shapefile in arcCatalog I see that the<br>
> projection is still WGS84. Any ideas? Also, is there a way to calculate the<br>
> 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\\<br>
> 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<br>
> +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()<br>
><br>
> _______________________________________________<br>
> gdal-dev mailing list<br>
> gdal-...@lists.osgeo.orghttp://<a href="http://lists.osgeo.org/mailman/listinfo/gdal-dev" target="_blank">lists.osgeo.org/mailman/listinfo/gdal-dev</a><br>
_______________________________________________<br>
gdal-dev mailing list<br>
<a href="mailto:gdal-dev@lists.osgeo.org" target="_blank">gdal-dev@lists.osgeo.org</a><br>
<a href="http://lists.osgeo.org/mailman/listinfo/gdal-dev" target="_blank">http://lists.osgeo.org/mailman/listinfo/gdal-dev</a><br>
</blockquote></div><br>