[gdal-dev] Re: Transform not working

[Frank Warmerdam]

deadpickle 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.

Jamie,

My suggestion in the following script would be to create the file with
the "target" srs instead of WGS84, and to transform the geometry
before even assigning it to the feature that will be written.

A few notes:
  * There is no way with OGR of changing the coordinate system of
    an existing layer.  The only point at which it can be set is
    layer creation time.
  * The CreateFeature() method on the layer creates the feature in
    the datastore (ie. shapefile) but there is no meaningful
    connection maintained between the ogr.Feature object and the
    feature in the file after this.  Alterations to the feature
    or it's geometry will have no effect on the file without doing
    another layer.SetFeature() call.

>> 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
> 


-- 
---------------------------------------+--------------------------------------
I set the clouds in motion - turn up   | Frank Warmerdam, warmerdam at pobox.com
light and sound - activate the windows | http://pobox.com/~warmerdam
and watch the world go round - Rush    | Geospatial Programmer for Rent