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 >
_______________________________________________ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev