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