[gdal-dev] Re: Transform not working

[Chris Garrard]

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 <deadpickle at 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... at 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-... at lists.osgeo.orghttp://lists.osgeo.org/mailman/listinfo/gdal-dev
> _______________________________________________
> gdal-dev mailing list
> gdal-dev at lists.osgeo.org
> http://lists.osgeo.org/mailman/listinfo/gdal-dev
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://lists.osgeo.org/pipermail/gdal-dev/attachments/20100301/458f897d/attachment-0001.html