Thanks for the response. I'm sorry I'm really new to this and am not sure where to transform the line geometry, I know it should be before layer.CreateFeature(feature) but after I add the points to the geometry. Is it before or after feature.SetGeometry(line)? And is the right way to call the transform <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 +datum=NAD83')<br>coordtrans = osr.CoordinateTransformation(source, target)<br>
line.Transform(coordtrans)?<br><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>source.SetFromUserInput('WGS84')<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>
#~ 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>#~ 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 +datum=NAD83')<br>
coordtrans = osr.CoordinateTransformation(source, target)<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><br><br><div class="gmail_quote">On Mon, Mar 1, 2010 at 6:49 PM, Chris Garrard <span dir="ltr"><<a href="mailto:garrard.chris%2Bgis@gmail.com">garrard.chris+gis@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;">Jamie, Try transforming the line geometry before adding it to the feature. That should do it.<br>
<font color="#888888"><br>chris</font><div><div></div><div class="h5"><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>
</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>