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 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 +datum=NAD83') coordtrans = osr.CoordinateTransformation(source, target) line.Transform(coordtrans)?
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') #~ 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') #~ 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])) #~ 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 +datum=NAD83') coordtrans = osr.CoordinateTransformation(source, target) 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 6:49 PM, Chris Garrard <garrard.chris+...@gmail.com<garrard.chris%2b...@gmail.com> > wrote: > Jamie, Try transforming the line geometry before adding it to the feature. > That should do it. > > chris > > > On Mon, Mar 1, 2010 at 4:24 PM, deadpickle <deadpic...@gmail.com> wrote: > >> Not much of a response going on here. In order to change the shapefile >> projection from geographic to projected do I need to reopen the file, >> get each feature, and transform the points individually? If so what >> calls do I need to invoke? >> >> On Mar 1, 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() >> > >> > _______________________________________________ >> > gdal-dev mailing list >> > gdal-...@lists.osgeo.orghttp:// >> lists.osgeo.org/mailman/listinfo/gdal-dev >> _______________________________________________ >> gdal-dev mailing list >> gdal-dev@lists.osgeo.org >> http://lists.osgeo.org/mailman/listinfo/gdal-dev >> > > -- 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