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.
import os try: from osgeo import ogr from osgeo import osr except: import ogr import osr #~ start driver driver = ogr.GetDriverByName('ESRI Shapefile') #~ set workspace ws = 'C:\\Users\\deadpickle\\Desktop\\case study\\verify\\' os.chdir(ws) #~ check if it exists if os.path.exists('test.shp'): driver.DeleteDataSource('test.shp') #~ create a new shp ds = driver.CreateDataSource('test.shp') #~ create spatref for new shp source = osr.SpatialReference() target = osr.SpatialReference() 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') target.SetFromUserInput('WGS84') coordtrans = osr.CoordinateTransformation(target, source) #~ create a new data layer layer = ds.CreateLayer('test',geom_type=ogr.wkbLineString,srs=source) #~ add an id field fieldDefn = ogr.FieldDefn('id', ogr.OFTInteger) layer.CreateField(fieldDefn) #~ create a geometry line = ogr.Geometry(ogr.wkbLineString) #~ get the def of the layer, scheme featureDefn = layer.GetLayerDefn() # create a new feature feature = ogr.Feature(featureDefn) #~ Read 2006track file file = open("2006track.csv","r") file.readline() lnlist = file.readlines() for lines in lnlist: value = lines.split(",") if float(value[1]) == 2: print float(value[3]),float(value[2]) line.AddPoint(float(value[3]),float(value[2])) #~ set the field 'ID' feature.SetField('id', float(value[1])) #~ trans from lat/lon input to lcc line.Transform(coordtrans) #~ set the geometry for this shp feature.SetGeometry(line) # add the feature to the output layer layer.CreateFeature(feature) # destroy the geometry and feature and close the data source ds.Destroy() line.Destroy() feature.Destroy() file.close() On Mon, Mar 1, 2010 at 12:54 PM, Jamie Lahowetz <deadpic...@gmail.com>wrote: > 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. > > import os > try: > from osgeo import ogr > from osgeo import osr > > except: > import ogr > import osr > > #~ start driver > driver = ogr.GetDriverByName('ESRI Shapefile') > #~ set workspace > ws = 'C:\\Users\\deadpickle\\ > Desktop\\case study\\verify\\' > os.chdir(ws) > > #~ check if it exists > if os.path.exists('test.shp'): > driver.DeleteDataSource('test.shp') > > #~ create a new shp > ds = driver.CreateDataSource('test.shp') > > #~ create spatref for new shp > source = osr.SpatialReference() > source.SetFromUserInput('WGS84') > > #~ create a new data layer > layer = ds.CreateLayer('test',geom_type=ogr.wkbLineString,srs=source) > > #~ add an id field > fieldDefn = ogr.FieldDefn('id', ogr.OFTInteger) > layer.CreateField(fieldDefn) > > #~ create a geometry > line = ogr.Geometry(ogr.wkbLineString) > > #~ get the def of the layer, scheme > featureDefn = layer.GetLayerDefn() > > # create a new feature > feature = ogr.Feature(featureDefn) > > #~ Read 2006track file > file = open("2006track.csv","r") > file.readline() > lnlist = file.readlines() > for lines in lnlist: > value = lines.split(",") > if float(value[1]) == 2: > print value > line.AddPoint(float(value[3]),float(value[2])) > #~ set the field 'ID' > feature.SetField('id', float(value[1])) > > #~ set the geometry for this shp > feature.SetGeometry(line) > > # add the feature to the output layer > layer.CreateFeature(feature) > > #~ feature = layer.GetFeature(0) > geom = feature.GetGeometryRef() > > #~ change projection to lcc > target = osr.SpatialReference() > 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') > coordtrans = osr.CoordinateTransformation(source, target) > > geom.Transform(coordtrans) > > # destroy the geometry and feature and close the data source > > ds.Destroy() > line.Destroy() > feature.Destroy() > file.close() > -- Jamie Ryan Lahowetz University of Nebraska - Lincoln Graduate Student - Geosciences 402.304.0766 jlahow...@gmail.com
_______________________________________________ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev